Actual source code: lobpcg.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepceps.h" I*/
33: typedef struct {
34: PetscInt bs; /* block size */
35: PetscBool lock; /* soft locking active/inactive */
36: PetscReal restart; /* restart parameter */
37: } EPS_LOBPCG;
39: PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
40: {
41: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
42: PetscInt k;
45: k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
46: if (*ncv) { /* ncv set */
47: if (*ncv<k) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv is not sufficiently large");
48: } else *ncv = k;
49: if (!*mpd) *mpd = 3*ctx->bs;
50: else if (*mpd!=3*ctx->bs) SETERRQ(PetscObjectComm((PetscObject)eps),1,"This solver does not allow a value of mpd different from 3*blocksize");
51: return(0);
52: }
54: PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
55: {
57: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
58: PetscBool istrivial;
61: if (!eps->ishermitian || (eps->isgeneralized && !eps->ispositive)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"LOBPCG only works for Hermitian problems");
62: if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
63: 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");
64: EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd);
65: if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
66: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
67: if (eps->which!=EPS_SMALLEST_REAL && eps->which!=EPS_LARGEST_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
68: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
69: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
70: RGIsTrivial(eps->rg,&istrivial);
71: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");
73: if (!ctx->restart) ctx->restart = 0.9;
75: EPSAllocateSolution(eps,0);
76: EPS_SetInnerProduct(eps);
77: DSSetType(eps->ds,DSGHEP);
78: DSAllocate(eps->ds,eps->mpd);
79: EPSSetWorkVecs(eps,1);
80: return(0);
81: }
83: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
84: {
86: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
87: PetscInt i,j,k,ld,nv,ini,nmat,nc,nconv,locked,guard,its;
88: PetscReal norm;
89: PetscScalar *eigr,dot;
90: PetscBool breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
91: Mat A,B,M;
92: Vec v,z,w=eps->work[0];
93: BV X,Y,Z,R,P,AX,BX;
94: SlepcSC sc;
97: DSGetLeadingDimension(eps->ds,&ld);
98: STGetNumMatrices(eps->st,&nmat);
99: STGetMatrix(eps->st,0,&A);
100: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
101: else B = NULL;
103: guard = (PetscInt)PetscRoundReal((1.0-ctx->restart)*ctx->bs); /* number of guard vectors */
105: if (eps->which==EPS_LARGEST_REAL) { /* flip spectrum */
106: flip = PETSC_TRUE;
107: DSGetSlepcSC(eps->ds,&sc);
108: sc->comparison = SlepcCompareSmallestReal;
109: }
111: /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
112: PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL);
114: /* 1. Allocate memory */
115: PetscCalloc1(3*ctx->bs,&eigr);
116: BVDuplicateResize(eps->V,3*ctx->bs,&Z);
117: BVDuplicateResize(eps->V,ctx->bs,&X);
118: BVDuplicateResize(eps->V,ctx->bs,&R);
119: BVDuplicateResize(eps->V,ctx->bs,&P);
120: BVDuplicateResize(eps->V,ctx->bs,&AX);
121: if (B) {
122: BVDuplicateResize(eps->V,ctx->bs,&BX);
123: }
124: nc = eps->nds;
125: if (nc>0 || eps->nev>ctx->bs-guard) {
126: BVDuplicateResize(eps->V,nc+eps->nev,&Y);
127: }
128: if (nc>0) {
129: for (j=0;j<nc;j++) {
130: BVGetColumn(eps->V,-nc+j,&v);
131: BVInsertVec(Y,j,v);
132: BVRestoreColumn(eps->V,-nc+j,&v);
133: }
134: BVSetActiveColumns(Y,0,nc);
135: }
137: /* 2. Apply the constraints to the initial vectors */
138: /* 3. B-orthogonalize initial vectors */
139: for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
140: BVSetRandomColumn(eps->V,k);
141: BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
142: }
143: nv = ctx->bs;
144: BVSetActiveColumns(eps->V,0,nv);
145: BVSetActiveColumns(Z,0,nv);
146: BVCopy(eps->V,Z);
147: BVCopy(Z,X);
149: /* 4. Compute initial Ritz vectors */
150: BVMatMult(X,A,AX);
151: DSSetDimensions(eps->ds,nv,0,0,0);
152: DSGetMat(eps->ds,DS_MAT_A,&M);
153: BVMatProject(AX,NULL,X,M);
154: if (flip) { MatScale(M,-1.0); }
155: DSRestoreMat(eps->ds,DS_MAT_A,&M);
156: DSSetIdentity(eps->ds,DS_MAT_B);
157: DSSetState(eps->ds,DS_STATE_RAW);
158: DSSolve(eps->ds,eigr,NULL);
159: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
160: DSSynchronize(eps->ds,eigr,NULL);
161: for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
162: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
163: DSGetMat(eps->ds,DS_MAT_X,&M);
164: BVMultInPlace(X,M,0,nv);
165: BVMultInPlace(AX,M,0,nv);
166: MatDestroy(&M);
168: /* 5. Initialize range of active iterates */
169: locked = 0; /* hard-locked vectors, the leading locked columns of V are eigenvectors */
170: nconv = 0; /* number of converged eigenvalues in the current block */
171: its = 0; /* iterations for the current block */
173: /* 6. Main loop */
174: while (eps->reason == EPS_CONVERGED_ITERATING) {
176: if (ctx->lock) {
177: BVSetActiveColumns(R,nconv,ctx->bs);
178: BVSetActiveColumns(AX,nconv,ctx->bs);
179: if (B) {
180: BVSetActiveColumns(BX,nconv,ctx->bs);
181: }
182: }
184: /* 7. Compute residuals */
185: ini = (ctx->lock)? nconv: 0;
186: BVCopy(AX,R);
187: if (B) { BVMatMult(X,B,BX); }
188: for (j=ini;j<ctx->bs;j++) {
189: BVGetColumn(R,j,&v);
190: BVGetColumn(B?BX:X,j,&z);
191: VecAXPY(v,-eps->eigr[locked+j],z);
192: BVRestoreColumn(R,j,&v);
193: BVRestoreColumn(B?BX:X,j,&z);
194: }
196: /* 8. Compute residual norms and update index set of active iterates */
197: k = ini;
198: countc = PETSC_TRUE;
199: for (j=ini;j<ctx->bs;j++) {
200: i = locked+j;
201: BVGetColumn(R,j,&v);
202: VecNorm(v,NORM_2,&norm);
203: BVRestoreColumn(R,j,&v);
204: (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx);
205: if (countc) {
206: if (eps->errest[i] < eps->tol) k++;
207: else countc = PETSC_FALSE;
208: }
209: if (!countc && !eps->trackall) break;
210: }
211: nconv = k;
212: eps->nconv = locked + nconv;
213: if (its) {
214: EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs);
215: }
216: (*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
217: if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-guard) {
218: BVSetActiveColumns(eps->V,locked,eps->nconv);
219: BVSetActiveColumns(X,0,nconv);
220: BVCopy(X,eps->V);
221: }
222: if (eps->reason != EPS_CONVERGED_ITERATING) {
223: break;
224: } else if (nconv >= ctx->bs-guard) {
225: eps->its += its-1;
226: its = 0;
227: } else its++;
229: if (nconv >= ctx->bs-guard) { /* force hard locking of vectors and compute new R */
231: /* extend constraints */
232: BVSetActiveColumns(Y,nc+locked,nc+locked+nconv);
233: BVCopy(X,Y);
234: BVSetActiveColumns(Y,0,nc+locked+nconv);
236: /* shift work BV's */
237: for (j=nconv;j<ctx->bs;j++) {
238: BVCopyColumn(X,j,j-nconv);
239: BVCopyColumn(R,j,j-nconv);
240: BVCopyColumn(P,j,j-nconv);
241: BVCopyColumn(AX,j,j-nconv);
242: if (B) {
243: BVCopyColumn(BX,j,j-nconv);
244: }
245: }
247: /* set new initial vectors */
248: BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv);
249: BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs);
250: BVCopy(eps->V,X);
251: for (j=ctx->bs-nconv;j<ctx->bs;j++) {
252: BVGetColumn(X,j,&v);
253: BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown);
254: if (norm>0.0 && !breakdown) {
255: VecScale(v,1.0/norm);
256: } else {
257: PetscInfo(eps,"Orthogonalization of initial vector failed\n");
258: eps->reason = EPS_DIVERGED_BREAKDOWN;
259: goto diverged;
260: }
261: BVRestoreColumn(X,j,&v);
262: }
263: locked += nconv;
264: nconv = 0;
265: BVSetActiveColumns(X,nconv,ctx->bs);
267: /* B-orthogonalize initial vectors */
268: BVOrthogonalize(X,NULL);
269: BVSetActiveColumns(Z,nconv,ctx->bs);
270: BVSetActiveColumns(AX,nconv,ctx->bs);
271: BVCopy(X,Z);
273: /* compute initial Ritz vectors */
274: nv = ctx->bs;
275: BVMatMult(X,A,AX);
276: DSSetDimensions(eps->ds,nv,0,0,0);
277: DSGetMat(eps->ds,DS_MAT_A,&M);
278: BVMatProject(AX,NULL,X,M);
279: if (flip) { MatScale(M,-1.0); }
280: DSRestoreMat(eps->ds,DS_MAT_A,&M);
281: DSSetIdentity(eps->ds,DS_MAT_B);
282: DSSetState(eps->ds,DS_STATE_RAW);
283: DSSolve(eps->ds,eigr,NULL);
284: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
285: DSSynchronize(eps->ds,eigr,NULL);
286: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
287: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
288: DSGetMat(eps->ds,DS_MAT_X,&M);
289: BVMultInPlace(X,M,0,nv);
290: BVMultInPlace(AX,M,0,nv);
291: MatDestroy(&M);
293: continue; /* skip the rest of the iteration */
294: }
296: ini = (ctx->lock)? nconv: 0;
297: if (ctx->lock) {
298: BVSetActiveColumns(R,nconv,ctx->bs);
299: BVSetActiveColumns(P,nconv,ctx->bs);
300: BVSetActiveColumns(AX,nconv,ctx->bs);
301: if (B) {
302: BVSetActiveColumns(BX,nconv,ctx->bs);
303: }
304: }
306: /* 9. Apply preconditioner to the residuals */
307: for (j=ini;j<ctx->bs;j++) {
308: BVGetColumn(R,j,&v);
309: STMatSolve(eps->st,v,w);
310: if (checkprecond) {
311: VecDot(v,w,&dot);
312: if (PetscRealPart(dot)<0.0) {
313: PetscInfo(eps,"The preconditioner is not positive-definite\n");
314: eps->reason = EPS_DIVERGED_BREAKDOWN;
315: goto diverged;
316: }
317: }
318: if (nc+locked>0) {
319: BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
320: if (norm>0.0 && !breakdown) {
321: VecScale(w,1.0/norm);
322: } else {
323: PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
324: eps->reason = EPS_DIVERGED_BREAKDOWN;
325: goto diverged;
326: }
327: }
328: VecCopy(w,v);
329: BVRestoreColumn(R,j,&v);
330: }
332: /* 11. B-orthonormalize preconditioned residuals */
333: BVOrthogonalize(R,NULL);
335: /* 13-16. B-orthonormalize conjugate directions */
336: if (its>1) {
337: BVOrthogonalize(P,NULL);
338: }
340: /* 17-23. Compute symmetric Gram matrices */
341: BVSetActiveColumns(Z,0,ctx->bs);
342: BVSetActiveColumns(X,0,ctx->bs);
343: BVCopy(X,Z);
344: BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
345: BVCopy(R,Z);
346: if (its>1) {
347: BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
348: BVCopy(P,Z);
349: }
351: if (its>1) nv = 3*ctx->bs-2*ini;
352: else nv = 2*ctx->bs-ini;
354: BVSetActiveColumns(Z,0,nv);
355: DSSetDimensions(eps->ds,nv,0,0,0);
356: DSGetMat(eps->ds,DS_MAT_A,&M);
357: BVMatProject(Z,A,Z,M);
358: if (flip) { MatScale(M,-1.0); }
359: DSRestoreMat(eps->ds,DS_MAT_A,&M);
360: DSGetMat(eps->ds,DS_MAT_B,&M);
361: BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
362: DSRestoreMat(eps->ds,DS_MAT_B,&M);
364: /* 24. Solve the generalized eigenvalue problem */
365: DSSetState(eps->ds,DS_STATE_RAW);
366: DSSolve(eps->ds,eigr,NULL);
367: DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
368: DSSynchronize(eps->ds,eigr,NULL);
369: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
370: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
372: /* 25-33. Compute Ritz vectors */
373: DSGetMat(eps->ds,DS_MAT_X,&M);
374: BVSetActiveColumns(Z,ctx->bs,nv);
375: if (ctx->lock) {
376: BVSetActiveColumns(P,0,ctx->bs);
377: }
378: BVMult(P,1.0,0.0,Z,M);
379: BVCopy(P,X);
380: if (ctx->lock) {
381: BVSetActiveColumns(P,nconv,ctx->bs);
382: }
383: BVSetActiveColumns(Z,0,ctx->bs);
384: BVMult(X,1.0,1.0,Z,M);
385: if (ctx->lock) {
386: BVSetActiveColumns(X,nconv,ctx->bs);
387: }
388: BVMatMult(X,A,AX);
389: MatDestroy(&M);
390: }
392: diverged:
393: eps->its += its;
395: if (flip) sc->comparison = SlepcCompareLargestReal;
396: PetscFree(eigr);
397: BVDestroy(&Z);
398: BVDestroy(&X);
399: BVDestroy(&R);
400: BVDestroy(&P);
401: BVDestroy(&AX);
402: if (B) {
403: BVDestroy(&BX);
404: }
405: if (nc>0 || eps->nev>ctx->bs-guard) {
406: BVDestroy(&Y);
407: }
408: return(0);
409: }
411: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
412: {
413: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
416: if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
417: if (bs <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %D",bs);
418: if (ctx->bs != bs) {
419: ctx->bs = bs;
420: eps->state = EPS_STATE_INITIAL;
421: }
422: return(0);
423: }
425: /*@
426: EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.
428: Logically Collective on EPS
430: Input Parameters:
431: + eps - the eigenproblem solver context
432: - bs - the block size
434: Options Database Key:
435: . -eps_lobpcg_blocksize - Sets the block size
437: Level: advanced
439: .seealso: EPSLOBPCGGetBlockSize()
440: @*/
441: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
442: {
448: PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
449: return(0);
450: }
452: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
453: {
454: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
457: *bs = ctx->bs;
458: return(0);
459: }
461: /*@
462: EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.
464: Not Collective
466: Input Parameter:
467: . eps - the eigenproblem solver context
469: Output Parameter:
470: . bs - the block size
472: Level: advanced
474: .seealso: EPSLOBPCGSetBlockSize()
475: @*/
476: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
477: {
483: PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
484: return(0);
485: }
487: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
488: {
489: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
492: if (restart==PETSC_DEFAULT) ctx->restart = 0.6;
493: else {
494: 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);
495: ctx->restart = restart;
496: }
497: return(0);
498: }
500: /*@
501: EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
502: The meaning of this parameter is the proportion of vectors within the
503: current block iterate that must have converged in order to force a
504: restart with hard locking.
506: Logically Collective on EPS
508: Input Parameters:
509: + eps - the eigenproblem solver context
510: - restart - the percentage of the block of vectors to force a restart
512: Options Database Key:
513: . -eps_lobpcg_restart - Sets the restart parameter
515: Notes:
516: Allowed values are in the range [0.1,1.0]. The default is 0.6.
518: Level: advanced
520: .seealso: EPSLOBPCGGetRestart()
521: @*/
522: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
523: {
529: PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
530: return(0);
531: }
533: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
534: {
535: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
538: *restart = ctx->restart;
539: return(0);
540: }
542: /*@
543: EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.
545: Not Collective
547: Input Parameter:
548: . eps - the eigenproblem solver context
550: Output Parameter:
551: . restart - the restart parameter
553: Level: advanced
555: .seealso: EPSLOBPCGSetRestart()
556: @*/
557: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
558: {
564: PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
565: return(0);
566: }
568: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
569: {
570: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
573: ctx->lock = lock;
574: return(0);
575: }
577: /*@
578: EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
579: the LOBPCG method.
581: Logically Collective on EPS
583: Input Parameters:
584: + eps - the eigenproblem solver context
585: - lock - true if the locking variant must be selected
587: Options Database Key:
588: . -eps_lobpcg_locking - Sets the locking flag
590: Notes:
591: This flag refers to soft locking (converged vectors within the current
592: block iterate), since hard locking is always used (when nev is larger
593: than the block size).
595: Level: advanced
597: .seealso: EPSLOBPCGGetLocking()
598: @*/
599: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
600: {
606: PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
607: return(0);
608: }
610: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
611: {
612: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
615: *lock = ctx->lock;
616: return(0);
617: }
619: /*@
620: EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.
622: Not Collective
624: Input Parameter:
625: . eps - the eigenproblem solver context
627: Output Parameter:
628: . lock - the locking flag
630: Level: advanced
632: .seealso: EPSLOBPCGSetLocking()
633: @*/
634: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
635: {
641: PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
642: return(0);
643: }
645: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
646: {
648: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
649: PetscBool isascii;
652: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
653: if (isascii) {
654: PetscViewerASCIIPrintf(viewer," block size %D\n",ctx->bs);
655: PetscViewerASCIIPrintf(viewer," restart parameter=%g (using %d guard vectors)\n",(double)ctx->restart,(int)PetscRoundReal((1.0-ctx->restart)*ctx->bs));
656: PetscViewerASCIIPrintf(viewer," soft locking %sactivated\n",ctx->lock?"":"de");
657: }
658: return(0);
659: }
661: PetscErrorCode EPSSetFromOptions_LOBPCG(PetscOptionItems *PetscOptionsObject,EPS eps)
662: {
664: PetscBool lock,flg;
665: PetscInt bs;
666: PetscReal restart;
669: PetscOptionsHead(PetscOptionsObject,"EPS LOBPCG Options");
671: PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg);
672: if (flg) { EPSLOBPCGSetBlockSize(eps,bs); }
674: PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg);
675: if (flg) { EPSLOBPCGSetRestart(eps,restart); }
677: PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg);
678: if (flg) { EPSLOBPCGSetLocking(eps,lock); }
680: PetscOptionsTail();
681: return(0);
682: }
684: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
685: {
689: PetscFree(eps->data);
690: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
691: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
692: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
693: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
694: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
695: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
696: return(0);
697: }
699: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
700: {
701: EPS_LOBPCG *lobpcg;
705: PetscNewLog(eps,&lobpcg);
706: eps->data = (void*)lobpcg;
707: lobpcg->lock = PETSC_TRUE;
709: eps->useds = PETSC_TRUE;
710: eps->categ = EPS_CATEGORY_PRECOND;
712: eps->ops->solve = EPSSolve_LOBPCG;
713: eps->ops->setup = EPSSetUp_LOBPCG;
714: eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
715: eps->ops->destroy = EPSDestroy_LOBPCG;
716: eps->ops->view = EPSView_LOBPCG;
717: eps->ops->backtransform = EPSBackTransform_Default;
718: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
720: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
721: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
722: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
723: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
724: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
725: PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
726: return(0);
727: }