Actual source code: lobpcg.c

slepc-3.12.1 2019-11-08
Report Typos and Errors
  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:   PetscInt  guard;     /* number of guard vectors */
 38: } EPS_LOBPCG;

 40: PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
 41: {
 42:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
 43:   PetscInt   k;

 46:   k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
 47:   if (*ncv) { /* ncv set */
 48:     if (*ncv<k) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv is not sufficiently large");
 49:   } else *ncv = k;
 50:   if (!*mpd) *mpd = 3*ctx->bs;
 51:   else if (*mpd!=3*ctx->bs) SETERRQ(PetscObjectComm((PetscObject)eps),1,"This solver does not allow a value of mpd different from 3*blocksize");
 52:   return(0);
 53: }

 55: PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
 56: {
 58:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
 59:   PetscBool      istrivial;

 62:   if (!eps->ishermitian || (eps->isgeneralized && !eps->ispositive)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"LOBPCG only works for Hermitian problems");
 63:   if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
 64:   if (eps->n-eps->nds<5*ctx->bs) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The problem size is too small relative to the block size");
 65:   EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd);
 66:   if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 67:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
 68:   if (eps->which!=EPS_SMALLEST_REAL && eps->which!=EPS_LARGEST_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
 69:   if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
 70:   if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
 71:   RGIsTrivial(eps->rg,&istrivial);
 72:   if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");

 74:   if (!ctx->restart) ctx->restart = 0.9;

 76:   /* number of guard vectors */
 77:   if (ctx->bs==1) ctx->guard = 0;
 78:   else ctx->guard = PetscMin((PetscInt)((1.0-ctx->restart)*ctx->bs+0.45),ctx->bs-1);

 80:   EPSAllocateSolution(eps,0);
 81:   EPS_SetInnerProduct(eps);
 82:   DSSetType(eps->ds,DSGHEP);
 83:   DSAllocate(eps->ds,eps->mpd);
 84:   EPSSetWorkVecs(eps,1);
 85:   return(0);
 86: }

 88: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
 89: {
 91:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
 92:   PetscInt       i,j,k,ld,nv,ini,nmat,nc,nconv,locked,its;
 93:   PetscReal      norm;
 94:   PetscScalar    *eigr,dot;
 95:   PetscBool      breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
 96:   Mat            A,B,M;
 97:   Vec            v,z,w=eps->work[0];
 98:   BV             X,Y=NULL,Z,R,P,AX,BX;
 99:   SlepcSC        sc;

102:   DSGetLeadingDimension(eps->ds,&ld);
103:   STGetNumMatrices(eps->st,&nmat);
104:   STGetMatrix(eps->st,0,&A);
105:   if (nmat>1) { STGetMatrix(eps->st,1,&B); }
106:   else B = NULL;

108:   if (eps->which==EPS_LARGEST_REAL) {  /* flip spectrum */
109:     flip = PETSC_TRUE;
110:     DSGetSlepcSC(eps->ds,&sc);
111:     sc->comparison = SlepcCompareSmallestReal;
112:   }

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

117:   /* 1. Allocate memory */
118:   PetscCalloc1(3*ctx->bs,&eigr);
119:   BVDuplicateResize(eps->V,3*ctx->bs,&Z);
120:   BVDuplicateResize(eps->V,ctx->bs,&X);
121:   BVDuplicateResize(eps->V,ctx->bs,&R);
122:   BVDuplicateResize(eps->V,ctx->bs,&P);
123:   BVDuplicateResize(eps->V,ctx->bs,&AX);
124:   if (B) {
125:     BVDuplicateResize(eps->V,ctx->bs,&BX);
126:   }
127:   nc = eps->nds;
128:   if (nc>0 || eps->nev>ctx->bs-ctx->guard) {
129:     BVDuplicateResize(eps->V,nc+eps->nev,&Y);
130:   }
131:   if (nc>0) {
132:     for (j=0;j<nc;j++) {
133:       BVGetColumn(eps->V,-nc+j,&v);
134:       BVInsertVec(Y,j,v);
135:       BVRestoreColumn(eps->V,-nc+j,&v);
136:     }
137:     BVSetActiveColumns(Y,0,nc);
138:   }

140:   /* 2. Apply the constraints to the initial vectors */
141:   /* 3. B-orthogonalize initial vectors */
142:   for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
143:     BVSetRandomColumn(eps->V,k);
144:     BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
145:   }
146:   nv = ctx->bs;
147:   BVSetActiveColumns(eps->V,0,nv);
148:   BVSetActiveColumns(Z,0,nv);
149:   BVCopy(eps->V,Z);
150:   BVCopy(Z,X);

152:   /* 4. Compute initial Ritz vectors */
153:   BVMatMult(X,A,AX);
154:   DSSetDimensions(eps->ds,nv,0,0,0);
155:   DSGetMat(eps->ds,DS_MAT_A,&M);
156:   BVMatProject(AX,NULL,X,M);
157:   if (flip) { MatScale(M,-1.0); }
158:   DSRestoreMat(eps->ds,DS_MAT_A,&M);
159:   DSSetIdentity(eps->ds,DS_MAT_B);
160:   DSSetState(eps->ds,DS_STATE_RAW);
161:   DSSolve(eps->ds,eigr,NULL);
162:   DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
163:   DSSynchronize(eps->ds,eigr,NULL);
164:   for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
165:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
166:   DSGetMat(eps->ds,DS_MAT_X,&M);
167:   BVMultInPlace(X,M,0,nv);
168:   BVMultInPlace(AX,M,0,nv);
169:   MatDestroy(&M);

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

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

179:     if (ctx->lock) {
180:       BVSetActiveColumns(R,nconv,ctx->bs);
181:       BVSetActiveColumns(AX,nconv,ctx->bs);
182:       if (B) {
183:         BVSetActiveColumns(BX,nconv,ctx->bs);
184:       }
185:     }

187:     /* 7. Compute residuals */
188:     ini = (ctx->lock)? nconv: 0;
189:     BVCopy(AX,R);
190:     if (B) { BVMatMult(X,B,BX); }
191:     for (j=ini;j<ctx->bs;j++) {
192:       BVGetColumn(R,j,&v);
193:       BVGetColumn(B?BX:X,j,&z);
194:       VecAXPY(v,-eps->eigr[locked+j],z);
195:       BVRestoreColumn(R,j,&v);
196:       BVRestoreColumn(B?BX:X,j,&z);
197:     }

199:     /* 8. Compute residual norms and update index set of active iterates */
200:     k = ini;
201:     countc = PETSC_TRUE;
202:     for (j=ini;j<ctx->bs;j++) {
203:       i = locked+j;
204:       BVGetColumn(R,j,&v);
205:       VecNorm(v,NORM_2,&norm);
206:       BVRestoreColumn(R,j,&v);
207:       (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx);
208:       if (countc) {
209:         if (eps->errest[i] < eps->tol) k++;
210:         else countc = PETSC_FALSE;
211:       }
212:       if (!countc && !eps->trackall) break;
213:     }
214:     nconv = k;
215:     eps->nconv = locked + nconv;
216:     if (its) {
217:       EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs);
218:     }
219:     (*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
220:     if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
221:       BVSetActiveColumns(eps->V,locked,eps->nconv);
222:       BVSetActiveColumns(X,0,nconv);
223:       BVCopy(X,eps->V);
224:     }
225:     if (eps->reason != EPS_CONVERGED_ITERATING) {
226:       break;
227:     } else if (nconv >= ctx->bs-ctx->guard) {
228:       eps->its += its-1;
229:       its = 0;
230:     } else its++;

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

234:       /* extend constraints */
235:       BVSetActiveColumns(Y,nc+locked,nc+locked+nconv);
236:       BVCopy(X,Y);
237:       BVSetActiveColumns(Y,0,nc+locked+nconv);

239:       /* shift work BV's */
240:       for (j=nconv;j<ctx->bs;j++) {
241:         BVCopyColumn(X,j,j-nconv);
242:         BVCopyColumn(R,j,j-nconv);
243:         BVCopyColumn(P,j,j-nconv);
244:         BVCopyColumn(AX,j,j-nconv);
245:         if (B) {
246:           BVCopyColumn(BX,j,j-nconv);
247:         }
248:       }

250:       /* set new initial vectors */
251:       BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv);
252:       BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs);
253:       BVCopy(eps->V,X);
254:       for (j=ctx->bs-nconv;j<ctx->bs;j++) {
255:         BVGetColumn(X,j,&v);
256:         BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown);
257:         if (norm>0.0 && !breakdown) {
258:           VecScale(v,1.0/norm);
259:         } else {
260:           PetscInfo(eps,"Orthogonalization of initial vector failed\n");
261:           eps->reason = EPS_DIVERGED_BREAKDOWN;
262:           goto diverged;
263:         }
264:         BVRestoreColumn(X,j,&v);
265:       }
266:       locked += nconv;
267:       nconv = 0;
268:       BVSetActiveColumns(X,nconv,ctx->bs);

270:       /* B-orthogonalize initial vectors */
271:       BVOrthogonalize(X,NULL);
272:       BVSetActiveColumns(Z,nconv,ctx->bs);
273:       BVSetActiveColumns(AX,nconv,ctx->bs);
274:       BVCopy(X,Z);

276:       /* compute initial Ritz vectors */
277:       nv = ctx->bs;
278:       BVMatMult(X,A,AX);
279:       DSSetDimensions(eps->ds,nv,0,0,0);
280:       DSGetMat(eps->ds,DS_MAT_A,&M);
281:       BVMatProject(AX,NULL,X,M);
282:       if (flip) { MatScale(M,-1.0); }
283:       DSRestoreMat(eps->ds,DS_MAT_A,&M);
284:       DSSetIdentity(eps->ds,DS_MAT_B);
285:       DSSetState(eps->ds,DS_STATE_RAW);
286:       DSSolve(eps->ds,eigr,NULL);
287:       DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
288:       DSSynchronize(eps->ds,eigr,NULL);
289:       for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
290:       DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
291:       DSGetMat(eps->ds,DS_MAT_X,&M);
292:       BVMultInPlace(X,M,0,nv);
293:       BVMultInPlace(AX,M,0,nv);
294:       MatDestroy(&M);

296:       continue;   /* skip the rest of the iteration */
297:     }

299:     ini = (ctx->lock)? nconv: 0;
300:     if (ctx->lock) {
301:       BVSetActiveColumns(R,nconv,ctx->bs);
302:       BVSetActiveColumns(P,nconv,ctx->bs);
303:       BVSetActiveColumns(AX,nconv,ctx->bs);
304:       if (B) {
305:         BVSetActiveColumns(BX,nconv,ctx->bs);
306:       }
307:     }

309:     /* 9. Apply preconditioner to the residuals */
310:     for (j=ini;j<ctx->bs;j++) {
311:       BVGetColumn(R,j,&v);
312:       STMatSolve(eps->st,v,w);
313:       if (checkprecond) {
314:         VecDot(v,w,&dot);
315:         if (PetscRealPart(dot)<0.0) {
316:           PetscInfo(eps,"The preconditioner is not positive-definite\n");
317:           eps->reason = EPS_DIVERGED_BREAKDOWN;
318:           goto diverged;
319:         }
320:       }
321:       if (nc+locked>0) {
322:         BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
323:         if (norm>0.0 && !breakdown) {
324:           VecScale(w,1.0/norm);
325:         } else {
326:           PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
327:           eps->reason = EPS_DIVERGED_BREAKDOWN;
328:           goto diverged;
329:         }
330:       }
331:       VecCopy(w,v);
332:       BVRestoreColumn(R,j,&v);
333:     }

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

338:     /* 13-16. B-orthonormalize conjugate directions */
339:     if (its>1) {
340:       BVOrthogonalize(P,NULL);
341:     }

343:     /* 17-23. Compute symmetric Gram matrices */
344:     BVSetActiveColumns(Z,0,ctx->bs);
345:     BVSetActiveColumns(X,0,ctx->bs);
346:     BVCopy(X,Z);
347:     BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
348:     BVCopy(R,Z);
349:     if (its>1) {
350:       BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
351:       BVCopy(P,Z);
352:     }

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

357:     BVSetActiveColumns(Z,0,nv);
358:     DSSetDimensions(eps->ds,nv,0,0,0);
359:     DSGetMat(eps->ds,DS_MAT_A,&M);
360:     BVMatProject(Z,A,Z,M);
361:     if (flip) { MatScale(M,-1.0); }
362:     DSRestoreMat(eps->ds,DS_MAT_A,&M);
363:     DSGetMat(eps->ds,DS_MAT_B,&M);
364:     BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
365:     DSRestoreMat(eps->ds,DS_MAT_B,&M);

367:     /* 24. Solve the generalized eigenvalue problem */
368:     DSSetState(eps->ds,DS_STATE_RAW);
369:     DSSolve(eps->ds,eigr,NULL);
370:     DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
371:     DSSynchronize(eps->ds,eigr,NULL);
372:     for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
373:     DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

375:     /* 25-33. Compute Ritz vectors */
376:     DSGetMat(eps->ds,DS_MAT_X,&M);
377:     BVSetActiveColumns(Z,ctx->bs,nv);
378:     if (ctx->lock) {
379:       BVSetActiveColumns(P,0,ctx->bs);
380:     }
381:     BVMult(P,1.0,0.0,Z,M);
382:     BVCopy(P,X);
383:     if (ctx->lock) {
384:       BVSetActiveColumns(P,nconv,ctx->bs);
385:     }
386:     BVSetActiveColumns(Z,0,ctx->bs);
387:     BVMult(X,1.0,1.0,Z,M);
388:     if (ctx->lock) {
389:       BVSetActiveColumns(X,nconv,ctx->bs);
390:     }
391:     BVMatMult(X,A,AX);
392:     MatDestroy(&M);
393:   }

395: diverged:
396:   eps->its += its;

398:   if (flip) sc->comparison = SlepcCompareLargestReal;
399:   PetscFree(eigr);
400:   BVDestroy(&Z);
401:   BVDestroy(&X);
402:   BVDestroy(&R);
403:   BVDestroy(&P);
404:   BVDestroy(&AX);
405:   if (B) {
406:     BVDestroy(&BX);
407:   }
408:   if (nc>0 || eps->nev>ctx->bs-ctx->guard) {
409:     BVDestroy(&Y);
410:   }
411:   return(0);
412: }

414: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
415: {
416:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

419:   if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
420:   if (bs <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %D",bs);
421:   if (ctx->bs != bs) {
422:     ctx->bs = bs;
423:     eps->state = EPS_STATE_INITIAL;
424:   }
425:   return(0);
426: }

428: /*@
429:    EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.

431:    Logically Collective on eps

433:    Input Parameters:
434: +  eps - the eigenproblem solver context
435: -  bs  - the block size

437:    Options Database Key:
438: .  -eps_lobpcg_blocksize - Sets the block size

440:    Level: advanced

442: .seealso: EPSLOBPCGGetBlockSize()
443: @*/
444: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
445: {

451:   PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
452:   return(0);
453: }

455: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
456: {
457:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

460:   *bs = ctx->bs;
461:   return(0);
462: }

464: /*@
465:    EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.

467:    Not Collective

469:    Input Parameter:
470: .  eps - the eigenproblem solver context

472:    Output Parameter:
473: .  bs - the block size

475:    Level: advanced

477: .seealso: EPSLOBPCGSetBlockSize()
478: @*/
479: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
480: {

486:   PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
487:   return(0);
488: }

490: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
491: {
492:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

495:   if (restart==PETSC_DEFAULT) restart = 0.9;
496:   if (restart<0.1 || restart>1.0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The restart argument %g must be in the range [0.1,1.0]",(double)restart);
497:   if (restart != ctx->restart) {
498:     ctx->restart = restart;
499:     eps->state = EPS_STATE_INITIAL;
500:   }
501:   return(0);
502: }

504: /*@
505:    EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
506:    The meaning of this parameter is the proportion of vectors within the
507:    current block iterate that must have converged in order to force a
508:    restart with hard locking.

510:    Logically Collective on eps

512:    Input Parameters:
513: +  eps - the eigenproblem solver context
514: -  restart - the percentage of the block of vectors to force a restart

516:    Options Database Key:
517: .  -eps_lobpcg_restart - Sets the restart parameter

519:    Notes:
520:    Allowed values are in the range [0.1,1.0]. The default is 0.6.

522:    Level: advanced

524: .seealso: EPSLOBPCGGetRestart()
525: @*/
526: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
527: {

533:   PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
534:   return(0);
535: }

537: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
538: {
539:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

542:   *restart = ctx->restart;
543:   return(0);
544: }

546: /*@
547:    EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.

549:    Not Collective

551:    Input Parameter:
552: .  eps - the eigenproblem solver context

554:    Output Parameter:
555: .  restart - the restart parameter

557:    Level: advanced

559: .seealso: EPSLOBPCGSetRestart()
560: @*/
561: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
562: {

568:   PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
569:   return(0);
570: }

572: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
573: {
574:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

577:   ctx->lock = lock;
578:   return(0);
579: }

581: /*@
582:    EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
583:    the LOBPCG method.

585:    Logically Collective on eps

587:    Input Parameters:
588: +  eps  - the eigenproblem solver context
589: -  lock - true if the locking variant must be selected

591:    Options Database Key:
592: .  -eps_lobpcg_locking - Sets the locking flag

594:    Notes:
595:    This flag refers to soft locking (converged vectors within the current
596:    block iterate), since hard locking is always used (when nev is larger
597:    than the block size).

599:    Level: advanced

601: .seealso: EPSLOBPCGGetLocking()
602: @*/
603: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
604: {

610:   PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
611:   return(0);
612: }

614: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
615: {
616:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

619:   *lock = ctx->lock;
620:   return(0);
621: }

623: /*@
624:    EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.

626:    Not Collective

628:    Input Parameter:
629: .  eps - the eigenproblem solver context

631:    Output Parameter:
632: .  lock - the locking flag

634:    Level: advanced

636: .seealso: EPSLOBPCGSetLocking()
637: @*/
638: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
639: {

645:   PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
646:   return(0);
647: }

649: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
650: {
652:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
653:   PetscBool      isascii;

656:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
657:   if (isascii) {
658:     PetscViewerASCIIPrintf(viewer,"  block size %D\n",ctx->bs);
659:     PetscViewerASCIIPrintf(viewer,"  restart parameter=%g (using %D guard vectors)\n",(double)ctx->restart,ctx->guard);
660:     PetscViewerASCIIPrintf(viewer,"  soft locking %sactivated\n",ctx->lock?"":"de");
661:   }
662:   return(0);
663: }

665: PetscErrorCode EPSSetFromOptions_LOBPCG(PetscOptionItems *PetscOptionsObject,EPS eps)
666: {
668:   PetscBool      lock,flg;
669:   PetscInt       bs;
670:   PetscReal      restart;

673:   PetscOptionsHead(PetscOptionsObject,"EPS LOBPCG Options");

675:     PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg);
676:     if (flg) { EPSLOBPCGSetBlockSize(eps,bs); }

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

681:     PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg);
682:     if (flg) { EPSLOBPCGSetLocking(eps,lock); }

684:   PetscOptionsTail();
685:   return(0);
686: }

688: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
689: {

693:   PetscFree(eps->data);
694:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
695:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
696:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
697:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
698:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
699:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
700:   return(0);
701: }

703: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
704: {
705:   EPS_LOBPCG     *lobpcg;

709:   PetscNewLog(eps,&lobpcg);
710:   eps->data = (void*)lobpcg;
711:   lobpcg->lock = PETSC_TRUE;

713:   eps->useds = PETSC_TRUE;
714:   eps->categ = EPS_CATEGORY_PRECOND;

716:   eps->ops->solve          = EPSSolve_LOBPCG;
717:   eps->ops->setup          = EPSSetUp_LOBPCG;
718:   eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
719:   eps->ops->destroy        = EPSDestroy_LOBPCG;
720:   eps->ops->view           = EPSView_LOBPCG;
721:   eps->ops->backtransform  = EPSBackTransform_Default;
722:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

724:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
725:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
726:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
727:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
728:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
729:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
730:   return(0);
731: }