Actual source code: lobpcg.c

slepc-3.15.0 2021-03-31
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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;

 46:   k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
 47:   if (*ncv!=PETSC_DEFAULT) { /* 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==PETSC_DEFAULT) *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;

 61:   EPSCheckHermitianDefinite(eps);
 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==PETSC_DEFAULT) 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),PETSC_ERR_SUP,"This solver supports only smallest real or largest real eigenvalues");
 68:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
 69:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);

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

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

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

 85: PetscErrorCode EPSSolve_LOBPCG(EPS eps)
 86: {
 88:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
 89:   PetscInt       i,j,k,ld,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;

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

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-ctx->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-ctx->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-ctx->guard) {
225:       eps->its += its-1;
226:       its = 0;
227:     } else its++;

229:     if (nconv >= ctx->bs-ctx->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:     BVGetMat(R,&V);
308:     if (prev != ctx->bs-ini) {
309:       prev = ctx->bs-ini;
310:       MatDestroy(&W);
311:       MatDuplicate(V,MAT_SHARE_NONZERO_PATTERN,&W);
312:     }
313:     STApplyMat(eps->st,V,W);
314:     if (checkprecond) {
315:       for (j=ini;j<ctx->bs;j++) {
316:         MatDenseGetColumnVecRead(V,j-ini,&v);
317:         MatDenseGetColumnVecRead(W,j-ini,&w);
318:         VecDot(v,w,&dot);
319:         MatDenseRestoreColumnVecRead(W,j-ini,&w);
320:         MatDenseRestoreColumnVecRead(V,j-ini,&v);
321:         if (PetscRealPart(dot)<0.0) {
322:           PetscInfo(eps,"The preconditioner is not positive-definite\n");
323:           eps->reason = EPS_DIVERGED_BREAKDOWN;
324:           goto diverged;
325:         }
326:       }
327:     }
328:     if (nc+locked>0) {
329:       for (j=ini;j<ctx->bs;j++) {
330:         MatDenseGetColumnVecWrite(W,j-ini,&w);
331:         BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown);
332:         if (norm>0.0 && !breakdown) {
333:           VecScale(w,1.0/norm);
334:         }
335:         MatDenseRestoreColumnVecWrite(W,j-ini,&w);
336:         if (norm<=0.0 || breakdown) {
337:           PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n");
338:           eps->reason = EPS_DIVERGED_BREAKDOWN;
339:           goto diverged;
340:         }
341:       }
342:     }
343:     MatCopy(W,V,SAME_NONZERO_PATTERN);
344:     BVRestoreMat(R,&V);

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

349:     /* 13-16. B-orthonormalize conjugate directions */
350:     if (its>1) {
351:       BVOrthogonalize(P,NULL);
352:     }

354:     /* 17-23. Compute symmetric Gram matrices */
355:     BVSetActiveColumns(Z,0,ctx->bs);
356:     BVSetActiveColumns(X,0,ctx->bs);
357:     BVCopy(X,Z);
358:     BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini);
359:     BVCopy(R,Z);
360:     if (its>1) {
361:       BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini);
362:       BVCopy(P,Z);
363:     }

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

368:     BVSetActiveColumns(Z,0,nv);
369:     DSSetDimensions(eps->ds,nv,0,0,0);
370:     DSGetMat(eps->ds,DS_MAT_A,&M);
371:     BVMatProject(Z,A,Z,M);
372:     if (flip) { MatScale(M,-1.0); }
373:     DSRestoreMat(eps->ds,DS_MAT_A,&M);
374:     DSGetMat(eps->ds,DS_MAT_B,&M);
375:     BVMatProject(Z,B,Z,M); /* covers also the case B=NULL */
376:     DSRestoreMat(eps->ds,DS_MAT_B,&M);

378:     /* 24. Solve the generalized eigenvalue problem */
379:     DSSetState(eps->ds,DS_STATE_RAW);
380:     DSSolve(eps->ds,eigr,NULL);
381:     DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL);
382:     DSSynchronize(eps->ds,eigr,NULL);
383:     for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
384:     DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

386:     /* 25-33. Compute Ritz vectors */
387:     DSGetMat(eps->ds,DS_MAT_X,&M);
388:     BVSetActiveColumns(Z,ctx->bs,nv);
389:     if (ctx->lock) {
390:       BVSetActiveColumns(P,0,ctx->bs);
391:     }
392:     BVMult(P,1.0,0.0,Z,M);
393:     BVCopy(P,X);
394:     if (ctx->lock) {
395:       BVSetActiveColumns(P,nconv,ctx->bs);
396:     }
397:     BVSetActiveColumns(Z,0,ctx->bs);
398:     BVMult(X,1.0,1.0,Z,M);
399:     if (ctx->lock) {
400:       BVSetActiveColumns(X,nconv,ctx->bs);
401:     }
402:     BVMatMult(X,A,AX);
403:     MatDestroy(&M);
404:   }

406: diverged:
407:   eps->its += its;

409:   if (flip) sc->comparison = SlepcCompareLargestReal;
410:   PetscFree(eigr);
411:   MatDestroy(&W);
412:   if (V) { /* only needed when goto diverged is reached */
413:     BVRestoreMat(R,&V);
414:   }
415:   BVDestroy(&Z);
416:   BVDestroy(&X);
417:   BVDestroy(&R);
418:   BVDestroy(&P);
419:   BVDestroy(&AX);
420:   if (B) {
421:     BVDestroy(&BX);
422:   }
423:   if (nc>0 || eps->nev>ctx->bs-ctx->guard) {
424:     BVDestroy(&Y);
425:   }
426:   return(0);
427: }

429: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
430: {
431:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

434:   if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 1;
435:   if (bs <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %D",bs);
436:   if (ctx->bs != bs) {
437:     ctx->bs = bs;
438:     eps->state = EPS_STATE_INITIAL;
439:   }
440:   return(0);
441: }

443: /*@
444:    EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.

446:    Logically Collective on eps

448:    Input Parameters:
449: +  eps - the eigenproblem solver context
450: -  bs  - the block size

452:    Options Database Key:
453: .  -eps_lobpcg_blocksize - Sets the block size

455:    Level: advanced

457: .seealso: EPSLOBPCGGetBlockSize()
458: @*/
459: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
460: {

466:   PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
467:   return(0);
468: }

470: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
471: {
472:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

475:   *bs = ctx->bs;
476:   return(0);
477: }

479: /*@
480:    EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.

482:    Not Collective

484:    Input Parameter:
485: .  eps - the eigenproblem solver context

487:    Output Parameter:
488: .  bs - the block size

490:    Level: advanced

492: .seealso: EPSLOBPCGSetBlockSize()
493: @*/
494: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
495: {

501:   PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
502:   return(0);
503: }

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

510:   if (restart==PETSC_DEFAULT) restart = 0.9;
511:   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);
512:   if (restart != ctx->restart) {
513:     ctx->restart = restart;
514:     eps->state = EPS_STATE_INITIAL;
515:   }
516:   return(0);
517: }

519: /*@
520:    EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
521:    The meaning of this parameter is the proportion of vectors within the
522:    current block iterate that must have converged in order to force a
523:    restart with hard locking.

525:    Logically Collective on eps

527:    Input Parameters:
528: +  eps - the eigenproblem solver context
529: -  restart - the percentage of the block of vectors to force a restart

531:    Options Database Key:
532: .  -eps_lobpcg_restart - Sets the restart parameter

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

537:    Level: advanced

539: .seealso: EPSLOBPCGGetRestart()
540: @*/
541: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
542: {

548:   PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
549:   return(0);
550: }

552: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
553: {
554:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

557:   *restart = ctx->restart;
558:   return(0);
559: }

561: /*@
562:    EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.

564:    Not Collective

566:    Input Parameter:
567: .  eps - the eigenproblem solver context

569:    Output Parameter:
570: .  restart - the restart parameter

572:    Level: advanced

574: .seealso: EPSLOBPCGSetRestart()
575: @*/
576: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
577: {

583:   PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
584:   return(0);
585: }

587: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
588: {
589:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

592:   ctx->lock = lock;
593:   return(0);
594: }

596: /*@
597:    EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
598:    the LOBPCG method.

600:    Logically Collective on eps

602:    Input Parameters:
603: +  eps  - the eigenproblem solver context
604: -  lock - true if the locking variant must be selected

606:    Options Database Key:
607: .  -eps_lobpcg_locking - Sets the locking flag

609:    Notes:
610:    This flag refers to soft locking (converged vectors within the current
611:    block iterate), since hard locking is always used (when nev is larger
612:    than the block size).

614:    Level: advanced

616: .seealso: EPSLOBPCGGetLocking()
617: @*/
618: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
619: {

625:   PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
626:   return(0);
627: }

629: static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
630: {
631:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

634:   *lock = ctx->lock;
635:   return(0);
636: }

638: /*@
639:    EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.

641:    Not Collective

643:    Input Parameter:
644: .  eps - the eigenproblem solver context

646:    Output Parameter:
647: .  lock - the locking flag

649:    Level: advanced

651: .seealso: EPSLOBPCGSetLocking()
652: @*/
653: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
654: {

660:   PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
661:   return(0);
662: }

664: PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
665: {
667:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
668:   PetscBool      isascii;

671:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
672:   if (isascii) {
673:     PetscViewerASCIIPrintf(viewer,"  block size %D\n",ctx->bs);
674:     PetscViewerASCIIPrintf(viewer,"  restart parameter=%g (using %D guard vectors)\n",(double)ctx->restart,ctx->guard);
675:     PetscViewerASCIIPrintf(viewer,"  soft locking %sactivated\n",ctx->lock?"":"de");
676:   }
677:   return(0);
678: }

680: PetscErrorCode EPSSetFromOptions_LOBPCG(PetscOptionItems *PetscOptionsObject,EPS eps)
681: {
683:   PetscBool      lock,flg;
684:   PetscInt       bs;
685:   PetscReal      restart;

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

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

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

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

699:   PetscOptionsTail();
700:   return(0);
701: }

703: PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
704: {

708:   PetscFree(eps->data);
709:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL);
710:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL);
711:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL);
712:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL);
713:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL);
714:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL);
715:   return(0);
716: }

718: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
719: {
720:   EPS_LOBPCG     *lobpcg;

724:   PetscNewLog(eps,&lobpcg);
725:   eps->data = (void*)lobpcg;
726:   lobpcg->lock = PETSC_TRUE;

728:   eps->useds = PETSC_TRUE;
729:   eps->categ = EPS_CATEGORY_PRECOND;

731:   eps->ops->solve          = EPSSolve_LOBPCG;
732:   eps->ops->setup          = EPSSetUp_LOBPCG;
733:   eps->ops->setupsort      = EPSSetUpSort_Default;
734:   eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
735:   eps->ops->destroy        = EPSDestroy_LOBPCG;
736:   eps->ops->view           = EPSView_LOBPCG;
737:   eps->ops->backtransform  = EPSBackTransform_Default;
738:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

740:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG);
741:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG);
742:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG);
743:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG);
744:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG);
745:   PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG);
746:   return(0);
747: }