Actual source code: lobpcg.c
slepc-main 2024-12-17
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: static PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
41: {
42: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
43: PetscInt k;
45: PetscFunctionBegin;
46: if (*nev==0) *nev = 1;
47: k = PetscMax(3*ctx->bs,((*nev-1)/ctx->bs+3)*ctx->bs);
48: if (*ncv!=PETSC_DETERMINE) { /* ncv set */
49: PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
50: } else *ncv = k;
51: if (*mpd==PETSC_DETERMINE) *mpd = 3*ctx->bs;
52: else PetscCheck(*mpd==3*ctx->bs,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"This solver does not allow a value of mpd different from 3*blocksize");
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
57: {
58: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
60: PetscFunctionBegin;
61: EPSCheckHermitianDefinite(eps);
62: EPSCheckNotStructured(eps);
63: if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev?eps->nev:1);
64: PetscCheck(eps->n-eps->nds>=5*ctx->bs,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The problem size is too small relative to the block size");
65: PetscCall(EPSSetDimensions_LOBPCG(eps,&eps->nev,&eps->ncv,&eps->mpd));
66: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
67: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
68: PetscCheck(eps->which==EPS_SMALLEST_REAL || eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real or largest real eigenvalues");
69: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD);
70: EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
72: if (!ctx->restart) ctx->restart = 0.9;
74: /* number of guard vectors */
75: if (ctx->bs==1) ctx->guard = 0;
76: else ctx->guard = PetscMin((PetscInt)((1.0-ctx->restart)*ctx->bs+0.45),ctx->bs-1);
78: PetscCall(EPSAllocateSolution(eps,0));
79: PetscCall(EPS_SetInnerProduct(eps));
80: PetscCall(DSSetType(eps->ds,DSGHEP));
81: PetscCall(DSAllocate(eps->ds,eps->mpd));
82: PetscCall(EPSSetWorkVecs(eps,1));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode EPSSolve_LOBPCG(EPS eps)
87: {
88: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
89: PetscInt i,j,k,nv,ini,nmat,nc,nconv,locked,its,prev=0;
90: PetscReal norm;
91: PetscScalar *eigr,dot;
92: PetscBool breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
93: Mat A,B,M,V=NULL,W=NULL;
94: Vec v,z,w=eps->work[0];
95: BV X,Y=NULL,Z,R,P,AX,BX;
96: SlepcSC sc;
98: PetscFunctionBegin;
99: PetscCall(STGetNumMatrices(eps->st,&nmat));
100: PetscCall(STGetMatrix(eps->st,0,&A));
101: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
102: else B = NULL;
104: if (eps->which==EPS_LARGEST_REAL) { /* flip spectrum */
105: flip = PETSC_TRUE;
106: PetscCall(DSGetSlepcSC(eps->ds,&sc));
107: sc->comparison = SlepcCompareSmallestReal;
108: }
110: /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
111: PetscCall(PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL));
113: /* 1. Allocate memory */
114: PetscCall(PetscCalloc1(3*ctx->bs,&eigr));
115: PetscCall(BVDuplicateResize(eps->V,3*ctx->bs,&Z));
116: PetscCall(BVDuplicateResize(eps->V,ctx->bs,&X));
117: PetscCall(BVDuplicateResize(eps->V,ctx->bs,&R));
118: PetscCall(BVDuplicateResize(eps->V,ctx->bs,&P));
119: PetscCall(BVDuplicateResize(eps->V,ctx->bs,&AX));
120: if (B) PetscCall(BVDuplicateResize(eps->V,ctx->bs,&BX));
121: nc = eps->nds;
122: if (nc>0 || eps->nev>ctx->bs-ctx->guard) PetscCall(BVDuplicateResize(eps->V,nc+eps->nev,&Y));
123: if (nc>0) {
124: for (j=0;j<nc;j++) {
125: PetscCall(BVGetColumn(eps->V,-nc+j,&v));
126: PetscCall(BVInsertVec(Y,j,v));
127: PetscCall(BVRestoreColumn(eps->V,-nc+j,&v));
128: }
129: PetscCall(BVSetActiveColumns(Y,0,nc));
130: }
132: /* 2. Apply the constraints to the initial vectors */
133: /* 3. B-orthogonalize initial vectors */
134: for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
135: PetscCall(BVSetRandomColumn(eps->V,k));
136: PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
137: }
138: nv = ctx->bs;
139: PetscCall(BVSetActiveColumns(eps->V,0,nv));
140: PetscCall(BVSetActiveColumns(Z,0,nv));
141: PetscCall(BVCopy(eps->V,Z));
142: PetscCall(BVCopy(Z,X));
144: /* 4. Compute initial Ritz vectors */
145: PetscCall(BVMatMult(X,A,AX));
146: PetscCall(DSSetDimensions(eps->ds,nv,0,0));
147: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
148: PetscCall(BVMatProject(AX,NULL,X,M));
149: if (flip) PetscCall(MatScale(M,-1.0));
150: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
151: PetscCall(DSSetIdentity(eps->ds,DS_MAT_B));
152: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
153: PetscCall(DSSolve(eps->ds,eigr,NULL));
154: PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
155: PetscCall(DSSynchronize(eps->ds,eigr,NULL));
156: for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
157: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
158: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
159: PetscCall(BVMultInPlace(X,M,0,nv));
160: PetscCall(BVMultInPlace(AX,M,0,nv));
161: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
163: /* 5. Initialize range of active iterates */
164: locked = 0; /* hard-locked vectors, the leading locked columns of V are eigenvectors */
165: nconv = 0; /* number of converged eigenvalues in the current block */
166: its = 0; /* iterations for the current block */
168: /* 6. Main loop */
169: while (eps->reason == EPS_CONVERGED_ITERATING) {
171: if (ctx->lock) {
172: PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
173: PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
174: if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
175: }
177: /* 7. Compute residuals */
178: ini = (ctx->lock)? nconv: 0;
179: PetscCall(BVCopy(AX,R));
180: if (B) PetscCall(BVMatMult(X,B,BX));
181: for (j=ini;j<ctx->bs;j++) {
182: PetscCall(BVGetColumn(R,j,&v));
183: PetscCall(BVGetColumn(B?BX:X,j,&z));
184: PetscCall(VecAXPY(v,-eps->eigr[locked+j],z));
185: PetscCall(BVRestoreColumn(R,j,&v));
186: PetscCall(BVRestoreColumn(B?BX:X,j,&z));
187: }
189: /* 8. Compute residual norms and update index set of active iterates */
190: k = ini;
191: countc = PETSC_TRUE;
192: for (j=ini;j<ctx->bs;j++) {
193: i = locked+j;
194: PetscCall(BVGetColumn(R,j,&v));
195: PetscCall(VecNorm(v,NORM_2,&norm));
196: PetscCall(BVRestoreColumn(R,j,&v));
197: PetscCall((*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx));
198: if (countc) {
199: if (eps->errest[i] < eps->tol) k++;
200: else countc = PETSC_FALSE;
201: }
202: if (!countc && !eps->trackall) break;
203: }
204: nconv = k;
205: eps->nconv = locked + nconv;
206: if (its) PetscCall(EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs));
207: PetscCall((*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
208: if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
209: PetscCall(BVSetActiveColumns(eps->V,locked,eps->nconv));
210: PetscCall(BVSetActiveColumns(X,0,nconv));
211: PetscCall(BVCopy(X,eps->V));
212: }
213: if (eps->reason != EPS_CONVERGED_ITERATING) {
214: break;
215: } else if (nconv >= ctx->bs-ctx->guard) {
216: eps->its += its-1;
217: its = 0;
218: } else its++;
220: if (nconv >= ctx->bs-ctx->guard) { /* force hard locking of vectors and compute new R */
222: /* extend constraints */
223: PetscCall(BVSetActiveColumns(Y,nc+locked,nc+locked+nconv));
224: PetscCall(BVCopy(X,Y));
225: PetscCall(BVSetActiveColumns(Y,0,nc+locked+nconv));
227: /* shift work BV's */
228: for (j=nconv;j<ctx->bs;j++) {
229: PetscCall(BVCopyColumn(X,j,j-nconv));
230: PetscCall(BVCopyColumn(R,j,j-nconv));
231: PetscCall(BVCopyColumn(P,j,j-nconv));
232: PetscCall(BVCopyColumn(AX,j,j-nconv));
233: if (B) PetscCall(BVCopyColumn(BX,j,j-nconv));
234: }
236: /* set new initial vectors */
237: PetscCall(BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv));
238: PetscCall(BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs));
239: PetscCall(BVCopy(eps->V,X));
240: for (j=ctx->bs-nconv;j<ctx->bs;j++) {
241: PetscCall(BVGetColumn(X,j,&v));
242: PetscCall(BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown));
243: if (norm>0.0 && !breakdown) PetscCall(VecScale(v,1.0/norm));
244: else {
245: PetscCall(PetscInfo(eps,"Orthogonalization of initial vector failed\n"));
246: eps->reason = EPS_DIVERGED_BREAKDOWN;
247: goto diverged;
248: }
249: PetscCall(BVRestoreColumn(X,j,&v));
250: }
251: locked += nconv;
252: nconv = 0;
253: PetscCall(BVSetActiveColumns(X,nconv,ctx->bs));
255: /* B-orthogonalize initial vectors */
256: PetscCall(BVOrthogonalize(X,NULL));
257: PetscCall(BVSetActiveColumns(Z,nconv,ctx->bs));
258: PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
259: PetscCall(BVCopy(X,Z));
261: /* compute initial Ritz vectors */
262: nv = ctx->bs;
263: PetscCall(BVMatMult(X,A,AX));
264: PetscCall(DSSetDimensions(eps->ds,nv,0,0));
265: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
266: PetscCall(BVMatProject(AX,NULL,X,M));
267: if (flip) PetscCall(MatScale(M,-1.0));
268: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
269: PetscCall(DSSetIdentity(eps->ds,DS_MAT_B));
270: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
271: PetscCall(DSSolve(eps->ds,eigr,NULL));
272: PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
273: PetscCall(DSSynchronize(eps->ds,eigr,NULL));
274: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
275: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
276: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
277: PetscCall(BVMultInPlace(X,M,0,nv));
278: PetscCall(BVMultInPlace(AX,M,0,nv));
279: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
281: continue; /* skip the rest of the iteration */
282: }
284: ini = (ctx->lock)? nconv: 0;
285: if (ctx->lock) {
286: PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
287: PetscCall(BVSetActiveColumns(P,nconv,ctx->bs));
288: PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
289: if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
290: }
292: /* 9. Apply preconditioner to the residuals */
293: PetscCall(BVGetMat(R,&V));
294: if (prev != ctx->bs-ini) {
295: prev = ctx->bs-ini;
296: PetscCall(MatDestroy(&W));
297: PetscCall(MatDuplicate(V,MAT_SHARE_NONZERO_PATTERN,&W));
298: }
299: PetscCall(STApplyMat(eps->st,V,W));
300: if (checkprecond) {
301: for (j=ini;j<ctx->bs;j++) {
302: PetscCall(MatDenseGetColumnVecRead(V,j-ini,&v));
303: PetscCall(MatDenseGetColumnVecRead(W,j-ini,&w));
304: PetscCall(VecDot(v,w,&dot));
305: PetscCall(MatDenseRestoreColumnVecRead(W,j-ini,&w));
306: PetscCall(MatDenseRestoreColumnVecRead(V,j-ini,&v));
307: if (PetscRealPart(dot)<0.0) {
308: PetscCall(PetscInfo(eps,"The preconditioner is not positive-definite\n"));
309: eps->reason = EPS_DIVERGED_BREAKDOWN;
310: goto diverged;
311: }
312: }
313: }
314: if (nc+locked>0) {
315: for (j=ini;j<ctx->bs;j++) {
316: PetscCall(MatDenseGetColumnVecWrite(W,j-ini,&w));
317: PetscCall(BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown));
318: if (norm>0.0 && !breakdown) PetscCall(VecScale(w,1.0/norm));
319: PetscCall(MatDenseRestoreColumnVecWrite(W,j-ini,&w));
320: if (norm<=0.0 || breakdown) {
321: PetscCall(PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n"));
322: eps->reason = EPS_DIVERGED_BREAKDOWN;
323: goto diverged;
324: }
325: }
326: }
327: PetscCall(MatCopy(W,V,SAME_NONZERO_PATTERN));
328: PetscCall(BVRestoreMat(R,&V));
330: /* 11. B-orthonormalize preconditioned residuals */
331: PetscCall(BVOrthogonalize(R,NULL));
333: /* 13-16. B-orthonormalize conjugate directions */
334: if (its>1) PetscCall(BVOrthogonalize(P,NULL));
336: /* 17-23. Compute symmetric Gram matrices */
337: PetscCall(BVSetActiveColumns(Z,0,ctx->bs));
338: PetscCall(BVSetActiveColumns(X,0,ctx->bs));
339: PetscCall(BVCopy(X,Z));
340: PetscCall(BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini));
341: PetscCall(BVCopy(R,Z));
342: if (its>1) {
343: PetscCall(BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini));
344: PetscCall(BVCopy(P,Z));
345: }
347: if (its>1) nv = 3*ctx->bs-2*ini;
348: else nv = 2*ctx->bs-ini;
350: PetscCall(BVSetActiveColumns(Z,0,nv));
351: PetscCall(DSSetDimensions(eps->ds,nv,0,0));
352: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
353: PetscCall(BVMatProject(Z,A,Z,M));
354: if (flip) PetscCall(MatScale(M,-1.0));
355: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
356: PetscCall(DSGetMat(eps->ds,DS_MAT_B,&M));
357: PetscCall(BVMatProject(Z,B,Z,M)); /* covers also the case B=NULL */
358: PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&M));
360: /* 24. Solve the generalized eigenvalue problem */
361: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
362: PetscCall(DSSolve(eps->ds,eigr,NULL));
363: PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
364: PetscCall(DSSynchronize(eps->ds,eigr,NULL));
365: for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
366: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
368: /* 25-33. Compute Ritz vectors */
369: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
370: PetscCall(BVSetActiveColumns(Z,ctx->bs,nv));
371: if (ctx->lock) PetscCall(BVSetActiveColumns(P,0,ctx->bs));
372: PetscCall(BVMult(P,1.0,0.0,Z,M));
373: PetscCall(BVCopy(P,X));
374: if (ctx->lock) PetscCall(BVSetActiveColumns(P,nconv,ctx->bs));
375: PetscCall(BVSetActiveColumns(Z,0,ctx->bs));
376: PetscCall(BVMult(X,1.0,1.0,Z,M));
377: if (ctx->lock) PetscCall(BVSetActiveColumns(X,nconv,ctx->bs));
378: PetscCall(BVMatMult(X,A,AX));
379: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
380: }
382: diverged:
383: eps->its += its;
385: if (flip) sc->comparison = SlepcCompareLargestReal;
386: PetscCall(PetscFree(eigr));
387: PetscCall(MatDestroy(&W));
388: if (V) PetscCall(BVRestoreMat(R,&V)); /* only needed when goto diverged is reached */
389: PetscCall(BVDestroy(&Z));
390: PetscCall(BVDestroy(&X));
391: PetscCall(BVDestroy(&R));
392: PetscCall(BVDestroy(&P));
393: PetscCall(BVDestroy(&AX));
394: if (B) PetscCall(BVDestroy(&BX));
395: if (nc>0 || eps->nev>ctx->bs-ctx->guard) PetscCall(BVDestroy(&Y));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
400: {
401: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
403: PetscFunctionBegin;
404: if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 0;
405: else PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %" PetscInt_FMT,bs);
406: if (ctx->bs != bs) {
407: ctx->bs = bs;
408: eps->state = EPS_STATE_INITIAL;
409: }
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.
416: Logically Collective
418: Input Parameters:
419: + eps - the eigenproblem solver context
420: - bs - the block size
422: Options Database Key:
423: . -eps_lobpcg_blocksize - Sets the block size
425: Level: advanced
427: .seealso: EPSLOBPCGGetBlockSize()
428: @*/
429: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
430: {
431: PetscFunctionBegin;
434: PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
439: {
440: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
442: PetscFunctionBegin;
443: *bs = ctx->bs;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*@
448: EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.
450: Not Collective
452: Input Parameter:
453: . eps - the eigenproblem solver context
455: Output Parameter:
456: . bs - the block size
458: Level: advanced
460: .seealso: EPSLOBPCGSetBlockSize()
461: @*/
462: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
463: {
464: PetscFunctionBegin;
466: PetscAssertPointer(bs,2);
467: PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
472: {
473: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
475: PetscFunctionBegin;
476: if (restart==(PetscReal)PETSC_DEFAULT || restart==(PetscReal)PETSC_DECIDE) restart = 0.9;
477: PetscCheck(restart>=0.1 && restart<=1.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The restart argument %g must be in the range [0.1,1.0]",(double)restart);
478: if (restart != ctx->restart) {
479: ctx->restart = restart;
480: eps->state = EPS_STATE_INITIAL;
481: }
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: /*@
486: EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
488: Logically Collective
490: Input Parameters:
491: + eps - the eigenproblem solver context
492: - restart - the percentage of the block of vectors to force a restart
494: Options Database Key:
495: . -eps_lobpcg_restart - Sets the restart parameter
497: Notes:
498: The meaning of this parameter is the proportion of vectors within the
499: current block iterate that must have converged in order to force a
500: restart with hard locking.
501: Allowed values are in the range [0.1,1.0]. The default is 0.9.
503: Level: advanced
505: .seealso: EPSLOBPCGGetRestart()
506: @*/
507: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
508: {
509: PetscFunctionBegin;
512: PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
517: {
518: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
520: PetscFunctionBegin;
521: *restart = ctx->restart;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*@
526: EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.
528: Not Collective
530: Input Parameter:
531: . eps - the eigenproblem solver context
533: Output Parameter:
534: . restart - the restart parameter
536: Level: advanced
538: .seealso: EPSLOBPCGSetRestart()
539: @*/
540: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
541: {
542: PetscFunctionBegin;
544: PetscAssertPointer(restart,2);
545: PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
550: {
551: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
553: PetscFunctionBegin;
554: ctx->lock = lock;
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
560: the LOBPCG method.
562: Logically Collective
564: Input Parameters:
565: + eps - the eigenproblem solver context
566: - lock - true if the locking variant must be selected
568: Options Database Key:
569: . -eps_lobpcg_locking - Sets the locking flag
571: Notes:
572: This flag refers to soft locking (converged vectors within the current
573: block iterate), since hard locking is always used (when nev is larger
574: than the block size).
576: Level: advanced
578: .seealso: EPSLOBPCGGetLocking()
579: @*/
580: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
581: {
582: PetscFunctionBegin;
585: PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
590: {
591: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
593: PetscFunctionBegin;
594: *lock = ctx->lock;
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: /*@
599: EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.
601: Not Collective
603: Input Parameter:
604: . eps - the eigenproblem solver context
606: Output Parameter:
607: . lock - the locking flag
609: Level: advanced
611: .seealso: EPSLOBPCGSetLocking()
612: @*/
613: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
614: {
615: PetscFunctionBegin;
617: PetscAssertPointer(lock,2);
618: PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: static PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
623: {
624: EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
625: PetscBool isascii;
627: PetscFunctionBegin;
628: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
629: if (isascii) {
630: PetscCall(PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs));
631: PetscCall(PetscViewerASCIIPrintf(viewer," restart parameter=%g (using %" PetscInt_FMT " guard vectors)\n",(double)ctx->restart,ctx->guard));
632: PetscCall(PetscViewerASCIIPrintf(viewer," soft locking %sactivated\n",ctx->lock?"":"de"));
633: }
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: static PetscErrorCode EPSSetFromOptions_LOBPCG(EPS eps,PetscOptionItems *PetscOptionsObject)
638: {
639: PetscBool lock,flg;
640: PetscInt bs;
641: PetscReal restart;
643: PetscFunctionBegin;
644: PetscOptionsHeadBegin(PetscOptionsObject,"EPS LOBPCG Options");
646: PetscCall(PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg));
647: if (flg) PetscCall(EPSLOBPCGSetBlockSize(eps,bs));
649: PetscCall(PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg));
650: if (flg) PetscCall(EPSLOBPCGSetRestart(eps,restart));
652: PetscCall(PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg));
653: if (flg) PetscCall(EPSLOBPCGSetLocking(eps,lock));
655: PetscOptionsHeadEnd();
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
660: {
661: PetscFunctionBegin;
662: PetscCall(PetscFree(eps->data));
663: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL));
664: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL));
665: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL));
666: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL));
667: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL));
668: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL));
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
673: {
674: EPS_LOBPCG *lobpcg;
676: PetscFunctionBegin;
677: PetscCall(PetscNew(&lobpcg));
678: eps->data = (void*)lobpcg;
679: lobpcg->lock = PETSC_TRUE;
681: eps->useds = PETSC_TRUE;
682: eps->categ = EPS_CATEGORY_PRECOND;
684: eps->ops->solve = EPSSolve_LOBPCG;
685: eps->ops->setup = EPSSetUp_LOBPCG;
686: eps->ops->setupsort = EPSSetUpSort_Default;
687: eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
688: eps->ops->destroy = EPSDestroy_LOBPCG;
689: eps->ops->view = EPSView_LOBPCG;
690: eps->ops->backtransform = EPSBackTransform_Default;
691: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
693: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG));
694: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG));
695: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG));
696: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG));
697: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG));
698: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }