Skip to content

Commit

Permalink
Missing a globalSum somewhere + need to change sketched restart in Da…
Browse files Browse the repository at this point in the history 
…vidson
  • Loading branch information
@Heatherms27
Heatherms27 committed Mar 25, 2024
1 parent 7c55148 commit 9d43bbc
Show file tree
Hide file tree
Showing 5 changed files with 46 additions and 43 deletions.
23 changes: 14 additions & 9 deletions src/eigs/lanczos.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,13 +155,16 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
/* Sketching Variables -----------------------------------------------------*/
SCALAR *V_temp; /* Copy of basis vectors for sketching */
SCALAR *SV; /* The sketched basis */
SCALAR *Q; /* The "Q" factor in the QR decomposition of SV */
SCALAR *T; /* The "R" factor in the QR decomposition of SV */
SCALAR *SW; /* The projected sketched basis */
SCALAR *S_vals; /* CSC Formatted Values */
PRIMME_INT *S_rows; /* CSC Formatted Rows */
PRIMME_INT nnzPerCol, ldSV, ldSW, ldT; /* Size and nnz of the sketching matrix */
PRIMME_INT nnzPerCol, ldSV, ldSW, ldQ, ldT; /* Size and nnz of the sketching matrix */
REAL *normalize_evecs;

nnzPerCol = ldSV = ldSW = ldQ = ldT = 0;

/* -------------------------------------------------------------- */
/* Allocate objects */
/* -------------------------------------------------------------- */
Expand Down Expand Up @@ -246,21 +249,23 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
if(primme->projectionParams.projection == primme_proj_sketched)
{
/* Default settings for sketch size and nnz per column. Based on Yuji and Tropp's manuscript */
ldSV = ldSW = primme->sketchingParams.sketchSize;
ldSV = ldSW = ldQ = primme->sketchingParams.sketchSize;
ldT = primme->maxBasisSize+maxBlockSize;
nnzPerCol = primme->sketchingParams.nnzPerCol;

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

CHKERR(Num_malloc_Sprimme(nnzPerCol*primme->nLocal, &S_vals, ctx));
CHKERR(Num_malloc_Sprimme(ldSV*(primme->maxBasisSize+maxBlockSize), &SV, ctx));
CHKERR(Num_malloc_Sprimme(ldQ*(primme->maxBasisSize+maxBlockSize), &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+maxBlockSize), &V_temp, ctx));

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

CHKERR(Num_zero_matrix_Sprimme(T, ldT, ldT, ldT, ctx));
CHKERR(Num_zero_matrix_Sprimme(SV, ldSV, primme->maxBasisSize+maxBlockSize, ldSV, ctx));

/* Build Sketch CSR Locally */
CHKERR(build_sketch_Sprimme(S_rows, S_vals, ctx));
Expand All @@ -277,7 +282,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

/* Update sketched basis if sketching is turned on */
if(primme->projectionParams.projection == primme_proj_sketched)
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, 0, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, 0, blockSize, S_rows, S_vals, ctx));

primme->stats.numOuterIterations++;

Expand All @@ -286,7 +291,6 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
* --------------------------------------------------------------- */
i = blockSize;
while(i < primme->maxBasisSize && (primme->maxOuterIterations == 0 || primme->stats.numOuterIterations < primme->maxOuterIterations) && primme->stats.numMatvecs < primme->maxMatvecs) {

blockSize = min(blockSize, primme->maxBasisSize - i); /* Adjust block size if needed */

CHKERR(ortho_Sprimme(&V[i*ldV], ldV, &H[(i-blockSize)*ldH + i], ldH, 0, blockSize-1, NULL, 0, 0, primme->nLocal, primme->iseed, ctx)); /* [V_i, b_i] = qr(V_i) */
Expand All @@ -309,7 +313,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

/* Update sketched basis if sketching is turned on */
if(primme->projectionParams.projection == primme_proj_sketched)
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, i, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, i, blockSize, S_rows, S_vals, ctx));

if(primme->printLevel >= 2 && (PRIMME_INT)(i+blockSize) % 100 == 0)
{
Expand All @@ -325,11 +329,11 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(ortho_Sprimme(&V_temp[ldV*(i+blockSize)], ldV, &H[i*ldH + (i+blockSize)], ldH, 0, blockSize-1, NULL, 0, 0, primme->nLocal, primme->iseed, ctx)); /* [V_i, b_i] = qr(V_i) */

/* SW = SV*H */
CHKERR(sketch_basis_Sprimme(V_temp, ldV, SV, ldSV, T, ldT, i+blockSize, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V_temp, ldV, SV, ldSV, Q, ldQ, T, ldT, i+blockSize, blockSize, S_rows, S_vals, ctx));
CHKERR(Num_gemm_Sprimme("N", "N", ldSV, i+blockSize, i+2*blockSize, 1.0, SV, ldSV, H, ldH, 0.0, SW, ldSW, ctx));

/* Getting our sketched basis and projected sketched basis */
CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, i+blockSize, hVals, i+blockSize, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, i+blockSize, hVals, i+blockSize, ctx));

} else { /* End sketching */
solve_timer = primme_wTimer();
Expand Down Expand Up @@ -431,11 +435,11 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(ortho_Sprimme(V, ldV, NULL, 0, primme->maxBasisSize, primme->maxBasisSize+blockSize-1, NULL, 0, 0, primme->nLocal, primme->iseed, ctx)); /* Orthogonalized the last block of V against the rest of the basis */

/* SW = SV*H */
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, primme->maxBasisSize, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, primme->maxBasisSize, blockSize, S_rows, S_vals, ctx));
CHKERR(Num_gemm_Sprimme("N", "N", ldSV, primme->maxBasisSize, primme->maxBasisSize+blockSize, 1.0, SV, ldSV, H, ldH, 0.0, SW, ldSW, ctx));

/* Performing sketched Rayleigh-Ritz */
CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, ldhVecs, hVals, primme->maxBasisSize, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, ldhVecs, hVals, primme->maxBasisSize, ctx));

} else { /* End sketched */
solve_timer = primme_wTimer();
Expand Down Expand Up @@ -495,6 +499,7 @@ int lanczos_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
{
CHKERR(Num_free_Sprimme(S_vals, 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_Rprimme(normalize_evecs, ctx));
Expand Down
34 changes: 19 additions & 15 deletions src/eigs/main_iter.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1593,14 +1593,15 @@ 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 *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, ldT, ldSW; /* Leading dimensions of SV, T, and SW */
PRIMME_INT ldSV, ldQ, ldT, ldSW; /* Leading dimensions of SV, T, and SW */


/* -------------------------------------------------------------- */
Expand Down Expand Up @@ -1653,19 +1654,21 @@ 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 = primme->sketchingParams.sketchSize;
ldSV = ldSW = ldQ = primme->sketchingParams.sketchSize;
ldT = primme->maxBasisSize;
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));


Expand Down Expand Up @@ -1756,10 +1759,10 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

targetShiftIndex = 0;

CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, NULL, 0, 0, basisSize, S_rows, S_vals, ctx));

CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));

numArbitraryVecs = 0;
maxRecentlyConverged = availableBlockSize = blockSize = 0;
Expand Down Expand Up @@ -1910,8 +1913,8 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
basisSize, blockSize, ctx));

/* Copy hVecs into prevhVecs */
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, basisSize, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, basisSize, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, basisSize, blockSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, NULL, 0, basisSize, blockSize, S_rows, S_vals, ctx));

CHKERR(Num_copy_matrix_SHprimme(hVecs, basisSize, basisSize,
basisSize, prevhVecs, primme->maxBasisSize, ctx));
Expand All @@ -1923,7 +1926,7 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
basisSize += blockSize;
blockSize = 0;

CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));

numArbitraryVecs = 0;
candidates_prepared = 0;
Expand Down Expand Up @@ -2123,7 +2126,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, 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, S_rows, S_vals, ctx));

restartsSinceReset++;

Expand Down Expand Up @@ -2163,11 +2166,11 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
CHKERR(matrixMatvec_Sprimme(V, primme->nLocal, ldV, W, ldW,
basisSize, numNew, ctx));

CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, basisSize, numNew, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, basisSize, numNew, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, basisSize, numNew, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, NULL, 0, basisSize, numNew, S_rows, S_vals, ctx));

basisSize += numNew;
CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));
}

primme->stats.numRestarts++;
Expand Down Expand Up @@ -2207,9 +2210,9 @@ int sketched_main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,

CHKERR(ortho_Sprimme(V, ldV, NULL, 0, 0, basisSize-1, NULL, 0, 0, primme->nLocal, primme->iseed, ctx));
CHKERR(matrixMatvec_Sprimme(V, primme->nLocal, ldV, W, ldW, 0, basisSize, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, NULL, 0, 0, basisSize, S_rows, S_vals, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, basisSize, hVals, basisSize, ctx));

CHKERR(verify_norms(V, ldV, W, ldW, NULL, 0, hVals,
restartLimitReached ? primme->numEvals : numConverged, resNorms,
Expand Down Expand Up @@ -2327,6 +2330,7 @@ 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));
Expand Down
8 changes: 4 additions & 4 deletions src/eigs/restart.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1773,7 +1773,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 *T, int ldT, SCALAR *SW, int ldSW, HSCALAR *hVecs,
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) {

Expand Down Expand Up @@ -1844,9 +1844,9 @@ STATIC int restart_sketched(SCALAR *V, int ldV, SCALAR *W, int ldW, SCALAR *SV,


// Update eiganpairs and residuals
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, T, ldT, 0, restartSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, 0, restartSize, S_rows, S_vals, ctx));
CHKERR(sketched_RR_Sprimme(SV, ldSV, T, ldT, SW, ldSW, hVecs, restartSize, hVals, restartSize, ctx));
CHKERR(sketch_basis_Sprimme(V, ldV, SV, ldSV, Q, ldQ, T, ldT, 0, restartSize, S_rows, S_vals, ctx));
CHKERR(sketch_basis_Sprimme(W, ldW, SW, ldSW, NULL, 0, NULL, 0, 0, restartSize, S_rows, S_vals, ctx));
CHKERR(sketched_RR_Sprimme(Q, ldQ, T, ldT, SW, ldSW, hVecs, restartSize, hVals, restartSize, ctx));

(*basisSize) = restartSize;

Expand Down
14 changes: 7 additions & 7 deletions src/eigs/sketch.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ int build_sketch_Sprimme(PRIMME_INT *S_rows, SCALAR *S_vals, primme_context ctx)
* - If "T" != NULL, then SV will be the "Q" factor
******************************************************************************/
TEMPLATE_PLEASE
int sketch_basis_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT, PRIMME_INT basisSize, PRIMME_INT blockSize, PRIMME_INT *S_rows, SCALAR *S_vals, primme_context ctx) {
int sketch_basis_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *SV, PRIMME_INT ldSV, SCALAR *Q, PRIMME_INT ldQ, SCALAR *T, PRIMME_INT ldT, PRIMME_INT basisSize, PRIMME_INT blockSize, PRIMME_INT *S_rows, SCALAR *S_vals, primme_context ctx) {

primme_params *primme = ctx.primme;

Expand All @@ -164,8 +164,9 @@ int sketch_basis_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *SV, PRIMME_INT ldSV,
CHKERR(globalSum_Sprimme(&SV[basisSize*ldSV], blockSize*ldSV, ctx));

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(SV, ldSV, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, ldSV, primme->iseed, ctx));
CHKERR(ortho_Sprimme(Q, ldQ, T, ldT, basisSize, basisSize+blockSize-1, NULL, 0, 0, ldQ, primme->iseed, ctx));
}

if (primme) primme->stats.timeSketchMatvec += primme_wTimer() - t0;
Expand Down Expand Up @@ -266,7 +267,7 @@ int sketched_RR_Sprimme(SCALAR *SV, PRIMME_INT ldSV, SCALAR *T, PRIMME_INT ldT,
REAL *sing_vals; /* Returned singular values from the SVD */
int *eval_perm; /* To sort the eigenpairs */

REAL cond_est; /* Estimation of the condition number of V - used to determine whether stabilization is needed */
REAL cond_est = 0; /* Estimation of the condition number of V - used to determine whether stabilization is needed */
PRIMME_INT tbasisSize = basisSize; /* Truncated basis size (for Stabilization) */

CHKERR(Num_malloc_Sprimme(basisSize*basisSize, &UtSW, ctx));
Expand All @@ -278,14 +279,13 @@ 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 */
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));

printf("Estimated Cond Num: %E\n", cond_est);
//printf("Estimated Cond Num: %E\n", 1/cond_est);

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

SCALAR *UtSWV; /* Left hand side of the generalized eigenvalue problem */
Expand Down Expand Up @@ -345,7 +345,7 @@ 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 */
assert(ldSV == ldSW);
CHKERR(Num_gemm_Sprimme("C", "N", basisSize, basisSize, ldSV, 1.0, SV, ldSV, SW, ldSW, 0.0, UtSW, basisSize, ctx));
Expand Down
10 changes: 2 additions & 8 deletions src/linalg/blaslapack.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1898,13 +1898,10 @@ int Num_trcon_Sprimme(const char *norm, const char *uplo, const char *diag, PRI
SCALAR *work;
PRIMME_BLASINT lwork = 2*n;

REAL lrcond = 0.0;

CHKERR(Num_malloc_Sprimme(lwork, &work, ctx));
CHKERR(Num_malloc_Rprimme(ln, &rwork, ctx));

XTRCON(norm, uplo, diag, &ln, a, &llda, &lrcond, work, rwork, &linfo);
*rcond = (REAL)lrcond;
XTRCON(norm, uplo, diag, &ln, a, &llda, rcond, work, rwork, &linfo);

CHKERR(Num_free_Sprimme(work, ctx));

Expand All @@ -1913,13 +1910,10 @@ int Num_trcon_Sprimme(const char *norm, const char *uplo, const char *diag, PRI
PRIMME_BLASINT *work;
PRIMME_BLASINT lwork = 3*n;

PRIMME_BLASINT lrcond = 0;

CHKERR(Num_malloc_Rprimme(lwork, &rwork, ctx));
CHKERR(Num_malloc_iblasprimme(ln, &work, ctx));

XTRCON(norm, uplo, diag, &ln, a, &llda, &lrcond, rwork, work, &linfo);
*rcond = (REAL)lrcond;
XTRCON(norm, uplo, diag, &ln, a, &llda, rcond, rwork, work, &linfo);

CHKERR(Num_free_iblasprimme(work, ctx));

Expand Down

0 comments on commit 9d43bbc

Please sign in to comment.