Actual source code: lobpcg.c

      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:   if (*nev==0) *nev = 1;
     47:   k = PetscMax(3*ctx->bs,((*nev-1)/ctx->bs+3)*ctx->bs);
     48:   if (*ncv!=PETSC_DETERMINE) { /* ncv set */
     49:     PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
     50:   } else *ncv = k;
     51:   if (*mpd==PETSC_DETERMINE) *mpd = 3*ctx->bs;
     52:   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");
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

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

     60:   PetscFunctionBegin;
     61:   EPSCheckHermitianDefinite(eps);
     62:   EPSCheckNotStructured(eps);
     63:   if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev?eps->nev:1);
     64:   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");
     65:   PetscCall(EPSSetDimensions_LOBPCG(eps,&eps->nev,&eps->ncv,&eps->mpd));
     66:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
     67:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
     68:   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");
     69:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD);
     70:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    416:    Logically Collective

    418:    Input Parameters:
    419: +  eps - the linear eigensolver context
    420: -  bs  - the block size

    422:    Options Database Key:
    423: .  -eps_lobpcg_blocksize bs - sets the block size

    425:    Level: advanced

    427: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGGetBlockSize()`
    428: @*/
    429: PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
    430: {
    431:   PetscFunctionBegin;
    434:   PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
    435:   PetscFunctionReturn(PETSC_SUCCESS);
    436: }

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

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

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

    450:    Not Collective

    452:    Input Parameter:
    453: .  eps - the linear eigensolver context

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

    458:    Level: advanced

    460: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGSetBlockSize()`
    461: @*/
    462: PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
    463: {
    464:   PetscFunctionBegin;
    466:   PetscAssertPointer(bs,2);
    467:   PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
    468:   PetscFunctionReturn(PETSC_SUCCESS);
    469: }

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

    475:   PetscFunctionBegin;
    476:   if (restart==(PetscReal)PETSC_DEFAULT || restart==(PetscReal)PETSC_DECIDE) restart = 0.9;
    477:   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);
    478:   if (restart != ctx->restart) {
    479:     ctx->restart = restart;
    480:     eps->state = EPS_STATE_INITIAL;
    481:   }
    482:   PetscFunctionReturn(PETSC_SUCCESS);
    483: }

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

    488:    Logically Collective

    490:    Input Parameters:
    491: +  eps - the linear eigensolver context
    492: -  restart - the percentage of the block of vectors to force a restart

    494:    Options Database Key:
    495: .  -eps_lobpcg_restart restart - sets the restart parameter

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

    503:    Level: advanced

    505: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGGetRestart()`
    506: @*/
    507: PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
    508: {
    509:   PetscFunctionBegin;
    512:   PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
    513:   PetscFunctionReturn(PETSC_SUCCESS);
    514: }

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

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

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

    528:    Not Collective

    530:    Input Parameter:
    531: .  eps - the linear eigensolver context

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

    536:    Level: advanced

    538: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGSetRestart()`
    539: @*/
    540: PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
    541: {
    542:   PetscFunctionBegin;
    544:   PetscAssertPointer(restart,2);
    545:   PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
    546:   PetscFunctionReturn(PETSC_SUCCESS);
    547: }

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

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

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

    562:    Logically Collective

    564:    Input Parameters:
    565: +  eps  - the linear eigensolver context
    566: -  lock - true if the locking variant must be selected

    568:    Options Database Key:
    569: .  -eps_lobpcg_locking (true|false) - sets the locking flag

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

    576:    Level: advanced

    578: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGGetLocking()`
    579: @*/
    580: PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
    581: {
    582:   PetscFunctionBegin;
    585:   PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
    586:   PetscFunctionReturn(PETSC_SUCCESS);
    587: }

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

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

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

    601:    Not Collective

    603:    Input Parameter:
    604: .  eps - the linear eigensolver context

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

    609:    Level: advanced

    611: .seealso: [](ch:eps), `EPSLOBPCG`, `EPSLOBPCGSetLocking()`
    612: @*/
    613: PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
    614: {
    615:   PetscFunctionBegin;
    617:   PetscAssertPointer(lock,2);
    618:   PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
    619:   PetscFunctionReturn(PETSC_SUCCESS);
    620: }

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

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

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

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

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

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

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

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

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

    672: /*MC
    673:    EPSLOBPCG - EPSLOBPCG = "lobpcg" - The Locally Optimal Block Preconditioned
    674:    Conjugate Gradient method.

    676:    Notes:
    677:    This is a preconditioned eigensolver, that is, a preconditioner can
    678:    be used to accelerate convergence, but as opposed to Davidson methods
    679:    it is intended for computing leftmost eigenvalues of Hermitian
    680:    problems, and not interior eigenvalues.

    682:    The implemented method is LOBPCG {cite:p}`Kny01`. Note that SLEPc
    683:    provides another implementation of LOBPCG via BLOPEX, see `EPSBLOPEX`.

    685:    The preconditioner is specified via the internal `ST` object and its
    686:    associated `KSP`.

    688:    Level: beginner

    690: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSGetST()`, `EPSBLOPEX`
    691: M*/
    692: SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
    693: {
    694:   EPS_LOBPCG     *lobpcg;

    696:   PetscFunctionBegin;
    697:   PetscCall(PetscNew(&lobpcg));
    698:   eps->data = (void*)lobpcg;
    699:   lobpcg->lock = PETSC_TRUE;

    701:   eps->useds = PETSC_TRUE;
    702:   eps->categ = EPS_CATEGORY_PRECOND;

    704:   eps->ops->solve          = EPSSolve_LOBPCG;
    705:   eps->ops->setup          = EPSSetUp_LOBPCG;
    706:   eps->ops->setupsort      = EPSSetUpSort_Default;
    707:   eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
    708:   eps->ops->destroy        = EPSDestroy_LOBPCG;
    709:   eps->ops->view           = EPSView_LOBPCG;
    710:   eps->ops->backtransform  = EPSBackTransform_Default;
    711:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

    713:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG));
    714:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG));
    715:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG));
    716:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG));
    717:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG));
    718:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG));
    719:   PetscFunctionReturn(PETSC_SUCCESS);
    720: }