From f8c144e63653a5215742de76981d3b037236b399 Mon Sep 17 00:00:00 2001 From: Heather M Switzer Date: Wed, 27 Mar 2024 18:58:40 -0400 Subject: [PATCH] Debugged Davidson - One *minor* issue with verification causing a neverending loop --- src/eigs/main_iter.c | 33 ++++++++++++++++++++------------- src/eigs/restart.c | 38 ++++++++++++++++++++++---------------- src/eigs/sketch.c | 1 - 3 files changed, 42 insertions(+), 30 deletions(-) diff --git a/src/eigs/main_iter.c b/src/eigs/main_iter.c index 30a82d8b..314ef1cb 100644 --- a/src/eigs/main_iter.c +++ b/src/eigs/main_iter.c @@ -1593,8 +1593,8 @@ 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 */ @@ -1602,7 +1602,7 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs, 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 */ @@ -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)); + } /* -------------------------------------------------------------- */ @@ -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++; @@ -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; } @@ -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. */ diff --git a/src/eigs/restart.c b/src/eigs/restart.c index e1439b90..dc0d52ad 100644 --- a/src/eigs/restart.c +++ b/src/eigs/restart.c @@ -1774,8 +1774,7 @@ 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; @@ -1783,13 +1782,11 @@ STATIC int restart_sketched(SCALAR *V, int ldV, SCALAR *W, int ldW, SCALAR *SV, 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)); @@ -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)); @@ -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)); diff --git a/src/eigs/sketch.c b/src/eigs/sketch.c index bf686e0d..68798bdd 100644 --- a/src/eigs/sketch.c +++ b/src/eigs/sketch.c @@ -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)); }