Actual source code: lobpcg.c

slepc-3.10.0 2018-09-18
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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: PETSC_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: }