diff --git a/src/eigs/lanczos.c b/src/eigs/lanczos.c index 7a20932f..cbdf820c 100644 --- a/src/eigs/lanczos.c +++ b/src/eigs/lanczos.c @@ -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)); diff --git a/src/eigs/main_iter.c b/src/eigs/main_iter.c index 101949ea..18f251b2 100644 --- a/src/eigs/main_iter.c +++ b/src/eigs/main_iter.c @@ -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 */ @@ -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)); @@ -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)); + /* -------------------------------------------------------------- */ diff --git a/src/eigs/sketch.c b/src/eigs/sketch.c index c125dbb7..f0c1a6d3 100644 --- a/src/eigs/sketch.c +++ b/src/eigs/sketch.c @@ -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; @@ -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) */ @@ -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]; @@ -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) */