Skip to content

Commit

Permalink
Fixed the QR issue, but not the condition number issue
Browse files Browse the repository at this point in the history
  • Loading branch information
@Heatherms27
Heatherms27 committed Mar 24, 2024
1 parent 3c0aaf3 commit 2d2e1c1
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 26 deletions.
2 changes: 1 addition & 1 deletion src/eigs/lanczos.c
View file
Expand Up @@ -260,7 +260,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

CHKERR(Num_malloc_Rprimme(primme->numEvals, &normalize_evecs, ctx));

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

/* Build Sketch CSR Locally */
CHKERR(build_sketch_Sprimme(S_rows, S_vals, ctx));
Expand Down
7 changes: 4 additions & 3 deletions src/eigs/main_iter.c
View file
Expand Up @@ -1574,7 +1574,6 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
HSCALAR *hVecs; /* Eigenvectors of H */
HSCALAR *prevhVecs=NULL; /* hVecs from previous iteration */

PRIMME_INT ldQ; /* The leading dimension of Q */
HSCALAR *hVecsRot = NULL; /* transformation of hVecs in arbitrary vectors */

HEVAL *hVals; /* Eigenvalues of H */
Expand Down Expand Up @@ -1610,9 +1609,9 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

maxEvecsSize = primme->numOrthoConst + primme->numEvals;

/* Use leading dimension ldOPs for the large dimension mats: V, W and Q */
/* Use leading dimension ldOPs for the large dimension mats: V, W */

ldV = ldW = ldBV = ldQ = primme->ldOPs;
ldV = ldW = ldBV = primme->ldOPs;
ldBevecs = primme->massMatrixMatvec ? primme->ldOPs : ldevecs;
if (primme->massMatrixMatvec) {
CHKERR(Num_malloc_Sprimme(ldBV*primme->maxBasisSize, &BV, ctx));
Expand Down Expand Up @@ -1667,6 +1666,8 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(Num_malloc_Sprimme(ldV*primme->maxBasisSize, &V_temp, ctx));

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



/* -------------------------------------------------------------- */
Expand Down
45 changes: 23 additions & 22 deletions src/eigs/sketch.c
View file
Expand Up @@ -162,7 +162,8 @@ int sketch_basis_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *SV, PRIMME_INT ldSV,

/* Find the sketched basis */
CHKERR(globalSum_Sprimme(&SV[basisSize*ldSV], blockSize*ldSV, ctx));
if (T) CHKERR(ortho_Sprimme(SV, ldSV, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, primme->maxBasisSize, primme->iseed, ctx));
CHKERR(Num_zero_matrix_Sprimme(&T[basisSize], blockSize, basisSize, ldT, ctx));
CHKERR(ortho_Sprimme(SV, ldSV, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, ldSV, primme->iseed, ctx));

if (primme) primme->stats.timeSketchMatvec += primme_wTimer() - t0;

Expand Down Expand Up @@ -274,28 +275,16 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT,
CHKERR(Num_malloc_iprimme(basisSize, &eval_perm, ctx));

/* Check condition number to see if stabilization is needed */
CHKERR(Num_copy_matrix_Sprimme(T, basisSize, basisSize, ldT, T_temp, basisSize, ctx)); /* Ensure T is not overwritten */
CHKERR(Num_gesvd_Sprimme("N", "N", basisSize, basisSize, T_temp, basisSize, sing_vals, NULL, basisSize, NULL, basisSize, ctx));
printf("Checking condition number\n");
CHKERR(Num_copy_matrix_Sprimme(&T[0], basisSize, basisSize, ldT, &T_temp[0], basisSize, ctx)); /* Ensure T is not overwritten */
CHKERR(Num_trcon_Sprimme("O", "U", "N", basisSize, T_temp, basisSize, &cond_est, ctx));

cond_est = sing_vals[0]/sing_vals[basisSize-1];
primme->stats.estimateLargestSVal = sing_vals[0];
printf("Estimated Cond Num: %E\n", cond_est);

/* XXX: Stabilization needed */
if(cond_est > 1/MACHINE_EPSILON) {
double stab_time = primme_wTimer(); // For monitoring

/* Truncate the basis for stabilization */
for(i = basisSize-1; i > primme->numEvals; i--)
{
if(sing_vals[0]/sing_vals[i] >= 1/MACHINE_EPSILON)
{
tbasisSize--;
} else {
break;
}
}

SCALAR *UtSWV; /* Left hand side of the generalized eigenvalue problem */
SCALAR *UVecs; /* Left singular vectors of the SVD */
SCALAR *VVecst; /* Right singular vectors of the SVD (transposed) */
Expand All @@ -306,20 +295,32 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT,

ldUVecs = ldSV;
ldVVecst = basisSize;
ldSigma = ldUtSW = ldUtSWV = ldtrunc_hVecs = tbasisSize;

CHKERR(Num_malloc_Sprimme(ldUVecs*basisSize, &UVecs, ctx));
CHKERR(Num_malloc_Sprimme(ldVVecst*basisSize, &VVecst, ctx));
CHKERR(Num_malloc_Sprimme(ldSV*basisSize, &temp_SV, ctx));
CHKERR(Num_malloc_Sprimme(ldUtSWV*tbasisSize, &UtSWV, ctx));
CHKERR(Num_malloc_Sprimme(ldtrunc_hVecs*tbasisSize, &trunc_hVecs, ctx));
CHKERR(Num_malloc_Sprimme(ldSigma*tbasisSize, &Sigma, ctx));

/* Take the SVD decomposition of SV */
CHKERR(Num_copy_matrix_Sprimme(&SV[0], ldSV, basisSize, ldSV, &temp_SV[0], ldSV, ctx));
CHKERR(Num_gesvd_Sprimme("S", "S", ldSV, basisSize, temp_SV, ldSV, sing_vals, UVecs, ldUVecs, VVecst, ldVVecst, ctx));
CHKERR(Num_free_Sprimme(temp_SV, ctx));

/* Truncate the basis for stabilization */
for(i = basisSize-1; i > primme->numEvals; i--)
{
if(sing_vals[0]/sing_vals[i] >= 1/MACHINE_EPSILON)
{
tbasisSize--;
} else {
break;
}
}

ldSigma = ldUtSW = ldUtSWV = ldtrunc_hVecs = tbasisSize;
CHKERR(Num_malloc_Sprimme(ldUtSWV*tbasisSize, &UtSWV, ctx));
CHKERR(Num_malloc_Sprimme(ldtrunc_hVecs*tbasisSize, &trunc_hVecs, ctx));
CHKERR(Num_malloc_Sprimme(ldSigma*tbasisSize, &Sigma, ctx));

/* Construct Mass Matrix (Diagonal matrix with singular values as entries) */
CHKERR(Num_zero_matrix_Sprimme(Sigma, tbasisSize, tbasisSize, ldSigma, ctx));
for(i = 0; i < tbasisSize; i++) Sigma[i*ldSigma+i] = sing_vals[i];
Expand Down Expand Up @@ -347,8 +348,8 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT,
CHKERR(Num_gemm_Sprimme("C", "N", basisSize, basisSize, ldSV, 1.0, SV, ldSV, SW, ldSW, 0.0, UtSW, basisSize, ctx));

/* Solve the eigenproblem */
CHKERR(Num_copy_matrix_Sprimme(T, basisSize, basisSize, ldT, T_temp, basisSize, ctx)); /* Ensure T is not overwritten */
CHKERR(Num_ggev_Sprimme("N", "V", basisSize, UtSW, basisSize, T_temp, ldT, hVals_a, NULL, hVals_b, NULL, basisSize, hVecs, ldhVecs, ctx));
CHKERR(Num_copy_matrix_Sprimme(&T[0], basisSize, basisSize, ldT, &T_temp[0], basisSize, ctx)); /* Ensure T is not overwritten */
CHKERR(Num_ggev_Sprimme("N", "V", basisSize, UtSW, basisSize, T_temp, basisSize, hVals_a, NULL, hVals_b, NULL, basisSize, hVecs, ldhVecs, ctx));

} /* End eigenproblem (with or without stabilization) */

Expand Down

0 comments on commit 2d2e1c1

Please sign in to comment.