Skip to content

Commit

Permalink
Debugged Davidson - One *minor* issue with verification causing a nev…
Browse files Browse the repository at this point in the history 
…erending loop
  • Loading branch information
@Heatherms27
Heatherms27 committed Mar 27, 2024
1 parent 5045c44 commit f8c144e
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 30 deletions.
33 changes: 20 additions & 13 deletions src/eigs/main_iter.c
View file
Expand Up @@ -1593,16 +1593,16 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

/* Variables for sketching */
SCALAR *SV; /* Sketched Basis vectors */
SCALAR *Q; /* The "Q" factor in the QR decomposition of SV */
SCALAR *T; /* The "R" factor in the QR decomposition of SV */
SCALAR *Q = NULL; /* The "Q" factor in the QR decomposition of SV */
SCALAR *T = NULL; /* The "R" factor in the QR decomposition of SV */
SCALAR *SW; /* Projected sketched Basis vectors */
SCALAR *V_temp; /* Work array for sketching */
SCALAR *S_vals; /* Holds the nonzero entries of the sketching matrix */
REAL *normalize_evecs; /* For normalizing the sketched Ritz vectors */
PRIMME_INT *S_rows; /* Holds the row numbers of the nonzero entries */
PRIMME_INT nnzPerCol; /* NNZ per column in the sketching matrix */
PRIMME_INT ldSV, ldQ, ldT, ldSW; /* Leading dimensions of SV, T, and SW */

ldQ = ldT = 0;

/* -------------------------------------------------------------- */
/* Allocate objects */
Expand Down Expand Up @@ -1654,23 +1654,28 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(Num_malloc_iprimme(maxEvecsSize, &ipivot, ctx));

/* Allocate space for the variables needed for sketching (if needed) */
ldSV = ldSW = ldQ = primme->sketchingParams.sketchSize;
ldT = primme->maxBasisSize;
ldSV = ldSW = primme->sketchingParams.sketchSize;
nnzPerCol = primme->sketchingParams.nnzPerCol;

S_rows = (PRIMME_INT*)malloc(nnzPerCol*primme->nLocal*sizeof(PRIMME_INT));

CHKERR(Num_malloc_Sprimme(primme->nLocal*nnzPerCol, &S_vals, ctx));
CHKERR(Num_malloc_Sprimme(ldSV*primme->maxBasisSize, &SV, ctx));
CHKERR(Num_malloc_Sprimme(ldQ*primme->maxBasisSize, &Q, ctx));
CHKERR(Num_malloc_Sprimme(ldT*ldT, &T, ctx));
CHKERR(Num_malloc_Sprimme(ldSW*primme->maxBasisSize, &SW, ctx));
CHKERR(Num_malloc_Sprimme(ldV*primme->maxBasisSize, &V_temp, ctx));

CHKERR(Num_malloc_Rprimme(primme->numEvals, &normalize_evecs, ctx));
CHKERR(Num_zero_matrix_Sprimme(SV, ldSV, primme->maxBasisSize, ldSV, ctx));
CHKERR(Num_zero_matrix_Sprimme(T, ldT, ldT, ldT, ctx));

if(primme->procID == 0) {
ldQ = primme->sketchingParams.sketchSize;
ldT = primme->maxBasisSize;

CHKERR(Num_malloc_Sprimme(ldQ*primme->maxBasisSize, &Q, ctx));
CHKERR(Num_malloc_Sprimme(ldT*ldT, &T, ctx));

CHKERR(Num_zero_matrix_Sprimme(T, ldT, ldT, ldT, ctx));
}


/* -------------------------------------------------------------- */
Expand Down Expand Up @@ -2126,7 +2131,7 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

assert(ldV == ldW); /* this function assumes ldV == ldW */

CHKERR(restart_sketched(V, ldV, W, ldW, SV, ldSV, Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, min(basisSize, primme->minRestartSize), &basisSize, S_rows, S_vals, ctx));
CHKERR(restart_sketched(V, ldV, W, ldW, SV, ldSV, Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, min(basisSize, primme->minRestartSize), &basisSize, ctx));

restartsSinceReset++;

Expand Down Expand Up @@ -2330,12 +2335,15 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(Num_free_Rprimme(normalize_evecs, ctx));

CHKERR(Num_free_Sprimme(SV, ctx));
CHKERR(Num_free_Sprimme(Q, ctx));
CHKERR(Num_free_Sprimme(T, ctx));
CHKERR(Num_free_Sprimme(SW, ctx));
CHKERR(Num_free_Sprimme(S_vals, ctx));
CHKERR(Num_free_Sprimme(V_temp, ctx));

if(primme->procID == 0){
CHKERR(Num_free_Sprimme(Q, ctx));
CHKERR(Num_free_Sprimme(T, ctx));
}

return 0;
}

Expand Down Expand Up @@ -2809,8 +2817,7 @@ STATIC int verify_norms(SCALAR *V, PRIMME_INT ldV, SCALAR *W, PRIMME_INT ldW,
}

CHKERR(globalSum_RHprimme(resNorms, basisSize, ctx));
for (i=0; i < basisSize; i++)
resNorms[i] = sqrt(resNorms[i]);
for (i=0; i < basisSize; i++) resNorms[i] = sqrt(resNorms[i]);

/* Check for convergence of the residual norms. */

Expand Down
38 changes: 22 additions & 16 deletions src/eigs/restart.c
View file
Expand Up @@ -1774,22 +1774,19 @@ STATIC int restart_RR(HSCALAR *H, int ldH, HSCALAR *VtBV, int ldVtBV,

STATIC int restart_sketched(SCALAR *V, int ldV, SCALAR *W, int ldW, SCALAR *SV,
int ldSV, SCALAR *Q, PRIMME_INT ldQ, SCALAR *T, int ldT, SCALAR *SW, int ldSW, HSCALAR *hVecs,
int ldhVecs, HEVAL *hVals, int restartSize, int *basisSize,
PRIMME_INT *S_rows, SCALAR *S_vals, primme_context ctx) {
int ldhVecs, HEVAL *hVals, int restartSize, int *basisSize, primme_context ctx) {

primme_params *primme = ctx.primme;

int i; /* Loop variable */
SCALAR *V_temp;
SCALAR *SV_temp;
SCALAR *VecNorms;
SCALAR *q;

int old_basisSize = *basisSize;
CHKERR(Num_malloc_Sprimme(ldV*restartSize, &V_temp, ctx));
CHKERR(Num_malloc_Sprimme(ldSV*restartSize, &SV_temp, ctx));
CHKERR(Num_malloc_Sprimme(restartSize, &VecNorms, ctx));
CHKERR(Num_malloc_Sprimme(old_basisSize*restartSize, &q, ctx));

//CHKERR(Num_malloc_Sprimme(old_basisSize*restartSize, &hVecs_temp, ctx));

Expand Down Expand Up @@ -1827,9 +1824,6 @@ STATIC int restart_sketched(SCALAR *V, int ldV, SCALAR *W, int ldW, SCALAR *SV,
CHKERR(Num_copy_matrix_Sprimme(SV_temp, ldSV, restartSize, ldSW, SW, ldSW, ctx)); // Copy the temporary matrix back into SW
CHKERR(Num_free_Sprimme(SV_temp, ctx));

// q = T * hVecs
CHKERR(Num_gemm_Sprimme("N", "N", old_basisSize, restartSize, old_basisSize, 1.0, T, ldT, hVecs, ldhVecs, 0.0, q, old_basisSize, ctx));

// Normalize the matrices
for(i = 0; i < restartSize; i++) VecNorms[i] = sqrt(Num_dot_Sprimme(primme->nLocal, &V[i*ldV], 1, &V[i*ldV], 1, ctx));
CHKERR(globalSum_Sprimme(VecNorms, restartSize, ctx));
Expand All @@ -1839,19 +1833,31 @@ STATIC int restart_sketched(SCALAR *V, int ldV, SCALAR *W, int ldW, SCALAR *SV,
CHKERR(Num_scal_Sprimme(ldW, 1.0/VecNorms[i], &W[i*ldW], 1, ctx));
CHKERR(Num_scal_Sprimme(ldSV, 1.0/VecNorms[i], &SV[i*ldSV], 1, ctx));
CHKERR(Num_scal_Sprimme(ldSW, 1.0/VecNorms[i], &SW[i*ldSW], 1, ctx));
CHKERR(Num_scal_Sprimme(old_basisSize, 1.0/VecNorms[i], &q[i*old_basisSize], 1, ctx));
}

// Recompute Q and R factors
CHKERR(ortho_Sprimme(q, old_basisSize, T, ldT, 0, restartSize-1, NULL, 0, 0, old_basisSize, primme->iseed, ctx)); // [q, r] = qr(T*hVecs*Vn)
SCALAR *Q_temp;
CHKERR(Num_malloc_Sprimme(ldQ*old_basisSize, &Q_temp, ctx));
CHKERR(Num_copy_matrix_Sprimme(Q, ldQ, old_basisSize, ldQ, Q_temp, ldQ, ctx)); // Temporary matrix to not overwrite Q
CHKERR(Num_gemm_SHprimme("N", "N", ldQ, restartSize, old_basisSize, 1.0, Q_temp, ldQ, q, old_basisSize, 0.0, Q, ldQ, ctx)); // Q = Q*q
// Recompute Q and R factors of the sketched basis, but only for process 0
if(primme->procID == 0){
SCALAR *q;
SCALAR *Q_temp;

CHKERR(Num_malloc_Sprimme(old_basisSize*restartSize, &q, ctx));
CHKERR(Num_malloc_Sprimme(ldQ*old_basisSize, &Q_temp, ctx));

CHKERR(Num_copy_matrix_Sprimme(Q, ldQ, old_basisSize, ldQ, Q_temp, ldQ, ctx)); // Temporary matrix to not overwrite Q

// q = T * hVecs * Vn
CHKERR(Num_gemm_Sprimme("N", "N", old_basisSize, restartSize, old_basisSize, 1.0, T, ldT, hVecs, ldhVecs, 0.0, q, old_basisSize, ctx)); // bar(q) = T * hVecs
for (i = 0; i < restartSize; i++) CHKERR(Num_scal_Sprimme(old_basisSize, 1.0/VecNorms[i], &q[i*old_basisSize], 1, ctx)); // bar(q) *= Vn

// [q, r] = qr(bar(q))
CHKERR(Bortho_local_Sprimme(q, old_basisSize, T, ldT, 0, restartSize-1, NULL, 0, 0, old_basisSize, NULL, 0, primme->iseed, ctx));
CHKERR(Num_gemm_Sprimme("N", "N", ldQ, restartSize, old_basisSize, 1.0, Q_temp, ldQ, q, old_basisSize, 0.0, Q, ldQ, ctx)); // Q = Q*q

CHKERR(Num_free_Sprimme(Q_temp, ctx));
CHKERR(Num_free_Sprimme(q, ctx));
}

CHKERR(Num_free_Sprimme(VecNorms, ctx));
CHKERR(Num_free_Sprimme(Q_temp, ctx));
CHKERR(Num_free_Sprimme(q, ctx));

// Update eiganpairs and residuals
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, restartSize, hVals, restartSize, ctx));
Expand Down
1 change: 0 additions & 1 deletion src/eigs/sketch.c
View file
Expand Up @@ -166,7 +166,6 @@ int sketch_basis_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *SV, PRIMME_INT ldSV,
if(T) {
CHKERR(Num_copy_matrix_Sprimme(&SV[basisSize*ldSV], ldSV, blockSize, ldSV, &Q[ldQ*basisSize], ldQ, ctx));
CHKERR(Num_zero_matrix_Sprimme(&T[basisSize], blockSize, basisSize, ldT, ctx));
//CHKERR(ortho_Sprimme(Q, ldQ, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, ldQ, primme->iseed, ctx));
CHKERR(Bortho_local_Sprimme(Q, ldQ, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, ldQ, NULL, 0, primme->iseed, ctx));
}

Expand Down

0 comments on commit f8c144e

Please sign in to comment.