From eab3ff37e128a56b27fee533351c56b40e3fa0f9 Mon Sep 17 00:00:00 2001 From: Heather M Switzer Date: Fri, 22 Mar 2024 16:26:35 -0400 Subject: [PATCH] Debugging why the sketched aigensolver with no stabilization is producing inaccurate eigenvalues --- src/eigs/lanczos.c | 2 ++ src/eigs/sketch.c | 16 ++++++++++------ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/eigs/lanczos.c b/src/eigs/lanczos.c index ac1be92a..7a20932f 100644 --- a/src/eigs/lanczos.c +++ b/src/eigs/lanczos.c @@ -260,6 +260,8 @@ 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)); + /* Build Sketch CSR Locally */ CHKERR(build_sketch_Sprimme(S_rows, S_vals, ctx)); diff --git a/src/eigs/sketch.c b/src/eigs/sketch.c index 09a6246c..adaa5967 100644 --- a/src/eigs/sketch.c +++ b/src/eigs/sketch.c @@ -162,7 +162,6 @@ 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)); if (primme) primme->stats.timeSketchMatvec += primme_wTimer() - t0; @@ -258,6 +257,7 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT, SCALAR *hVals_b; /* The "beta" eigenvalues returned from LaPack's GGEV */ SCALAR normalize_hvecs; /* For normalizing the returned eigenvectors from LaPack's GGEV*/ SCALAR *UtSW; /* For solving the generalized eigenvalue problem */ + SCALAR *T_temp; /* For taking the SVD of T */ REAL *sing_vals; /* Returned singular values from the SVD */ int *eval_perm; /* To sort the eigenpairs */ @@ -268,20 +268,21 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT, CHKERR(Num_malloc_Sprimme(basisSize*basisSize, &UtSW, ctx)); CHKERR(Num_malloc_Sprimme(basisSize, &hVals_a, ctx)); CHKERR(Num_malloc_Sprimme(basisSize, &hVals_b, ctx)); + CHKERR(Num_malloc_Sprimme(basisSize*basisSize, &T_temp, ctx)); CHKERR(Num_malloc_Rprimme(basisSize, &sing_vals, ctx)); CHKERR(Num_malloc_iprimme(basisSize, &eval_perm, ctx)); /* Check condition number to see if stabilization is needed */ - CHKERR(Num_gesvd_Sprimme("N", "N", basisSize, basisSize, T, ldT, sing_vals, NULL, basisSize, NULL, basisSize, ctx)); + 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)); + 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) { - if(1) { - + if(cond_est > 1/MACHINE_EPSILON) { double stab_time = primme_wTimer(); // For monitoring /* Truncate the basis for stabilization */ @@ -341,10 +342,12 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT, /* XXX: Stabilization NOT needed */ } else { - /* Compute left-hand-side of the eigenproblem */ + /* Compute left-hand-side of the eigenproblem */ + assert(ldSV == ldSW); 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, ldT, hVals_a, NULL, hVals_b, NULL, basisSize, hVecs, ldhVecs, ctx)); } /* End eigenproblem (with or without stabilization) */ @@ -371,6 +374,7 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT, CHKERR(Num_free_Sprimme(hVals_a, ctx)); CHKERR(Num_free_Sprimme(hVals_b, ctx)); CHKERR(Num_free_iprimme(eval_perm, ctx)); + CHKERR(Num_free_Sprimme(T_temp, ctx)); } /* End if procID == 0 */