Actual source code: lobpcg.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  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:   return 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:   return 0;
 80: }

 82: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
 83: {
 84:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
 85:   PetscInt       i,j,k,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:   STGetNumMatrices(eps->st,&nmat);
 95:   STGetMatrix(eps->st,0,&A);
 96:   if (nmat>1) STGetMatrix(eps->st,1,&B);
 97:   else B = NULL;

 99:   if (eps->which==EPS_LARGEST_REAL) {  /* flip spectrum */
100:     flip = PETSC_TRUE;
101:     DSGetSlepcSC(eps->ds,&sc);
102:     sc->comparison = SlepcCompareSmallestReal;
103:   }

105:   /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
106:   PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL);

108:   /* 1. Allocate memory */
109:   PetscCalloc1(3*ctx->bs,&eigr);
110:   BVDuplicateResize(eps->V,3*ctx->bs,&Z);
111:   BVDuplicateResize(eps->V,ctx->bs,&X);
112:   BVDuplicateResize(eps->V,ctx->bs,&R);
113:   BVDuplicateResize(eps->V,ctx->bs,&P);
114:   BVDuplicateResize(eps->V,ctx->bs,&AX);
115:   if (B) BVDuplicateResize(eps->V,ctx->bs,&BX);
116:   nc = eps->nds;
117:   if (nc>0 || eps->nev>ctx->bs-ctx->guard) BVDuplicateResize(eps->V,nc+eps->nev,&Y);
118:   if (nc>0) {
119:     for (j=0;j<nc;j++) {
120:       BVGetColumn(eps->V,-nc+j,&v);
121:       BVInsertVec(Y,j,v);
122:       BVRestoreColumn(eps->V,-nc+j,&v);
123:     }
124:     BVSetActiveColumns(Y,0,nc);
125:   }

127:   /* 2. Apply the constraints to the initial vectors */
128:   /* 3. B-orthogonalize initial vectors */
129:   for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
130:     BVSetRandomColumn(eps->V,k);
131:     BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
132:   }
133:   nv = ctx->bs;
134:   BVSetActiveColumns(eps->V,0,nv);
135:   BVSetActiveColumns(Z,0,nv);
136:   BVCopy(eps->V,Z);
137:   BVCopy(Z,X);

139:   /* 4. Compute initial Ritz vectors */
140:   BVMatMult(X,A,AX);
141:   DSSetDimensions(eps->ds,nv,0,0);
142:   DSGetMat(eps->ds,DS_MAT_A,&M);
143:   BVMatProject(AX,NULL,X,M);
144:   if (flip) MatScale(M,-1.0);
145:   DSRestoreMat(eps->ds,DS_MAT_A,&M);
146:   DSSetIdentity(eps->ds,DS_MAT_B);
147:   DSSetState(eps->ds,DS_STATE_RAW);
148:   DSSolve(eps->ds,eigr,NULL);
149:   DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
150:   DSSynchronize(eps->ds,eigr,NULL);
151:   for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
152:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
153:   DSGetMat(eps->ds,DS_MAT_X,&M);
154:   BVMultInPlace(X,M,0,nv);
155:   BVMultInPlace(AX,M,0,nv);
156:   DSRestoreMat(eps->ds,DS_MAT_X,&M);

158:   /* 5. Initialize range of active iterates */
159:   locked = 0;  /* hard-locked vectors, the leading locked columns of V are eigenvectors */
160:   nconv  = 0;  /* number of converged eigenvalues in the current block */
161:   its    = 0;  /* iterations for the current block */

163:   /* 6. Main loop */
164:   while (eps->reason == EPS_CONVERGED_ITERATING) {

166:     if (ctx->lock) {
167:       BVSetActiveColumns(R,nconv,ctx->bs);
168:       BVSetActiveColumns(AX,nconv,ctx->bs);
169:       if (B) BVSetActiveColumns(BX,nconv,ctx->bs);
170:     }

172:     /* 7. Compute residuals */
173:     ini = (ctx->lock)? nconv: 0;
174:     BVCopy(AX,R);
175:     if (B) BVMatMult(X,B,BX);
176:     for (j=ini;j<ctx->bs;j++) {
177:       BVGetColumn(R,j,&v);
178:       BVGetColumn(B?BX:X,j,&z);
179:       VecAXPY(v,-eps->eigr[locked+j],z);
180:       BVRestoreColumn(R,j,&v);
181:       BVRestoreColumn(B?BX:X,j,&z);
182:     }

184:     /* 8. Compute residual norms and update index set of active iterates */
185:     k = ini;
186:     countc = PETSC_TRUE;
187:     for (j=ini;j<ctx->bs;j++) {
188:       i = locked+j;
189:       BVGetColumn(R,j,&v);
190:       VecNorm(v,NORM_2,&norm);
191:       BVRestoreColumn(R,j,&v);
192:       (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx);
193:       if (countc) {
194:         if (eps->errest[i] < eps->tol) k++;
195:         else countc = PETSC_FALSE;
196:       }
197:       if (!countc && !eps->trackall) break;
198:     }
199:     nconv = k;
200:     eps->nconv = locked + nconv;
201:     if (its) EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs);
202:     (*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
203:     if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
204:       BVSetActiveColumns(eps->V,locked,eps->nconv);
205:       BVSetActiveColumns(X,0,nconv);
206:       BVCopy(X,eps->V);
207:     }
208:     if (eps->reason != EPS_CONVERGED_ITERATING) {
209:       break;
210:     } else if (nconv >= ctx->bs-ctx->guard) {
211:       eps->its += its-1;
212:       its = 0;
213:     } else its++;

215:     if (nconv >= ctx->bs-ctx->guard) {  /* force hard locking of vectors and compute new R */

217:       /* extend constraints */
218:       BVSetActiveColumns(Y,nc+locked,nc+locked+nconv);
219:       BVCopy(X,Y);
220:       BVSetActiveColumns(Y,0,nc+locked+nconv);

222:       /* shift work BV's */
223:       for (j=nconv;j<ctx->bs;j++) {
224:         BVCopyColumn(X,j,j-nconv);
225:         BVCopyColumn(R,j,j-nconv);
226:         BVCopyColumn(P,j,j-nconv);
227:         BVCopyColumn(AX,j,j-nconv);
228:         if (B) BVCopyColumn(BX,j,j-nconv);
229:       }

231:       /* set new initial vectors */
232:       BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv);
233:       BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs);
234:       BVCopy(eps->V,X);
235:       for (j=ctx->bs-nconv;j<ctx->bs;j++) {
236:         BVGetColumn(X,j,&v);
237:         BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown);
238:         if (norm>0.0 && !breakdown) VecScale(v,1.0/norm);
239:         else {
240:           PetscInfo(eps,"Orthogonalization of initial vector failed\n");
241:           eps->reason = EPS_DIVERGED_BREAKDOWN;
242:           goto diverged;
243:         }
244:         BVRestoreColumn(X,j,&v);
245:       }
246:       locked += nconv;
247:       nconv = 0;
248:       BVSetActiveColumns(X,nconv,ctx->bs);

250:       /* B-orthogonalize initial vectors */
251:       BVOrthogonalize(X,NULL);
252:       BVSetActiveColumns(Z,nconv,ctx->bs);
253:       BVSetActiveColumns(AX,nconv,ctx->bs);
254:       BVCopy(X,Z);

256:       /* compute initial Ritz vectors */
257:       nv = ctx->bs;
258:       BVMatMult(X,A,AX);
259:       DSSetDimensions(eps->ds,nv,0,0);
260:       DSGetMat(eps->ds,DS_MAT_A,&M);
261:       BVMatProject(AX,NULL,X,M);
262:       if (flip) MatScale(M,-1.0);
263:       DSRestoreMat(eps->ds,DS_MAT_A,&M);
264:       DSSetIdentity(eps->ds,DS_MAT_B);
265:       DSSetState(eps->ds,DS_STATE_RAW);
266:       DSSolve(eps->ds,eigr,NULL);
267:       DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
268:       DSSynchronize(eps->ds,eigr,NULL);
269:       for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
270:       DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
271:       DSGetMat(eps->ds,DS_MAT_X,&M);
272:       BVMultInPlace(X,M,0,nv);
273:       BVMultInPlace(AX,M,0,nv);
274:       DSRestoreMat(eps->ds,DS_MAT_X,&M);

276:       continue;   /* skip the rest of the iteration */
277:     }

279:     ini = (ctx->lock)? nconv: 0;
280:     if (ctx->lock) {
281:       BVSetActiveColumns(R,nconv,ctx->bs);
282:       BVSetActiveColumns(P,nconv,ctx->bs);
283:       BVSetActiveColumns(AX,nconv,ctx->bs);
284:       if (B) BVSetActiveColumns(BX,nconv,ctx->bs);
285:     }

287:     /* 9. Apply preconditioner to the residuals */
288:     BVGetMat(R,&V);
289:     if (prev != ctx->bs-ini) {
290:       prev = ctx->bs-ini;
291:       MatDestroy(&W);
292:       MatDuplicate(V,MAT_SHARE_NONZERO_PATTERN,&W);
293:     }
294:     STApplyMat(eps->st,V,W);
295:     if (checkprecond) {
296:       for (j=ini;j<ctx->bs;j++) {
297:         MatDenseGetColumnVecRead(V,j-ini,&v);
298:         MatDenseGetColumnVecRead(W,j-ini,&w);
299:         VecDot(v,w,&dot);
300:         MatDenseRestoreColumnVecRead(W,j-ini,&w);
301:         MatDenseRestoreColumnVecRead(V,j-ini,&v);
302:         if (PetscRealPart(dot)<0.0) {
303:           PetscInfo(eps,"The preconditioner is not positive-definite\n");
304:           eps->reason = EPS_DIVERGED_BREAKDOWN;
305:           goto diverged;
306:         }
307:       }
308:     }
309:     if (nc+locked>0) {
310:       for (j=ini;j<ctx->bs;j++) {
311:         MatDenseGetColumnVecWrite(W,j-ini,&w);
312:         BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
313:         if (norm>0.0 && !breakdown) VecScale(w,1.0/norm);
314:         MatDenseRestoreColumnVecWrite(W,j-ini,&w);
315:         if (norm<=0.0 || breakdown) {
316:           PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
317:           eps->reason = EPS_DIVERGED_BREAKDOWN;
318:           goto diverged;
319:         }
320:       }
321:     }
322:     MatCopy(W,V,SAME_NONZERO_PATTERN);
323:     BVRestoreMat(R,&V);

325:     /* 11. B-orthonormalize preconditioned residuals */
326:     BVOrthogonalize(R,NULL);

328:     /* 13-16. B-orthonormalize conjugate directions */
329:     if (its>1) BVOrthogonalize(P,NULL);

331:     /* 17-23. Compute symmetric Gram matrices */
332:     BVSetActiveColumns(Z,0,ctx->bs);
333:     BVSetActiveColumns(X,0,ctx->bs);
334:     BVCopy(X,Z);
335:     BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
336:     BVCopy(R,Z);
337:     if (its>1) {
338:       BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
339:       BVCopy(P,Z);
340:     }

342:     if (its>1) nv = 3*ctx->bs-2*ini;
343:     else nv = 2*ctx->bs-ini;

345:     BVSetActiveColumns(Z,0,nv);
346:     DSSetDimensions(eps->ds,nv,0,0);
347:     DSGetMat(eps->ds,DS_MAT_A,&M);
348:     BVMatProject(Z,A,Z,M);
349:     if (flip) MatScale(M,-1.0);
350:     DSRestoreMat(eps->ds,DS_MAT_A,&M);
351:     DSGetMat(eps->ds,DS_MAT_B,&M);
352:     BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
353:     DSRestoreMat(eps->ds,DS_MAT_B,&M);

355:     /* 24. Solve the generalized eigenvalue problem */
356:     DSSetState(eps->ds,DS_STATE_RAW);
357:     DSSolve(eps->ds,eigr,NULL);
358:     DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
359:     DSSynchronize(eps->ds,eigr,NULL);
360:     for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
361:     DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

363:     /* 25-33. Compute Ritz vectors */
364:     DSGetMat(eps->ds,DS_MAT_X,&M);
365:     BVSetActiveColumns(Z,ctx->bs,nv);
366:     if (ctx->lock) BVSetActiveColumns(P,0,ctx->bs);
367:     BVMult(P,1.0,0.0,Z,M);
368:     BVCopy(P,X);
369:     if (ctx->lock) BVSetActiveColumns(P,nconv,ctx->bs);
370:     BVSetActiveColumns(Z,0,ctx->bs);
371:     BVMult(X,1.0,1.0,Z,M);
372:     if (ctx->lock) BVSetActiveColumns(X,nconv,ctx->bs);
373:     BVMatMult(X,A,AX);
374:     DSRestoreMat(eps->ds,DS_MAT_X,&M);
375:   }

377: diverged:
378:   eps->its += its;

380:   if (flip) sc->comparison = SlepcCompareLargestReal;
381:   PetscFree(eigr);
382:   MatDestroy(&W);
383:   if (V) BVRestoreMat(R,&V); /* only needed when goto diverged is reached */
384:   BVDestroy(&Z);
385:   BVDestroy(&X);
386:   BVDestroy(&R);
387:   BVDestroy(&P);
388:   BVDestroy(&AX);
389:   if (B) BVDestroy(&BX);
390:   if (nc>0 || eps->nev>ctx->bs-ctx->guard) BVDestroy(&Y);
391:   return 0;
392: }

394: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
395: {
396:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

398:   if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
400:   if (ctx->bs != bs) {
401:     ctx->bs = bs;
402:     eps->state = EPS_STATE_INITIAL;
403:   }
404:   return 0;
405: }

407: /*@
408:    EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.

410:    Logically Collective on eps

412:    Input Parameters:
413: +  eps - the eigenproblem solver context
414: -  bs  - the block size

416:    Options Database Key:
417: .  -eps_lobpcg_blocksize - Sets the block size

419:    Level: advanced

421: .seealso: EPSLOBPCGGetBlockSize()
422: @*/
423: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
424: {
427:   PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
428:   return 0;
429: }

431: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
432: {
433:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

435:   *bs = ctx->bs;
436:   return 0;
437: }

439: /*@
440:    EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.

442:    Not Collective

444:    Input Parameter:
445: .  eps - the eigenproblem solver context

447:    Output Parameter:
448: .  bs - the block size

450:    Level: advanced

452: .seealso: EPSLOBPCGSetBlockSize()
453: @*/
454: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
455: {
458:   PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
459:   return 0;
460: }

462: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
463: {
464:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

466:   if (restart==PETSC_DEFAULT) restart = 0.9;
468:   if (restart != ctx->restart) {
469:     ctx->restart = restart;
470:     eps->state = EPS_STATE_INITIAL;
471:   }
472:   return 0;
473: }

475: /*@
476:    EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.

478:    Logically Collective on eps

480:    Input Parameters:
481: +  eps - the eigenproblem solver context
482: -  restart - the percentage of the block of vectors to force a restart

484:    Options Database Key:
485: .  -eps_lobpcg_restart - Sets the restart parameter

487:    Notes:
488:    The meaning of this parameter is the proportion of vectors within the
489:    current block iterate that must have converged in order to force a
490:    restart with hard locking.
491:    Allowed values are in the range [0.1,1.0]. The default is 0.9.

493:    Level: advanced

495: .seealso: EPSLOBPCGGetRestart()
496: @*/
497: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
498: {
501:   PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
502:   return 0;
503: }

505: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
506: {
507:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

509:   *restart = ctx->restart;
510:   return 0;
511: }

513: /*@
514:    EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.

516:    Not Collective

518:    Input Parameter:
519: .  eps - the eigenproblem solver context

521:    Output Parameter:
522: .  restart - the restart parameter

524:    Level: advanced

526: .seealso: EPSLOBPCGSetRestart()
527: @*/
528: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
529: {
532:   PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
533:   return 0;
534: }

536: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
537: {
538:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

540:   ctx->lock = lock;
541:   return 0;
542: }

544: /*@
545:    EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
546:    the LOBPCG method.

548:    Logically Collective on eps

550:    Input Parameters:
551: +  eps  - the eigenproblem solver context
552: -  lock - true if the locking variant must be selected

554:    Options Database Key:
555: .  -eps_lobpcg_locking - Sets the locking flag

557:    Notes:
558:    This flag refers to soft locking (converged vectors within the current
559:    block iterate), since hard locking is always used (when nev is larger
560:    than the block size).

562:    Level: advanced

564: .seealso: EPSLOBPCGGetLocking()
565: @*/
566: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
567: {
570:   PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
571:   return 0;
572: }

574: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
575: {
576:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

578:   *lock = ctx->lock;
579:   return 0;
580: }

582: /*@
583:    EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.

585:    Not Collective

587:    Input Parameter:
588: .  eps - the eigenproblem solver context

590:    Output Parameter:
591: .  lock - the locking flag

593:    Level: advanced

595: .seealso: EPSLOBPCGSetLocking()
596: @*/
597: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
598: {
601:   PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
602:   return 0;
603: }

605: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
606: {
607:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
608:   PetscBool      isascii;

610:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
611:   if (isascii) {
612:     PetscViewerASCIIPrintf(viewer,"  block size %" PetscInt_FMT "\n",ctx->bs);
613:     PetscViewerASCIIPrintf(viewer,"  restart parameter=%g (using %" PetscInt_FMT " guard vectors)\n",(double)ctx->restart,ctx->guard);
614:     PetscViewerASCIIPrintf(viewer,"  soft locking %sactivated\n",ctx->lock?"":"de");
615:   }
616:   return 0;
617: }

619: PetscErrorCode EPSSetFromOptions_LOBPCG(EPS eps,PetscOptionItems *PetscOptionsObject)
620: {
621:   PetscBool      lock,flg;
622:   PetscInt       bs;
623:   PetscReal      restart;

625:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS LOBPCG Options");

627:     PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg);
628:     if (flg) EPSLOBPCGSetBlockSize(eps,bs);

630:     PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg);
631:     if (flg) EPSLOBPCGSetRestart(eps,restart);

633:     PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg);
634:     if (flg) EPSLOBPCGSetLocking(eps,lock);

636:   PetscOptionsHeadEnd();
637:   return 0;
638: }

640: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
641: {
642:   PetscFree(eps->data);
643:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
644:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
645:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
646:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
647:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
648:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
649:   return 0;
650: }

652: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
653: {
654:   EPS_LOBPCG     *lobpcg;

656:   PetscNew(&lobpcg);
657:   eps->data = (void*)lobpcg;
658:   lobpcg->lock = PETSC_TRUE;

660:   eps->useds = PETSC_TRUE;
661:   eps->categ = EPS_CATEGORY_PRECOND;

663:   eps->ops->solve          = EPSSolve_LOBPCG;
664:   eps->ops->setup          = EPSSetUp_LOBPCG;
665:   eps->ops->setupsort      = EPSSetUpSort_Default;
666:   eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
667:   eps->ops->destroy        = EPSDestroy_LOBPCG;
668:   eps->ops->view           = EPSView_LOBPCG;
669:   eps->ops->backtransform  = EPSBackTransform_Default;
670:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

672:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
673:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
674:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
675:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
676:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
677:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
678:   return 0;
679: }