Actual source code: lobpcg.c
slepc-3.13.4 2020-09-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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;
46: k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
47: if (*ncv!=PETSC_DEFAULT) { /* ncv set */
48: if (*ncv<k) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv is not sufficiently large");
49: } else *ncv = k;
50: if (*mpd==PETSC_DEFAULT) *mpd = 3*ctx->bs;
51: else if (*mpd!=3*ctx->bs) SETERRQ(PetscObjectComm((PetscObject)eps),1,"This solver does not allow a value of mpd different from 3*blocksize");
52: return(0);
53: }
55: PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
56: {
58: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
59: PetscBool istrivial;
62: if (!eps->ishermitian || (eps->isgeneralized && !eps->ispositive)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"LOBPCG only works for Hermitian problems");
63: if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
64: if (eps->n-eps->nds<5*ctx->bs) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The problem size is too small relative to the block size");
65: EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd);
66: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
67: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
68: if (eps->which!=EPS_SMALLEST_REAL && eps->which!=EPS_LARGEST_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
69: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
70: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
71: RGIsTrivial(eps->rg,&istrivial);
72: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");
74: if (!ctx->restart) ctx->restart = 0.9;
76: /* number of guard vectors */
77: if (ctx->bs==1) ctx->guard = 0;
78: else ctx->guard = PetscMin((PetscInt)((1.0-ctx->restart)*ctx->bs+0.45),ctx->bs-1);
80: EPSAllocateSolution(eps,0);
81: EPS_SetInnerProduct(eps);
82: DSSetType(eps->ds,DSGHEP);
83: DSAllocate(eps->ds,eps->mpd);
84: EPSSetWorkVecs(eps,1);
85: return(0);
86: }
88: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
89: {
91: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
92: PetscInt i,j,k,ld,nv,ini,nmat,nc,nconv,locked,its;
93: PetscReal norm;
94: PetscScalar *eigr,dot;
95: PetscBool breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
96: Mat A,B,M;
97: Vec v,z,w=eps->work[0];
98: BV X,Y=NULL,Z,R,P,AX,BX;
99: SlepcSC sc;
102: DSGetLeadingDimension(eps->ds,&ld);
103: STGetNumMatrices(eps->st,&nmat);
104: STGetMatrix(eps->st,0,&A);
105: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
106: else B = NULL;
108: if (eps->which==EPS_LARGEST_REAL) { /* flip spectrum */
109: flip = PETSC_TRUE;
110: DSGetSlepcSC(eps->ds,&sc);
111: sc->comparison = SlepcCompareSmallestReal;
112: }
114: /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
115: PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL);
117: /* 1. Allocate memory */
118: PetscCalloc1(3*ctx->bs,&eigr);
119: BVDuplicateResize(eps->V,3*ctx->bs,&Z);
120: BVDuplicateResize(eps->V,ctx->bs,&X);
121: BVDuplicateResize(eps->V,ctx->bs,&R);
122: BVDuplicateResize(eps->V,ctx->bs,&P);
123: BVDuplicateResize(eps->V,ctx->bs,&AX);
124: if (B) {
125: BVDuplicateResize(eps->V,ctx->bs,&BX);
126: }
127: nc = eps->nds;
128: if (nc>0 || eps->nev>ctx->bs-ctx->guard) {
129: BVDuplicateResize(eps->V,nc+eps->nev,&Y);
130: }
131: if (nc>0) {
132: for (j=0;j<nc;j++) {
133: BVGetColumn(eps->V,-nc+j,&v);
134: BVInsertVec(Y,j,v);
135: BVRestoreColumn(eps->V,-nc+j,&v);
136: }
137: BVSetActiveColumns(Y,0,nc);
138: }
140: /* 2. Apply the constraints to the initial vectors */
141: /* 3. B-orthogonalize initial vectors */
142: for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
143: BVSetRandomColumn(eps->V,k);
144: BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
145: }
146: nv = ctx->bs;
147: BVSetActiveColumns(eps->V,0,nv);
148: BVSetActiveColumns(Z,0,nv);
149: BVCopy(eps->V,Z);
150: BVCopy(Z,X);
152: /* 4. Compute initial Ritz vectors */
153: BVMatMult(X,A,AX);
154: DSSetDimensions(eps->ds,nv,0,0,0);
155: DSGetMat(eps->ds,DS_MAT_A,&M);
156: BVMatProject(AX,NULL,X,M);
157: if (flip) { MatScale(M,-1.0); }
158: DSRestoreMat(eps->ds,DS_MAT_A,&M);
159: DSSetIdentity(eps->ds,DS_MAT_B);
160: DSSetState(eps->ds,DS_STATE_RAW);
161: DSSolve(eps->ds,eigr,NULL);
162: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
163: DSSynchronize(eps->ds,eigr,NULL);
164: for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
165: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
166: DSGetMat(eps->ds,DS_MAT_X,&M);
167: BVMultInPlace(X,M,0,nv);
168: BVMultInPlace(AX,M,0,nv);
169: MatDestroy(&M);
171: /* 5. Initialize range of active iterates */
172: locked = 0; /* hard-locked vectors, the leading locked columns of V are eigenvectors */
173: nconv = 0; /* number of converged eigenvalues in the current block */
174: its = 0; /* iterations for the current block */
176: /* 6. Main loop */
177: while (eps->reason == EPS_CONVERGED_ITERATING) {
179: if (ctx->lock) {
180: BVSetActiveColumns(R,nconv,ctx->bs);
181: BVSetActiveColumns(AX,nconv,ctx->bs);
182: if (B) {
183: BVSetActiveColumns(BX,nconv,ctx->bs);
184: }
185: }
187: /* 7. Compute residuals */
188: ini = (ctx->lock)? nconv: 0;
189: BVCopy(AX,R);
190: if (B) { BVMatMult(X,B,BX); }
191: for (j=ini;j<ctx->bs;j++) {
192: BVGetColumn(R,j,&v);
193: BVGetColumn(B?BX:X,j,&z);
194: VecAXPY(v,-eps->eigr[locked+j],z);
195: BVRestoreColumn(R,j,&v);
196: BVRestoreColumn(B?BX:X,j,&z);
197: }
199: /* 8. Compute residual norms and update index set of active iterates */
200: k = ini;
201: countc = PETSC_TRUE;
202: for (j=ini;j<ctx->bs;j++) {
203: i = locked+j;
204: BVGetColumn(R,j,&v);
205: VecNorm(v,NORM_2,&norm);
206: BVRestoreColumn(R,j,&v);
207: (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx);
208: if (countc) {
209: if (eps->errest[i] < eps->tol) k++;
210: else countc = PETSC_FALSE;
211: }
212: if (!countc && !eps->trackall) break;
213: }
214: nconv = k;
215: eps->nconv = locked + nconv;
216: if (its) {
217: EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs);
218: }
219: (*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
220: if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
221: BVSetActiveColumns(eps->V,locked,eps->nconv);
222: BVSetActiveColumns(X,0,nconv);
223: BVCopy(X,eps->V);
224: }
225: if (eps->reason != EPS_CONVERGED_ITERATING) {
226: break;
227: } else if (nconv >= ctx->bs-ctx->guard) {
228: eps->its += its-1;
229: its = 0;
230: } else its++;
232: if (nconv >= ctx->bs-ctx->guard) { /* force hard locking of vectors and compute new R */
234: /* extend constraints */
235: BVSetActiveColumns(Y,nc+locked,nc+locked+nconv);
236: BVCopy(X,Y);
237: BVSetActiveColumns(Y,0,nc+locked+nconv);
239: /* shift work BV's */
240: for (j=nconv;j<ctx->bs;j++) {
241: BVCopyColumn(X,j,j-nconv);
242: BVCopyColumn(R,j,j-nconv);
243: BVCopyColumn(P,j,j-nconv);
244: BVCopyColumn(AX,j,j-nconv);
245: if (B) {
246: BVCopyColumn(BX,j,j-nconv);
247: }
248: }
250: /* set new initial vectors */
251: BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv);
252: BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs);
253: BVCopy(eps->V,X);
254: for (j=ctx->bs-nconv;j<ctx->bs;j++) {
255: BVGetColumn(X,j,&v);
256: BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown);
257: if (norm>0.0 && !breakdown) {
258: VecScale(v,1.0/norm);
259: } else {
260: PetscInfo(eps,"Orthogonalization of initial vector failed\n");
261: eps->reason = EPS_DIVERGED_BREAKDOWN;
262: goto diverged;
263: }
264: BVRestoreColumn(X,j,&v);
265: }
266: locked += nconv;
267: nconv = 0;
268: BVSetActiveColumns(X,nconv,ctx->bs);
270: /* B-orthogonalize initial vectors */
271: BVOrthogonalize(X,NULL);
272: BVSetActiveColumns(Z,nconv,ctx->bs);
273: BVSetActiveColumns(AX,nconv,ctx->bs);
274: BVCopy(X,Z);
276: /* compute initial Ritz vectors */
277: nv = ctx->bs;
278: BVMatMult(X,A,AX);
279: DSSetDimensions(eps->ds,nv,0,0,0);
280: DSGetMat(eps->ds,DS_MAT_A,&M);
281: BVMatProject(AX,NULL,X,M);
282: if (flip) { MatScale(M,-1.0); }
283: DSRestoreMat(eps->ds,DS_MAT_A,&M);
284: DSSetIdentity(eps->ds,DS_MAT_B);
285: DSSetState(eps->ds,DS_STATE_RAW);
286: DSSolve(eps->ds,eigr,NULL);
287: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
288: DSSynchronize(eps->ds,eigr,NULL);
289: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
290: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
291: DSGetMat(eps->ds,DS_MAT_X,&M);
292: BVMultInPlace(X,M,0,nv);
293: BVMultInPlace(AX,M,0,nv);
294: MatDestroy(&M);
296: continue; /* skip the rest of the iteration */
297: }
299: ini = (ctx->lock)? nconv: 0;
300: if (ctx->lock) {
301: BVSetActiveColumns(R,nconv,ctx->bs);
302: BVSetActiveColumns(P,nconv,ctx->bs);
303: BVSetActiveColumns(AX,nconv,ctx->bs);
304: if (B) {
305: BVSetActiveColumns(BX,nconv,ctx->bs);
306: }
307: }
309: /* 9. Apply preconditioner to the residuals */
310: for (j=ini;j<ctx->bs;j++) {
311: BVGetColumn(R,j,&v);
312: STApply(eps->st,v,w);
313: if (checkprecond) {
314: VecDot(v,w,&dot);
315: if (PetscRealPart(dot)<0.0) {
316: PetscInfo(eps,"The preconditioner is not positive-definite\n");
317: eps->reason = EPS_DIVERGED_BREAKDOWN;
318: goto diverged;
319: }
320: }
321: if (nc+locked>0) {
322: BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
323: if (norm>0.0 && !breakdown) {
324: VecScale(w,1.0/norm);
325: } else {
326: PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
327: eps->reason = EPS_DIVERGED_BREAKDOWN;
328: goto diverged;
329: }
330: }
331: VecCopy(w,v);
332: BVRestoreColumn(R,j,&v);
333: }
335: /* 11. B-orthonormalize preconditioned residuals */
336: BVOrthogonalize(R,NULL);
338: /* 13-16. B-orthonormalize conjugate directions */
339: if (its>1) {
340: BVOrthogonalize(P,NULL);
341: }
343: /* 17-23. Compute symmetric Gram matrices */
344: BVSetActiveColumns(Z,0,ctx->bs);
345: BVSetActiveColumns(X,0,ctx->bs);
346: BVCopy(X,Z);
347: BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
348: BVCopy(R,Z);
349: if (its>1) {
350: BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
351: BVCopy(P,Z);
352: }
354: if (its>1) nv = 3*ctx->bs-2*ini;
355: else nv = 2*ctx->bs-ini;
357: BVSetActiveColumns(Z,0,nv);
358: DSSetDimensions(eps->ds,nv,0,0,0);
359: DSGetMat(eps->ds,DS_MAT_A,&M);
360: BVMatProject(Z,A,Z,M);
361: if (flip) { MatScale(M,-1.0); }
362: DSRestoreMat(eps->ds,DS_MAT_A,&M);
363: DSGetMat(eps->ds,DS_MAT_B,&M);
364: BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
365: DSRestoreMat(eps->ds,DS_MAT_B,&M);
367: /* 24. Solve the generalized eigenvalue problem */
368: DSSetState(eps->ds,DS_STATE_RAW);
369: DSSolve(eps->ds,eigr,NULL);
370: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
371: DSSynchronize(eps->ds,eigr,NULL);
372: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
373: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
375: /* 25-33. Compute Ritz vectors */
376: DSGetMat(eps->ds,DS_MAT_X,&M);
377: BVSetActiveColumns(Z,ctx->bs,nv);
378: if (ctx->lock) {
379: BVSetActiveColumns(P,0,ctx->bs);
380: }
381: BVMult(P,1.0,0.0,Z,M);
382: BVCopy(P,X);
383: if (ctx->lock) {
384: BVSetActiveColumns(P,nconv,ctx->bs);
385: }
386: BVSetActiveColumns(Z,0,ctx->bs);
387: BVMult(X,1.0,1.0,Z,M);
388: if (ctx->lock) {
389: BVSetActiveColumns(X,nconv,ctx->bs);
390: }
391: BVMatMult(X,A,AX);
392: MatDestroy(&M);
393: }
395: diverged:
396: eps->its += its;
398: if (flip) sc->comparison = SlepcCompareLargestReal;
399: PetscFree(eigr);
400: BVDestroy(&Z);
401: BVDestroy(&X);
402: BVDestroy(&R);
403: BVDestroy(&P);
404: BVDestroy(&AX);
405: if (B) {
406: BVDestroy(&BX);
407: }
408: if (nc>0 || eps->nev>ctx->bs-ctx->guard) {
409: BVDestroy(&Y);
410: }
411: return(0);
412: }
414: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
415: {
416: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
419: if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
420: if (bs <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %D",bs);
421: if (ctx->bs != bs) {
422: ctx->bs = bs;
423: eps->state = EPS_STATE_INITIAL;
424: }
425: return(0);
426: }
428: /*@
429: EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.
431: Logically Collective on eps
433: Input Parameters:
434: + eps - the eigenproblem solver context
435: - bs - the block size
437: Options Database Key:
438: . -eps_lobpcg_blocksize - Sets the block size
440: Level: advanced
442: .seealso: EPSLOBPCGGetBlockSize()
443: @*/
444: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
445: {
451: PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
452: return(0);
453: }
455: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
456: {
457: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
460: *bs = ctx->bs;
461: return(0);
462: }
464: /*@
465: EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.
467: Not Collective
469: Input Parameter:
470: . eps - the eigenproblem solver context
472: Output Parameter:
473: . bs - the block size
475: Level: advanced
477: .seealso: EPSLOBPCGSetBlockSize()
478: @*/
479: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
480: {
486: PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
487: return(0);
488: }
490: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
491: {
492: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
495: if (restart==PETSC_DEFAULT) restart = 0.9;
496: if (restart<0.1 || restart>1.0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The restart argument %g must be in the range [0.1,1.0]",(double)restart);
497: if (restart != ctx->restart) {
498: ctx->restart = restart;
499: eps->state = EPS_STATE_INITIAL;
500: }
501: return(0);
502: }
504: /*@
505: EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
506: The meaning of this parameter is the proportion of vectors within the
507: current block iterate that must have converged in order to force a
508: restart with hard locking.
510: Logically Collective on eps
512: Input Parameters:
513: + eps - the eigenproblem solver context
514: - restart - the percentage of the block of vectors to force a restart
516: Options Database Key:
517: . -eps_lobpcg_restart - Sets the restart parameter
519: Notes:
520: Allowed values are in the range [0.1,1.0]. The default is 0.6.
522: Level: advanced
524: .seealso: EPSLOBPCGGetRestart()
525: @*/
526: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
527: {
533: PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
534: return(0);
535: }
537: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
538: {
539: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
542: *restart = ctx->restart;
543: return(0);
544: }
546: /*@
547: EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.
549: Not Collective
551: Input Parameter:
552: . eps - the eigenproblem solver context
554: Output Parameter:
555: . restart - the restart parameter
557: Level: advanced
559: .seealso: EPSLOBPCGSetRestart()
560: @*/
561: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
562: {
568: PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
569: return(0);
570: }
572: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
573: {
574: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
577: ctx->lock = lock;
578: return(0);
579: }
581: /*@
582: EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
583: the LOBPCG method.
585: Logically Collective on eps
587: Input Parameters:
588: + eps - the eigenproblem solver context
589: - lock - true if the locking variant must be selected
591: Options Database Key:
592: . -eps_lobpcg_locking - Sets the locking flag
594: Notes:
595: This flag refers to soft locking (converged vectors within the current
596: block iterate), since hard locking is always used (when nev is larger
597: than the block size).
599: Level: advanced
601: .seealso: EPSLOBPCGGetLocking()
602: @*/
603: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
604: {
610: PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
611: return(0);
612: }
614: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
615: {
616: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
619: *lock = ctx->lock;
620: return(0);
621: }
623: /*@
624: EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.
626: Not Collective
628: Input Parameter:
629: . eps - the eigenproblem solver context
631: Output Parameter:
632: . lock - the locking flag
634: Level: advanced
636: .seealso: EPSLOBPCGSetLocking()
637: @*/
638: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
639: {
645: PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
646: return(0);
647: }
649: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
650: {
652: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
653: PetscBool isascii;
656: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
657: if (isascii) {
658: PetscViewerASCIIPrintf(viewer," block size %D\n",ctx->bs);
659: PetscViewerASCIIPrintf(viewer," restart parameter=%g (using %D guard vectors)\n",(double)ctx->restart,ctx->guard);
660: PetscViewerASCIIPrintf(viewer," soft locking %sactivated\n",ctx->lock?"":"de");
661: }
662: return(0);
663: }
665: PetscErrorCode EPSSetFromOptions_LOBPCG(PetscOptionItems *PetscOptionsObject,EPS eps)
666: {
668: PetscBool lock,flg;
669: PetscInt bs;
670: PetscReal restart;
673: PetscOptionsHead(PetscOptionsObject,"EPS LOBPCG Options");
675: PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg);
676: if (flg) { EPSLOBPCGSetBlockSize(eps,bs); }
678: PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg);
679: if (flg) { EPSLOBPCGSetRestart(eps,restart); }
681: PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg);
682: if (flg) { EPSLOBPCGSetLocking(eps,lock); }
684: PetscOptionsTail();
685: return(0);
686: }
688: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
689: {
693: PetscFree(eps->data);
694: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
695: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
696: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
697: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
698: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
699: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
700: return(0);
701: }
703: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
704: {
705: EPS_LOBPCG *lobpcg;
709: PetscNewLog(eps,&lobpcg);
710: eps->data = (void*)lobpcg;
711: lobpcg->lock = PETSC_TRUE;
713: eps->useds = PETSC_TRUE;
714: eps->categ = EPS_CATEGORY_PRECOND;
716: eps->ops->solve = EPSSolve_LOBPCG;
717: eps->ops->setup = EPSSetUp_LOBPCG;
718: eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
719: eps->ops->destroy = EPSDestroy_LOBPCG;
720: eps->ops->view = EPSView_LOBPCG;
721: eps->ops->backtransform = EPSBackTransform_Default;
722: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
724: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
725: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
726: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
727: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
728: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
729: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
730: return(0);
731: }