Restarting doesn't work but non-restarted Davidson seems to produce a…

…ccurate results (sort of). Saving for now.
primme · Oct 25, 2023 · 775a1d8 · 775a1d8
1 parent cbcaee1
commit 775a1d8
Show file tree

Hide file tree

Showing 8 changed files with 1,297 additions and 310 deletions.
diff --git a/examples/ex_eigs_dparTests.c b/examples/ex_eigs_dparTests.c
diff --git a/include/primme_eigs.h b/include/primme_eigs.h
@@ -61,7 +61,7 @@ typedef enum {
    primme_expansion_davidson
 } primme_expansion;
 
-/* TODO: Remove after I'm done testing sketched Lanczos */
+/* Used for testing sketched Lanczos */
 typedef enum {
    primme_residual_sketched,
    primme_residual_RQ,
@@ -176,6 +176,11 @@ typedef struct residual_params {
    primme_residual residual;
 } residual_params;
 
+typedef struct sketching_params {
+   PRIMME_INT sketchSize;
+   PRIMME_INT nnzPerCol;
+} sketching_params;
+
 typedef struct correction_params {
    int precondition;
    int robustShifts;
@@ -266,6 +271,7 @@ typedef struct primme_params {
    struct projection_params projectionParams; 
    struct expansion_params expansionParams; 
    struct residual_params residualParams; 
+   struct sketching_params sketchingParams; 
    struct restarting_params restartingParams;
    struct correction_params correctionParams;
    struct primme_stats stats;

diff --git a/src/eigs/lanczos.c b/src/eigs/lanczos.c
@@ -222,7 +222,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
    SCALAR *SW;                /* The projected sketched basis                  */
    SCALAR *S_vals;                                       /* CSC Formatted Values */
    PRIMME_INT *S_rows;                                          /* CSC Formatted Rows */
-   PRIMME_INT nnz_per_col, ldSV, ldSW;  /* Size and nnz of the sketching matrix */
+   PRIMME_INT nnzPerCol, ldSV, ldSW;  /* Size and nnz of the sketching matrix */
    REAL *normalize_evecs;
 
    /* -------------------------------------------------------------- */
@@ -327,11 +327,11 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
    if(primme->projectionParams.projection == primme_proj_sketched)
    {
       /* Default settings for sketch size and nnz per column. Based on Yuji and Tropp's manuscript */
-      ldSV = ldSW = 4*maxBasisSize;
-      nnz_per_col = (int)(ceil(2*log(maxBasisSize+1)));   
+      ldSV = ldSW = primme->sketchingParams.sketchSize;
+      nnzPerCol = primme->sketchingParams.nnzPerCol; 
 
-      CHKERR(Num_malloc_Sprimme(nnz_per_col*nLocal, &S_vals, ctx));
-      S_rows = (PRIMME_INT*)malloc(nnz_per_col*nLocal*sizeof(PRIMME_INT));
+      CHKERR(Num_malloc_Sprimme(nnzPerCol*nLocal, &S_vals, ctx));
+      S_rows = (PRIMME_INT*)malloc(nnzPerCol*nLocal*sizeof(PRIMME_INT));
 
       CHKERR(Num_malloc_Sprimme(ldSV*(maxBasisSize+maxBlockSize), &SV, ctx));
       CHKERR(Num_malloc_Sprimme(ldSW*maxBasisSize, &SW, ctx));
@@ -340,7 +340,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
       CHKERR(Num_malloc_Sprimme(ldV*(maxBasisSize+maxBlockSize), &V_temp, ctx));
 
       /* Build Sketch CSR Locally */
-      CHKERR(build_sketch_Sprimme(S_rows, S_vals, ldSV, nnz_per_col, ctx));
+      CHKERR(build_sketch_Sprimme(S_rows, S_vals, ctx));
 
    } /* End sketching matrix build */
 
@@ -354,7 +354,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
 
    /* Update sketched basis if sketching is turned on */
    if(primme->projectionParams.projection == primme_proj_sketched)
-      CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, 0, blockSize, nnz_per_col, S_rows, S_vals, ctx));
+      CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, 0, blockSize, S_rows, S_vals, ctx));
 
    t0 = primme_wTimer(); /* Start timing our basis build */
 
@@ -383,7 +383,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
 
       /* Update sketched basis if sketching is turned on */
       if(primme->projectionParams.projection == primme_proj_sketched)
-         CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, i, blockSize, nnz_per_col, S_rows, S_vals, ctx));
+         CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, i, blockSize, S_rows, S_vals, ctx));
 
       if(primme->printLevel == 4 && (PRIMME_INT)(i+blockSize) % 10 == 0)
       {
@@ -400,7 +400,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
             CHKERR(ortho_Sprimme(&V_temp[ldV*(i+blockSize)], ldV, &H[i*ldH + (i+blockSize)], ldH, 0, blockSize-1, NULL, 0, 0, nLocal, primme->iseed, ctx));   /* [V_i, b_i] = qr(V_i) */
 
             /* SW = SV*H */
-            CHKERR(sketch_basis_Sprimme(V_temp, ldV, SV, ldSV, i+blockSize, blockSize, nnz_per_col, S_rows, S_vals, ctx));
+            CHKERR(sketch_basis_Sprimme(V_temp, ldV, SV, ldSV, i+blockSize, blockSize, S_rows, S_vals, ctx));
             CHKERR(Num_gemm_Sprimme("N", "N", ldSV, i+blockSize, i+2*blockSize, 1.0, SV, ldSV, H, ldH, 0.0, SW, ldSW, ctx));
 
             /* Getting our sketched basis and projected sketched basis */
@@ -497,7 +497,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
          CHKERR(ortho_Sprimme(V, ldV, NULL, 0, maxBasisSize, maxBasisSize+blockSize-1, NULL, 0, 0, nLocal, primme->iseed, ctx));   /* Orthogonalized the last block of V against the rest of the basis */
 
          /* SW = SV*H */
-         CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, maxBasisSize, blockSize, nnz_per_col, S_rows, S_vals, ctx));
+         CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, maxBasisSize, blockSize, S_rows, S_vals, ctx));
          CHKERR(Num_gemm_Sprimme("N", "N", ldSV, maxBasisSize, maxBasisSize+blockSize, 1.0, SV, ldSV, H, ldH, 0.0, SW, ldSW, ctx));
 
          /* Getting our sketched basis and projected sketched basis */