Actual source code: lobpcg.c
slepc-3.17.0 2022-03-31
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "lobpcg"
13: Method: Locally Optimal Block Preconditioned Conjugate Gradient
15: Algorithm:
17: LOBPCG with soft and hard locking. Follows the implementation
18: in BLOPEX [2].
20: References:
22: [1] A. V. Knyazev, "Toward the optimal preconditioned eigensolver:
23: locally optimal block preconditioned conjugate gradient method",
24: SIAM J. Sci. Comput. 23(2):517-541, 2001.
26: [2] A. V. Knyazev et al., "Block Locally Optimal Preconditioned
27: Eigenvalue Xolvers (BLOPEX) in Hypre and PETSc", SIAM J. Sci.
28: Comput. 29(5):2224-2239, 2007.
29: */
31: #include <slepc/private/epsimpl.h>
33: typedef struct {
34: PetscInt bs; /* block size */
35: PetscBool lock; /* soft locking active/inactive */
36: PetscReal restart; /* restart parameter */
37: PetscInt guard; /* number of guard vectors */
38: } EPS_LOBPCG;
40: PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
41: {
42: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
43: PetscInt k;
45: k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
46: if (*ncv!=PETSC_DEFAULT) { /* ncv set */
48: } else *ncv = k;
49: if (*mpd==PETSC_DEFAULT) *mpd = 3*ctx->bs;
51: PetscFunctionReturn(0);
52: }
54: PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
55: {
56: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
58: EPSCheckHermitianDefinite(eps);
59: if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
61: EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd);
62: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
63: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
65: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
66: EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
68: if (!ctx->restart) ctx->restart = 0.9;
70: /* number of guard vectors */
71: if (ctx->bs==1) ctx->guard = 0;
72: else ctx->guard = PetscMin((PetscInt)((1.0-ctx->restart)*ctx->bs+0.45),ctx->bs-1);
74: EPSAllocateSolution(eps,0);
75: EPS_SetInnerProduct(eps);
76: DSSetType(eps->ds,DSGHEP);
77: DSAllocate(eps->ds,eps->mpd);
78: EPSSetWorkVecs(eps,1);
79: PetscFunctionReturn(0);
80: }
82: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
83: {
84: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
85: PetscInt i,j,k,ld,nv,ini,nmat,nc,nconv,locked,its,prev=0;
86: PetscReal norm;
87: PetscScalar *eigr,dot;
88: PetscBool breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
89: Mat A,B,M,V=NULL,W=NULL;
90: Vec v,z,w=eps->work[0];
91: BV X,Y=NULL,Z,R,P,AX,BX;
92: SlepcSC sc;
94: DSGetLeadingDimension(eps->ds,&ld);
95: STGetNumMatrices(eps->st,&nmat);
96: STGetMatrix(eps->st,0,&A);
97: if (nmat>1) STGetMatrix(eps->st,1,&B);
98: else B = NULL;
100: if (eps->which==EPS_LARGEST_REAL) { /* flip spectrum */
101: flip = PETSC_TRUE;
102: DSGetSlepcSC(eps->ds,&sc);
103: sc->comparison = SlepcCompareSmallestReal;
104: }
106: /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
107: PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL);
109: /* 1. Allocate memory */
110: PetscCalloc1(3*ctx->bs,&eigr);
111: BVDuplicateResize(eps->V,3*ctx->bs,&Z);
112: BVDuplicateResize(eps->V,ctx->bs,&X);
113: BVDuplicateResize(eps->V,ctx->bs,&R);
114: BVDuplicateResize(eps->V,ctx->bs,&P);
115: BVDuplicateResize(eps->V,ctx->bs,&AX);
116: if (B) BVDuplicateResize(eps->V,ctx->bs,&BX);
117: nc = eps->nds;
118: if (nc>0 || eps->nev>ctx->bs-ctx->guard) BVDuplicateResize(eps->V,nc+eps->nev,&Y);
119: if (nc>0) {
120: for (j=0;j<nc;j++) {
121: BVGetColumn(eps->V,-nc+j,&v);
122: BVInsertVec(Y,j,v);
123: BVRestoreColumn(eps->V,-nc+j,&v);
124: }
125: BVSetActiveColumns(Y,0,nc);
126: }
128: /* 2. Apply the constraints to the initial vectors */
129: /* 3. B-orthogonalize initial vectors */
130: for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
131: BVSetRandomColumn(eps->V,k);
132: BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
133: }
134: nv = ctx->bs;
135: BVSetActiveColumns(eps->V,0,nv);
136: BVSetActiveColumns(Z,0,nv);
137: BVCopy(eps->V,Z);
138: BVCopy(Z,X);
140: /* 4. Compute initial Ritz vectors */
141: BVMatMult(X,A,AX);
142: DSSetDimensions(eps->ds,nv,0,0);
143: DSGetMat(eps->ds,DS_MAT_A,&M);
144: BVMatProject(AX,NULL,X,M);
145: if (flip) MatScale(M,-1.0);
146: DSRestoreMat(eps->ds,DS_MAT_A,&M);
147: DSSetIdentity(eps->ds,DS_MAT_B);
148: DSSetState(eps->ds,DS_STATE_RAW);
149: DSSolve(eps->ds,eigr,NULL);
150: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
151: DSSynchronize(eps->ds,eigr,NULL);
152: for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
153: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
154: DSGetMat(eps->ds,DS_MAT_X,&M);
155: BVMultInPlace(X,M,0,nv);
156: BVMultInPlace(AX,M,0,nv);
157: MatDestroy(&M);
159: /* 5. Initialize range of active iterates */
160: locked = 0; /* hard-locked vectors, the leading locked columns of V are eigenvectors */
161: nconv = 0; /* number of converged eigenvalues in the current block */
162: its = 0; /* iterations for the current block */
164: /* 6. Main loop */
165: while (eps->reason == EPS_CONVERGED_ITERATING) {
167: if (ctx->lock) {
168: BVSetActiveColumns(R,nconv,ctx->bs);
169: BVSetActiveColumns(AX,nconv,ctx->bs);
170: if (B) BVSetActiveColumns(BX,nconv,ctx->bs);
171: }
173: /* 7. Compute residuals */
174: ini = (ctx->lock)? nconv: 0;
175: BVCopy(AX,R);
176: if (B) BVMatMult(X,B,BX);
177: for (j=ini;j<ctx->bs;j++) {
178: BVGetColumn(R,j,&v);
179: BVGetColumn(B?BX:X,j,&z);
180: VecAXPY(v,-eps->eigr[locked+j],z);
181: BVRestoreColumn(R,j,&v);
182: BVRestoreColumn(B?BX:X,j,&z);
183: }
185: /* 8. Compute residual norms and update index set of active iterates */
186: k = ini;
187: countc = PETSC_TRUE;
188: for (j=ini;j<ctx->bs;j++) {
189: i = locked+j;
190: BVGetColumn(R,j,&v);
191: VecNorm(v,NORM_2,&norm);
192: BVRestoreColumn(R,j,&v);
193: (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx);
194: if (countc) {
195: if (eps->errest[i] < eps->tol) k++;
196: else countc = PETSC_FALSE;
197: }
198: if (!countc && !eps->trackall) break;
199: }
200: nconv = k;
201: eps->nconv = locked + nconv;
202: if (its) EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs);
203: (*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
204: if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
205: BVSetActiveColumns(eps->V,locked,eps->nconv);
206: BVSetActiveColumns(X,0,nconv);
207: BVCopy(X,eps->V);
208: }
209: if (eps->reason != EPS_CONVERGED_ITERATING) {
210: break;
211: } else if (nconv >= ctx->bs-ctx->guard) {
212: eps->its += its-1;
213: its = 0;
214: } else its++;
216: if (nconv >= ctx->bs-ctx->guard) { /* force hard locking of vectors and compute new R */
218: /* extend constraints */
219: BVSetActiveColumns(Y,nc+locked,nc+locked+nconv);
220: BVCopy(X,Y);
221: BVSetActiveColumns(Y,0,nc+locked+nconv);
223: /* shift work BV's */
224: for (j=nconv;j<ctx->bs;j++) {
225: BVCopyColumn(X,j,j-nconv);
226: BVCopyColumn(R,j,j-nconv);
227: BVCopyColumn(P,j,j-nconv);
228: BVCopyColumn(AX,j,j-nconv);
229: if (B) BVCopyColumn(BX,j,j-nconv);
230: }
232: /* set new initial vectors */
233: BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv);
234: BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs);
235: BVCopy(eps->V,X);
236: for (j=ctx->bs-nconv;j<ctx->bs;j++) {
237: BVGetColumn(X,j,&v);
238: BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown);
239: if (norm>0.0 && !breakdown) VecScale(v,1.0/norm);
240: else {
241: PetscInfo(eps,"Orthogonalization of initial vector failed\n");
242: eps->reason = EPS_DIVERGED_BREAKDOWN;
243: goto diverged;
244: }
245: BVRestoreColumn(X,j,&v);
246: }
247: locked += nconv;
248: nconv = 0;
249: BVSetActiveColumns(X,nconv,ctx->bs);
251: /* B-orthogonalize initial vectors */
252: BVOrthogonalize(X,NULL);
253: BVSetActiveColumns(Z,nconv,ctx->bs);
254: BVSetActiveColumns(AX,nconv,ctx->bs);
255: BVCopy(X,Z);
257: /* compute initial Ritz vectors */
258: nv = ctx->bs;
259: BVMatMult(X,A,AX);
260: DSSetDimensions(eps->ds,nv,0,0);
261: DSGetMat(eps->ds,DS_MAT_A,&M);
262: BVMatProject(AX,NULL,X,M);
263: if (flip) MatScale(M,-1.0);
264: DSRestoreMat(eps->ds,DS_MAT_A,&M);
265: DSSetIdentity(eps->ds,DS_MAT_B);
266: DSSetState(eps->ds,DS_STATE_RAW);
267: DSSolve(eps->ds,eigr,NULL);
268: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
269: DSSynchronize(eps->ds,eigr,NULL);
270: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
271: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
272: DSGetMat(eps->ds,DS_MAT_X,&M);
273: BVMultInPlace(X,M,0,nv);
274: BVMultInPlace(AX,M,0,nv);
275: MatDestroy(&M);
277: continue; /* skip the rest of the iteration */
278: }
280: ini = (ctx->lock)? nconv: 0;
281: if (ctx->lock) {
282: BVSetActiveColumns(R,nconv,ctx->bs);
283: BVSetActiveColumns(P,nconv,ctx->bs);
284: BVSetActiveColumns(AX,nconv,ctx->bs);
285: if (B) BVSetActiveColumns(BX,nconv,ctx->bs);
286: }
288: /* 9. Apply preconditioner to the residuals */
289: BVGetMat(R,&V);
290: if (prev != ctx->bs-ini) {
291: prev = ctx->bs-ini;
292: MatDestroy(&W);
293: MatDuplicate(V,MAT_SHARE_NONZERO_PATTERN,&W);
294: }
295: STApplyMat(eps->st,V,W);
296: if (checkprecond) {
297: for (j=ini;j<ctx->bs;j++) {
298: MatDenseGetColumnVecRead(V,j-ini,&v);
299: MatDenseGetColumnVecRead(W,j-ini,&w);
300: VecDot(v,w,&dot);
301: MatDenseRestoreColumnVecRead(W,j-ini,&w);
302: MatDenseRestoreColumnVecRead(V,j-ini,&v);
303: if (PetscRealPart(dot)<0.0) {
304: PetscInfo(eps,"The preconditioner is not positive-definite\n");
305: eps->reason = EPS_DIVERGED_BREAKDOWN;
306: goto diverged;
307: }
308: }
309: }
310: if (nc+locked>0) {
311: for (j=ini;j<ctx->bs;j++) {
312: MatDenseGetColumnVecWrite(W,j-ini,&w);
313: BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
314: if (norm>0.0 && !breakdown) VecScale(w,1.0/norm);
315: MatDenseRestoreColumnVecWrite(W,j-ini,&w);
316: if (norm<=0.0 || breakdown) {
317: PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
318: eps->reason = EPS_DIVERGED_BREAKDOWN;
319: goto diverged;
320: }
321: }
322: }
323: MatCopy(W,V,SAME_NONZERO_PATTERN);
324: BVRestoreMat(R,&V);
326: /* 11. B-orthonormalize preconditioned residuals */
327: BVOrthogonalize(R,NULL);
329: /* 13-16. B-orthonormalize conjugate directions */
330: if (its>1) BVOrthogonalize(P,NULL);
332: /* 17-23. Compute symmetric Gram matrices */
333: BVSetActiveColumns(Z,0,ctx->bs);
334: BVSetActiveColumns(X,0,ctx->bs);
335: BVCopy(X,Z);
336: BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
337: BVCopy(R,Z);
338: if (its>1) {
339: BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
340: BVCopy(P,Z);
341: }
343: if (its>1) nv = 3*ctx->bs-2*ini;
344: else nv = 2*ctx->bs-ini;
346: BVSetActiveColumns(Z,0,nv);
347: DSSetDimensions(eps->ds,nv,0,0);
348: DSGetMat(eps->ds,DS_MAT_A,&M);
349: BVMatProject(Z,A,Z,M);
350: if (flip) MatScale(M,-1.0);
351: DSRestoreMat(eps->ds,DS_MAT_A,&M);
352: DSGetMat(eps->ds,DS_MAT_B,&M);
353: BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
354: DSRestoreMat(eps->ds,DS_MAT_B,&M);
356: /* 24. Solve the generalized eigenvalue problem */
357: DSSetState(eps->ds,DS_STATE_RAW);
358: DSSolve(eps->ds,eigr,NULL);
359: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
360: DSSynchronize(eps->ds,eigr,NULL);
361: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
362: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
364: /* 25-33. Compute Ritz vectors */
365: DSGetMat(eps->ds,DS_MAT_X,&M);
366: BVSetActiveColumns(Z,ctx->bs,nv);
367: if (ctx->lock) BVSetActiveColumns(P,0,ctx->bs);
368: BVMult(P,1.0,0.0,Z,M);
369: BVCopy(P,X);
370: if (ctx->lock) BVSetActiveColumns(P,nconv,ctx->bs);
371: BVSetActiveColumns(Z,0,ctx->bs);
372: BVMult(X,1.0,1.0,Z,M);
373: if (ctx->lock) BVSetActiveColumns(X,nconv,ctx->bs);
374: BVMatMult(X,A,AX);
375: MatDestroy(&M);
376: }
378: diverged:
379: eps->its += its;
381: if (flip) sc->comparison = SlepcCompareLargestReal;
382: PetscFree(eigr);
383: MatDestroy(&W);
384: if (V) BVRestoreMat(R,&V); /* only needed when goto diverged is reached */
385: BVDestroy(&Z);
386: BVDestroy(&X);
387: BVDestroy(&R);
388: BVDestroy(&P);
389: BVDestroy(&AX);
390: if (B) BVDestroy(&BX);
391: if (nc>0 || eps->nev>ctx->bs-ctx->guard) BVDestroy(&Y);
392: PetscFunctionReturn(0);
393: }
395: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
396: {
397: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
399: if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
401: if (ctx->bs != bs) {
402: ctx->bs = bs;
403: eps->state = EPS_STATE_INITIAL;
404: }
405: PetscFunctionReturn(0);
406: }
408: /*@
409: EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.
411: Logically Collective on eps
413: Input Parameters:
414: + eps - the eigenproblem solver context
415: - bs - the block size
417: Options Database Key:
418: . -eps_lobpcg_blocksize - Sets the block size
420: Level: advanced
422: .seealso: EPSLOBPCGGetBlockSize()
423: @*/
424: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
425: {
428: PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
429: PetscFunctionReturn(0);
430: }
432: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
433: {
434: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
436: *bs = ctx->bs;
437: PetscFunctionReturn(0);
438: }
440: /*@
441: EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.
443: Not Collective
445: Input Parameter:
446: . eps - the eigenproblem solver context
448: Output Parameter:
449: . bs - the block size
451: Level: advanced
453: .seealso: EPSLOBPCGSetBlockSize()
454: @*/
455: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
456: {
459: PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
460: PetscFunctionReturn(0);
461: }
463: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
464: {
465: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
467: if (restart==PETSC_DEFAULT) restart = 0.9;
469: if (restart != ctx->restart) {
470: ctx->restart = restart;
471: eps->state = EPS_STATE_INITIAL;
472: }
473: PetscFunctionReturn(0);
474: }
476: /*@
477: EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
479: Logically Collective on eps
481: Input Parameters:
482: + eps - the eigenproblem solver context
483: - restart - the percentage of the block of vectors to force a restart
485: Options Database Key:
486: . -eps_lobpcg_restart - Sets the restart parameter
488: Notes:
489: The meaning of this parameter is the proportion of vectors within the
490: current block iterate that must have converged in order to force a
491: restart with hard locking.
492: Allowed values are in the range [0.1,1.0]. The default is 0.9.
494: Level: advanced
496: .seealso: EPSLOBPCGGetRestart()
497: @*/
498: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
499: {
502: PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
503: PetscFunctionReturn(0);
504: }
506: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
507: {
508: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
510: *restart = ctx->restart;
511: PetscFunctionReturn(0);
512: }
514: /*@
515: EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.
517: Not Collective
519: Input Parameter:
520: . eps - the eigenproblem solver context
522: Output Parameter:
523: . restart - the restart parameter
525: Level: advanced
527: .seealso: EPSLOBPCGSetRestart()
528: @*/
529: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
530: {
533: PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
534: PetscFunctionReturn(0);
535: }
537: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
538: {
539: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
541: ctx->lock = lock;
542: PetscFunctionReturn(0);
543: }
545: /*@
546: EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
547: the LOBPCG method.
549: Logically Collective on eps
551: Input Parameters:
552: + eps - the eigenproblem solver context
553: - lock - true if the locking variant must be selected
555: Options Database Key:
556: . -eps_lobpcg_locking - Sets the locking flag
558: Notes:
559: This flag refers to soft locking (converged vectors within the current
560: block iterate), since hard locking is always used (when nev is larger
561: than the block size).
563: Level: advanced
565: .seealso: EPSLOBPCGGetLocking()
566: @*/
567: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
568: {
571: PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
572: PetscFunctionReturn(0);
573: }
575: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
576: {
577: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
579: *lock = ctx->lock;
580: PetscFunctionReturn(0);
581: }
583: /*@
584: EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.
586: Not Collective
588: Input Parameter:
589: . eps - the eigenproblem solver context
591: Output Parameter:
592: . lock - the locking flag
594: Level: advanced
596: .seealso: EPSLOBPCGSetLocking()
597: @*/
598: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
599: {
602: PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
603: PetscFunctionReturn(0);
604: }
606: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
607: {
608: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
609: PetscBool isascii;
611: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
612: if (isascii) {
613: PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs);
614: PetscViewerASCIIPrintf(viewer," restart parameter=%g (using %" PetscInt_FMT " guard vectors)\n",(double)ctx->restart,ctx->guard);
615: PetscViewerASCIIPrintf(viewer," soft locking %sactivated\n",ctx->lock?"":"de");
616: }
617: PetscFunctionReturn(0);
618: }
620: PetscErrorCode EPSSetFromOptions_LOBPCG(PetscOptionItems *PetscOptionsObject,EPS eps)
621: {
622: PetscBool lock,flg;
623: PetscInt bs;
624: PetscReal restart;
626: PetscOptionsHead(PetscOptionsObject,"EPS LOBPCG Options");
628: PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg);
629: if (flg) EPSLOBPCGSetBlockSize(eps,bs);
631: PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg);
632: if (flg) EPSLOBPCGSetRestart(eps,restart);
634: PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg);
635: if (flg) EPSLOBPCGSetLocking(eps,lock);
637: PetscOptionsTail();
638: PetscFunctionReturn(0);
639: }
641: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
642: {
643: PetscFree(eps->data);
644: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
645: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
646: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
647: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
648: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
649: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
650: PetscFunctionReturn(0);
651: }
653: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
654: {
655: EPS_LOBPCG *lobpcg;
657: PetscNewLog(eps,&lobpcg);
658: eps->data = (void*)lobpcg;
659: lobpcg->lock = PETSC_TRUE;
661: eps->useds = PETSC_TRUE;
662: eps->categ = EPS_CATEGORY_PRECOND;
664: eps->ops->solve = EPSSolve_LOBPCG;
665: eps->ops->setup = EPSSetUp_LOBPCG;
666: eps->ops->setupsort = EPSSetUpSort_Default;
667: eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
668: eps->ops->destroy = EPSDestroy_LOBPCG;
669: eps->ops->view = EPSView_LOBPCG;
670: eps->ops->backtransform = EPSBackTransform_Default;
671: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
673: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
674: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
675: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
676: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
677: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
678: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
679: PetscFunctionReturn(0);
680: }