Actual source code: lobpcg.c

slepc-main 2024-11-15
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: static PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
 41: {
 42:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
 43:   PetscInt   k;

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

 55: static PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
 56: {
 57:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;

 59:   PetscFunctionBegin;
 60:   EPSCheckHermitianDefinite(eps);
 61:   EPSCheckNotStructured(eps);
 62:   if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
 63:   PetscCheck(eps->n-eps->nds>=5*ctx->bs,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The problem size is too small relative to the block size");
 64:   PetscCall(EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd));
 65:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 66:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
 67:   PetscCheck(eps->which==EPS_SMALLEST_REAL || eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real or largest real eigenvalues");
 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:   PetscCall(EPSAllocateSolution(eps,0));
 78:   PetscCall(EPS_SetInnerProduct(eps));
 79:   PetscCall(DSSetType(eps->ds,DSGHEP));
 80:   PetscCall(DSAllocate(eps->ds,eps->mpd));
 81:   PetscCall(EPSSetWorkVecs(eps,1));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

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

 97:   PetscFunctionBegin;
 98:   PetscCall(STGetNumMatrices(eps->st,&nmat));
 99:   PetscCall(STGetMatrix(eps->st,0,&A));
100:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
101:   else B = NULL;

103:   if (eps->which==EPS_LARGEST_REAL) {  /* flip spectrum */
104:     flip = PETSC_TRUE;
105:     PetscCall(DSGetSlepcSC(eps->ds,&sc));
106:     sc->comparison = SlepcCompareSmallestReal;
107:   }

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

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

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

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

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

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

170:     if (ctx->lock) {
171:       PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
172:       PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
173:       if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
174:     }

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

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

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

221:       /* extend constraints */
222:       PetscCall(BVSetActiveColumns(Y,nc+locked,nc+locked+nconv));
223:       PetscCall(BVCopy(X,Y));
224:       PetscCall(BVSetActiveColumns(Y,0,nc+locked+nconv));

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

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

254:       /* B-orthogonalize initial vectors */
255:       PetscCall(BVOrthogonalize(X,NULL));
256:       PetscCall(BVSetActiveColumns(Z,nconv,ctx->bs));
257:       PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
258:       PetscCall(BVCopy(X,Z));

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

280:       continue;   /* skip the rest of the iteration */
281:     }

283:     ini = (ctx->lock)? nconv: 0;
284:     if (ctx->lock) {
285:       PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
286:       PetscCall(BVSetActiveColumns(P,nconv,ctx->bs));
287:       PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
288:       if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
289:     }

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

329:     /* 11. B-orthonormalize preconditioned residuals */
330:     PetscCall(BVOrthogonalize(R,NULL));

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

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

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

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

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

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

381: diverged:
382:   eps->its += its;

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

398: static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
399: {
400:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

402:   PetscFunctionBegin;
403:   if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 0;
404:   else PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %" PetscInt_FMT,bs);
405:   if (ctx->bs != bs) {
406:     ctx->bs = bs;
407:     eps->state = EPS_STATE_INITIAL;
408:   }
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: /*@
413:    EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.

415:    Logically Collective

417:    Input Parameters:
418: +  eps - the eigenproblem solver context
419: -  bs  - the block size

421:    Options Database Key:
422: .  -eps_lobpcg_blocksize - Sets the block size

424:    Level: advanced

426: .seealso: EPSLOBPCGGetBlockSize()
427: @*/
428: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
429: {
430:   PetscFunctionBegin;
433:   PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
438: {
439:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

441:   PetscFunctionBegin;
442:   *bs = ctx->bs;
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*@
447:    EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.

449:    Not Collective

451:    Input Parameter:
452: .  eps - the eigenproblem solver context

454:    Output Parameter:
455: .  bs - the block size

457:    Level: advanced

459: .seealso: EPSLOBPCGSetBlockSize()
460: @*/
461: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
462: {
463:   PetscFunctionBegin;
465:   PetscAssertPointer(bs,2);
466:   PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
471: {
472:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

474:   PetscFunctionBegin;
475:   if (restart==(PetscReal)PETSC_DEFAULT || restart==(PetscReal)PETSC_DECIDE) restart = 0.9;
476:   PetscCheck(restart>=0.1 && restart<=1.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The restart argument %g must be in the range [0.1,1.0]",(double)restart);
477:   if (restart != ctx->restart) {
478:     ctx->restart = restart;
479:     eps->state = EPS_STATE_INITIAL;
480:   }
481:   PetscFunctionReturn(PETSC_SUCCESS);
482: }

484: /*@
485:    EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.

487:    Logically Collective

489:    Input Parameters:
490: +  eps - the eigenproblem solver context
491: -  restart - the percentage of the block of vectors to force a restart

493:    Options Database Key:
494: .  -eps_lobpcg_restart - Sets the restart parameter

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

502:    Level: advanced

504: .seealso: EPSLOBPCGGetRestart()
505: @*/
506: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
507: {
508:   PetscFunctionBegin;
511:   PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
516: {
517:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

519:   PetscFunctionBegin;
520:   *restart = ctx->restart;
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: /*@
525:    EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.

527:    Not Collective

529:    Input Parameter:
530: .  eps - the eigenproblem solver context

532:    Output Parameter:
533: .  restart - the restart parameter

535:    Level: advanced

537: .seealso: EPSLOBPCGSetRestart()
538: @*/
539: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
540: {
541:   PetscFunctionBegin;
543:   PetscAssertPointer(restart,2);
544:   PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
549: {
550:   EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;

552:   PetscFunctionBegin;
553:   ctx->lock = lock;
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*@
558:    EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
559:    the LOBPCG method.

561:    Logically Collective

563:    Input Parameters:
564: +  eps  - the eigenproblem solver context
565: -  lock - true if the locking variant must be selected

567:    Options Database Key:
568: .  -eps_lobpcg_locking - Sets the locking flag

570:    Notes:
571:    This flag refers to soft locking (converged vectors within the current
572:    block iterate), since hard locking is always used (when nev is larger
573:    than the block size).

575:    Level: advanced

577: .seealso: EPSLOBPCGGetLocking()
578: @*/
579: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
580: {
581:   PetscFunctionBegin;
584:   PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

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

592:   PetscFunctionBegin;
593:   *lock = ctx->lock;
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

597: /*@
598:    EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.

600:    Not Collective

602:    Input Parameter:
603: .  eps - the eigenproblem solver context

605:    Output Parameter:
606: .  lock - the locking flag

608:    Level: advanced

610: .seealso: EPSLOBPCGSetLocking()
611: @*/
612: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
613: {
614:   PetscFunctionBegin;
616:   PetscAssertPointer(lock,2);
617:   PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
618:   PetscFunctionReturn(PETSC_SUCCESS);
619: }

621: static PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
622: {
623:   EPS_LOBPCG     *ctx = (EPS_LOBPCG*)eps->data;
624:   PetscBool      isascii;

626:   PetscFunctionBegin;
627:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
628:   if (isascii) {
629:     PetscCall(PetscViewerASCIIPrintf(viewer,"  block size %" PetscInt_FMT "\n",ctx->bs));
630:     PetscCall(PetscViewerASCIIPrintf(viewer,"  restart parameter=%g (using %" PetscInt_FMT " guard vectors)\n",(double)ctx->restart,ctx->guard));
631:     PetscCall(PetscViewerASCIIPrintf(viewer,"  soft locking %sactivated\n",ctx->lock?"":"de"));
632:   }
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: static PetscErrorCode EPSSetFromOptions_LOBPCG(EPS eps,PetscOptionItems *PetscOptionsObject)
637: {
638:   PetscBool      lock,flg;
639:   PetscInt       bs;
640:   PetscReal      restart;

642:   PetscFunctionBegin;
643:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS LOBPCG Options");

645:     PetscCall(PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg));
646:     if (flg) PetscCall(EPSLOBPCGSetBlockSize(eps,bs));

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

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

654:   PetscOptionsHeadEnd();
655:   PetscFunctionReturn(PETSC_SUCCESS);
656: }

658: static PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
659: {
660:   PetscFunctionBegin;
661:   PetscCall(PetscFree(eps->data));
662:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL));
663:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL));
664:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL));
665:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL));
666:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL));
667:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL));
668:   PetscFunctionReturn(PETSC_SUCCESS);
669: }

671: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
672: {
673:   EPS_LOBPCG     *lobpcg;

675:   PetscFunctionBegin;
676:   PetscCall(PetscNew(&lobpcg));
677:   eps->data = (void*)lobpcg;
678:   lobpcg->lock = PETSC_TRUE;

680:   eps->useds = PETSC_TRUE;
681:   eps->categ = EPS_CATEGORY_PRECOND;

683:   eps->ops->solve          = EPSSolve_LOBPCG;
684:   eps->ops->setup          = EPSSetUp_LOBPCG;
685:   eps->ops->setupsort      = EPSSetUpSort_Default;
686:   eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
687:   eps->ops->destroy        = EPSDestroy_LOBPCG;
688:   eps->ops->view           = EPSView_LOBPCG;
689:   eps->ops->backtransform  = EPSBackTransform_Default;
690:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

692:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG));
693:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG));
694:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG));
695:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG));
696:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG));
697:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG));
698:   PetscFunctionReturn(PETSC_SUCCESS);
699: }