Actual source code: power.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: "power"

     13:    Method: Power Iteration

     15:    Algorithm:

     17:        This solver implements the power iteration for finding dominant
     18:        eigenpairs. It also includes the following well-known methods:
     19:        - Inverse Iteration: when used in combination with shift-and-invert
     20:          spectral transformation.
     21:        - Rayleigh Quotient Iteration (RQI): also with shift-and-invert plus
     22:          a variable shift.

     24:        It can also be used for nonlinear inverse iteration on the problem
     25:        A(x)*x=lambda*B(x)*x, where A and B are not constant but depend on x.

     27:    References:

     29:        [1] "Single Vector Iteration Methods in SLEPc", SLEPc Technical Report
     30:            STR-2, available at https://slepc.upv.es.
     31: */

     33: #include <slepc/private/epsimpl.h>
     34: #include <slepcblaslapack.h>
     35: /* PETSc headers */
     36: #include <petscdm.h>
     37: #include <petscsnes.h>

     39: static PetscErrorCode EPSPowerFormFunction_Update(SNES,Vec,Vec,void*);
     40: static PetscErrorCode EPSSolve_Power(EPS);
     41: static PetscErrorCode EPSSolve_TS_Power(EPS);

     43: typedef struct {
     44:   EPSPowerShiftType shift_type;
     45:   PetscBool         nonlinear;
     46:   PetscBool         update;
     47:   PetscBool         sign_normalization;
     48:   SNES              snes;
     49:   PetscErrorCode    (*formFunctionB)(SNES,Vec,Vec,void*);
     50:   void              *formFunctionBctx;
     51:   PetscErrorCode    (*formFunctionA)(SNES,Vec,Vec,void*);
     52:   void              *formFunctionActx;
     53:   PetscErrorCode    (*formFunctionAB)(SNES,Vec,Vec,Vec,void*);
     54:   PetscErrorCode    (*formNorm)(SNES,Vec,PetscReal*,void*);
     55:   void              *formNormCtx;
     56:   PetscInt          idx;  /* index of the first nonzero entry in the iteration vector */
     57:   PetscMPIInt       p;    /* process id of the owner of idx */
     58:   PetscReal         norm0; /* norm of initial vector */
     59: } EPS_POWER;

     61: static PetscErrorCode SNESMonitor_PowerUpdate(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
     62: {
     63:   EPS            eps;

     65:   PetscFunctionBegin;
     66:   PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
     67:   PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
     68:   /* Call EPS monitor on each SNES iteration */
     69:   PetscCall(EPSMonitor(eps,its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: static PetscErrorCode EPSSetUp_Power(EPS eps)
     74: {
     75:   EPS_POWER      *power = (EPS_POWER*)eps->data;
     76:   STMatMode      mode;
     77:   Mat            A,B,P;
     78:   Vec            res;
     79:   PetscContainer container;
     80:   PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
     81:   PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
     82:   PetscErrorCode (*formNorm)(SNES,Vec,PetscReal*,void*);
     83:   void           *ctx;

     85:   PetscFunctionBegin;
     86:   EPSCheckNotStructured(eps);
     87:   if (eps->nev==0) eps->nev = 1;
     88:   if (eps->ncv!=PETSC_DETERMINE) {
     89:     PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
     90:   } else eps->ncv = eps->nev;
     91:   if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
     92:   if (eps->max_it==PETSC_DETERMINE) {
     93:     /* SNES will directly return the solution for us, and we need to do only one iteration */
     94:     if (power->nonlinear && power->update) eps->max_it = 1;
     95:     else eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
     96:   }
     97:   if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
     98:   PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
     99:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
    100:     PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
    101:     EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," (with variable shifts)");
    102:     PetscCall(STGetMatMode(eps->st,&mode));
    103:     PetscCheck(mode!=ST_MATMODE_INPLACE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
    104:   }
    105:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_THRESHOLD);
    106:   EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
    107:   PetscCall(EPSAllocateSolution(eps,0));
    108:   PetscCall(EPS_SetInnerProduct(eps));

    110:   if (power->nonlinear) {
    111:     PetscCheck(eps->nev==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
    112:     PetscCall(EPSSetWorkVecs(eps,3));
    113:     PetscCheck(!power->update || eps->max_it==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"More than one iteration is not allowed for Newton eigensolver (SNES)");

    115:     /* Set up SNES solver */
    116:     /* If SNES was setup, we need to reset it so that it will be consistent with the current EPS */
    117:     if (power->snes) PetscCall(SNESReset(power->snes));
    118:     else PetscCall(EPSPowerGetSNES(eps,&power->snes));

    120:     /* For nonlinear solver (Newton), we should scale the initial vector back.
    121:        The initial vector will be scaled in EPSSetUp. */
    122:     if (eps->IS) PetscCall(VecNorm(eps->IS[0],NORM_2,&power->norm0));

    124:     PetscCall(EPSGetOperators(eps,&A,&B));

    126:     /* form function */
    127:     PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA));
    128:     PetscCheck(formFunctionA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
    129:     PetscCall(PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container));
    130:     if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
    131:     else ctx = NULL;
    132:     PetscCall(MatCreateVecs(A,&res,NULL));
    133:     power->formFunctionA = formFunctionA;
    134:     power->formFunctionActx = ctx;
    135:     if (power->update) {
    136:       PetscCall(SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx));
    137:       PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB));
    138:       PetscCall(SNESMonitorSet(power->snes,SNESMonitor_PowerUpdate,NULL,NULL));
    139:     }
    140:     else PetscCall(SNESSetFunction(power->snes,res,formFunctionA,ctx));
    141:     PetscCall(VecDestroy(&res));

    143:     /* form Jacobian */
    144:     PetscCall(PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA));
    145:     PetscCheck(formJacobianA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
    146:     PetscCall(PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container));
    147:     if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
    148:     else ctx = NULL;
    149:     /* This allows users to compute a different preconditioning matrix than A.
    150:      * It is useful when A and B are shell matrices.
    151:      */
    152:     PetscCall(STGetPreconditionerMat(eps->st,&P));
    153:     /* If users do not provide a matrix, we simply use A */
    154:     PetscCall(SNESSetJacobian(power->snes,A,P? P:A,formJacobianA,ctx));
    155:     PetscCall(SNESSetFromOptions(power->snes));
    156:     PetscCall(SNESSetUp(power->snes));

    158:     ctx = NULL;
    159:     formNorm = NULL;
    160:     if (B) {
    161:       PetscCall(PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB));
    162:       PetscCall(PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container));
    163:       if (power->formFunctionB && container) PetscCall(PetscContainerGetPointer(container,&power->formFunctionBctx));
    164:       else power->formFunctionBctx = NULL;

    166:       /* form norm */
    167:       PetscCall(PetscObjectQueryFunction((PetscObject)B,"formNorm",&formNorm));
    168:       if (formNorm) {
    169:         PetscCall(PetscObjectQuery((PetscObject)B,"formNormCtx",(PetscObject*)&container));
    170:         if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
    171:       }
    172:     }
    173:     power->formNorm = formNorm;
    174:     power->formNormCtx = ctx;
    175:   } else {
    176:     if (eps->twosided) PetscCall(EPSSetWorkVecs(eps,3));
    177:     else PetscCall(EPSSetWorkVecs(eps,2));
    178:     PetscCall(DSSetType(eps->ds,DSNHEP));
    179:     PetscCall(DSAllocate(eps->ds,eps->nev));
    180:   }
    181:   /* dispatch solve method */
    182:   if (eps->twosided) {
    183:     PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
    184:     PetscCheck(power->shift_type!=EPS_POWER_SHIFT_WILKINSON,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
    185:     eps->ops->solve = EPSSolve_TS_Power;
    186:   } else eps->ops->solve = EPSSolve_Power;
    187:   PetscFunctionReturn(PETSC_SUCCESS);
    188: }

    190: /*
    191:    Find the index of the first nonzero entry of x, and the process that owns it.
    192: */
    193: static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscMPIInt *p)
    194: {
    195:   PetscInt          i,first,last,N;
    196:   PetscLayout       map;
    197:   const PetscScalar *xx;

    199:   PetscFunctionBegin;
    200:   PetscCall(VecGetSize(x,&N));
    201:   PetscCall(VecGetOwnershipRange(x,&first,&last));
    202:   PetscCall(VecGetArrayRead(x,&xx));
    203:   for (i=first;i<last;i++) {
    204:     if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
    205:   }
    206:   PetscCall(VecRestoreArrayRead(x,&xx));
    207:   if (i==last) i=N;
    208:   PetscCallMPI(MPIU_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x)));
    209:   PetscCheck(*idx!=N,PetscObjectComm((PetscObject)x),PETSC_ERR_PLIB,"Zero vector found");
    210:   PetscCall(VecGetLayout(x,&map));
    211:   PetscCall(PetscLayoutFindOwner(map,*idx,p));
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: /*
    216:    Normalize a vector x with respect to a given norm as well as, optionally, the
    217:    sign of the first nonzero entry (assumed to be idx in process p).
    218: */
    219: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscMPIInt p,PetscBool sign_normalization,PetscScalar *sign)
    220: {
    221:   PetscScalar       alpha=1.0;
    222:   PetscInt          first,last;
    223:   PetscReal         absv;
    224:   const PetscScalar *xx;

    226:   PetscFunctionBegin;
    227:   if (sign_normalization) {
    228:     PetscCall(VecGetOwnershipRange(x,&first,&last));
    229:     if (idx>=first && idx<last) {
    230:       PetscCall(VecGetArrayRead(x,&xx));
    231:       absv = PetscAbsScalar(xx[idx-first]);
    232:       if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
    233:       PetscCall(VecRestoreArrayRead(x,&xx));
    234:     }
    235:     PetscCallMPI(MPI_Bcast(&alpha,1,MPIU_SCALAR,p,PetscObjectComm((PetscObject)x)));
    236:   }
    237:   if (sign) *sign = alpha;
    238:   alpha *= norm;
    239:   PetscCall(VecScale(x,1.0/alpha));
    240:   PetscFunctionReturn(PETSC_SUCCESS);
    241: }

    243: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
    244: {
    245:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    246:   Mat            B;

    248:   PetscFunctionBegin;
    249:   PetscCall(STResetMatrixState(eps->st));
    250:   PetscCall(EPSGetOperators(eps,NULL,&B));
    251:   if (B) {
    252:     if (power->formFunctionB) PetscCall((*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx));
    253:     else PetscCall(MatMult(B,x,Bx));
    254:   } else PetscCall(VecCopy(x,Bx));
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
    259: {
    260:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    261:   Mat            A;

    263:   PetscFunctionBegin;
    264:   PetscCall(STResetMatrixState(eps->st));
    265:   PetscCall(EPSGetOperators(eps,&A,NULL));
    266:   PetscCheck(A,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
    267:   if (power->formFunctionA) PetscCall((*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx));
    268:   else PetscCall(MatMult(A,x,Ax));
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
    273: {
    274:   EPS            eps;
    275:   PetscReal      bx;
    276:   Vec            Bx;
    277:   PetscScalar    sign;
    278:   EPS_POWER      *power;

    280:   PetscFunctionBegin;
    281:   PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
    282:   PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
    283:   PetscCall(STResetMatrixState(eps->st));
    284:   power = (EPS_POWER*)eps->data;
    285:   Bx = eps->work[2];
    286:   if (power->formFunctionAB) PetscCall((*power->formFunctionAB)(snes,x,y,Bx,ctx));
    287:   else {
    288:     PetscCall(EPSPowerUpdateFunctionA(eps,x,y));
    289:     PetscCall(EPSPowerUpdateFunctionB(eps,x,Bx));
    290:   }
    291:   if (power->formNorm) PetscCall((*power->formNorm)(snes,Bx,&bx,power->formNormCtx));
    292:   else PetscCall(VecNorm(Bx,NORM_2,&bx));
    293:   PetscCall(Normalize(Bx,bx,power->idx,power->p,power->sign_normalization,&sign));
    294:   PetscCall(VecAXPY(y,-1.0,Bx));
    295:   /* Keep tracking eigenvalue update. It would be useful when we want to monitor solver progress via snes monitor. */
    296:   eps->eigr[(eps->nconv < eps->nev)? eps->nconv:(eps->nconv-1)] = 1.0/((PetscScalar)bx*sign);
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: /*
    301:    Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
    302: */
    303: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
    304: {
    305:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    306:   Vec            Bx;

    308:   PetscFunctionBegin;
    309:   PetscCall(VecCopy(x,y));
    310:   if (power->update) PetscCall(SNESSolve(power->snes,NULL,y));
    311:   else {
    312:     Bx = eps->work[2];
    313:     PetscCall(SNESSolve(power->snes,Bx,y));
    314:   }
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: /*
    319:    Use nonlinear inverse power to compute an initial guess
    320: */
    321: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
    322: {
    323:   EPS            powereps;
    324:   Mat            A,B,P;
    325:   Vec            v1,v2;
    326:   SNES           snes;
    327:   DM             dm,newdm;
    328:   PetscBool      sign_normalization;

    330:   PetscFunctionBegin;
    331:   PetscCall(EPSCreate(PetscObjectComm((PetscObject)eps),&powereps));
    332:   PetscCall(EPSGetOperators(eps,&A,&B));
    333:   PetscCall(EPSSetType(powereps,EPSPOWER));
    334:   PetscCall(EPSSetOperators(powereps,A,B));
    335:   PetscCall(EPSSetTolerances(powereps,1e-6,4));
    336:   PetscCall(EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix));
    337:   PetscCall(EPSAppendOptionsPrefix(powereps,"init_"));
    338:   PetscCall(EPSSetProblemType(powereps,EPS_GNHEP));
    339:   PetscCall(EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE));
    340:   PetscCall(EPSPowerSetNonlinear(powereps,PETSC_TRUE));
    341:   PetscCall(STGetPreconditionerMat(eps->st,&P));
    342:   /* attach dm to initial solve */
    343:   PetscCall(EPSPowerGetSNES(eps,&snes));
    344:   PetscCall(SNESGetDM(snes,&dm));
    345:   /* use  dmshell to temporarily store snes context */
    346:   PetscCall(DMCreate(PetscObjectComm((PetscObject)eps),&newdm));
    347:   PetscCall(DMSetType(newdm,DMSHELL));
    348:   PetscCall(DMSetUp(newdm));
    349:   PetscCall(DMCopyDMSNES(dm,newdm));
    350:   PetscCall(EPSPowerGetSNES(powereps,&snes));
    351:   PetscCall(SNESSetDM(snes,dm));
    352:   PetscCall(EPSPowerGetSignNormalization(eps,&sign_normalization));
    353:   PetscCall(EPSPowerSetSignNormalization(powereps,sign_normalization));
    354:   PetscCall(EPSSetFromOptions(powereps));
    355:   if (P) PetscCall(STSetPreconditionerMat(powereps->st,P));
    356:   PetscCall(EPSSolve(powereps));
    357:   PetscCall(BVGetColumn(eps->V,0,&v2));
    358:   PetscCall(BVGetColumn(powereps->V,0,&v1));
    359:   PetscCall(VecCopy(v1,v2));
    360:   PetscCall(BVRestoreColumn(powereps->V,0,&v1));
    361:   PetscCall(BVRestoreColumn(eps->V,0,&v2));
    362:   PetscCall(EPSDestroy(&powereps));
    363:   /* restore context back to the old nonlinear solver */
    364:   PetscCall(DMCopyDMSNES(newdm,dm));
    365:   PetscCall(DMDestroy(&newdm));
    366:   PetscFunctionReturn(PETSC_SUCCESS);
    367: }

    369: static PetscErrorCode EPSSolve_Power(EPS eps)
    370: {
    371:   EPS_POWER           *power = (EPS_POWER*)eps->data;
    372:   PetscInt            k,ld;
    373:   Vec                 v,y,e,Bx;
    374:   Mat                 A;
    375:   KSP                 ksp;
    376:   PetscReal           relerr,norm,norm1,rt1,rt2,cs1;
    377:   PetscScalar         theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
    378:   PetscBool           breakdown;
    379:   KSPConvergedReason  reason;
    380:   SNESConvergedReason snesreason;

    382:   PetscFunctionBegin;
    383:   e = eps->work[0];
    384:   y = eps->work[1];
    385:   if (power->nonlinear) Bx = eps->work[2];
    386:   else Bx = NULL;

    388:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
    389:   if (power->nonlinear) {
    390:     PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps));
    391:     /* Compute an initial guess only when users do not provide one */
    392:     if (power->update && !eps->nini) PetscCall(EPSPowerComputeInitialGuess_Update(eps));
    393:   } else PetscCall(DSGetLeadingDimension(eps->ds,&ld));
    394:   if (!power->update) PetscCall(EPSGetStartVector(eps,0,NULL));
    395:   if (power->nonlinear) {
    396:     PetscCall(BVGetColumn(eps->V,0,&v));
    397:     if (eps->nini) {
    398:       /* We scale the initial vector back if the initial vector was provided by users */
    399:       PetscCall(VecScale(v,power->norm0));
    400:     }
    401:     /* Do a couple things:
    402:      * 1) Determine the first non-zero index for Bx and the proc that owns it. This will be
    403:      *    used if performing normalization by the sign of the first nonzero element of Bx.
    404:      * 2) Scale Bx by its norm. For non-update power iterations, Bx (in code terms) is used
    405:      *    as the RHS argument to SNESSolve. And recall that the formula for generalized
    406:      *    inverse power iterations in this case is: (Ax)_n = (Bx)_{n-1}/|(Bx)_{n-1}| (in
    407:      *    math terms)
    408:      */
    409:     PetscCall(EPSPowerUpdateFunctionB(eps,v,Bx));
    410:     PetscCall(BVRestoreColumn(eps->V,0,&v));
    411:     if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
    412:     else PetscCall(VecNorm(Bx,NORM_2,&norm));
    413:     PetscCall(FirstNonzeroIdx(Bx,&power->idx,&power->p));
    414:     PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,NULL));
    415:   }

    417:   PetscCall(STGetShift(eps->st,&sigma));    /* original shift */
    418:   rho = sigma;

    420:   while (eps->reason == EPS_CONVERGED_ITERATING) {
    421:     eps->its++;
    422:     k = eps->nconv;

    424:     /* y = OP v */
    425:     PetscCall(BVGetColumn(eps->V,k,&v));
    426:     if (power->nonlinear) PetscCall(EPSPowerApply_SNES(eps,v,y));
    427:     else PetscCall(STApply(eps->st,v,y));
    428:     PetscCall(BVRestoreColumn(eps->V,k,&v));

    430:     /* purge previously converged eigenvectors */
    431:     if (PetscUnlikely(power->nonlinear)) {
    432:       /* We do not need to call this for Newton eigensolver since eigenvalue is
    433:        * updated in function evaluations.
    434:        */
    435:       if (!power->update) {
    436:         PetscCall(EPSPowerUpdateFunctionB(eps,y,Bx));
    437:         if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
    438:         else PetscCall(VecNorm(Bx,NORM_2,&norm));
    439:         PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,&sign));
    440:       }
    441:     } else {
    442:       PetscCall(DSGetArray(eps->ds,DS_MAT_A,&T));
    443:       PetscCall(BVSetActiveColumns(eps->V,0,k));
    444:       PetscCall(BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL));
    445:     }

    447:     /* theta = (v,y)_B */
    448:     PetscCall(BVSetActiveColumns(eps->V,k,k+1));
    449:     PetscCall(BVDotVec(eps->V,y,&theta));
    450:     if (!power->nonlinear) {
    451:       T[k+k*ld] = theta;
    452:       PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&T));
    453:     }

    455:     /* Eigenvalue is already stored in function evaluations.
    456:      * Assign eigenvalue to theta to make the rest of the code consistent
    457:      */
    458:     if (power->update) theta = eps->eigr[eps->nconv];
    459:     else if (power->nonlinear) theta = 1.0/((PetscScalar)norm*sign); /* Eigenvalue: 1/|Bx| */

    461:     if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */

    463:       /* approximate eigenvalue is the Rayleigh quotient */
    464:       eps->eigr[eps->nconv] = theta;

    466:       /**
    467:        * If the Newton method (update, SNES) is used, we do not compute "relerr"
    468:        * since SNES determines the convergence.
    469:        */
    470:       if (PetscUnlikely(power->update)) relerr = 0.;
    471:       else {
    472:         /* compute relative error as ||y-theta v||_2/|theta| */
    473:         PetscCall(VecCopy(y,e));
    474:         PetscCall(BVGetColumn(eps->V,k,&v));
    475:         PetscCall(VecAXPY(e,power->nonlinear?-1.0:-theta,v));
    476:         PetscCall(BVRestoreColumn(eps->V,k,&v));
    477:         PetscCall(VecNorm(e,NORM_2,&relerr));
    478:         if (PetscUnlikely(power->nonlinear)) relerr *= PetscAbsScalar(theta);
    479:         else relerr /= PetscAbsScalar(theta);
    480:       }

    482:     } else {  /* RQI */

    484:       /* delta = ||y||_B */
    485:       delta = norm;

    487:       /* compute relative error */
    488:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
    489:       else relerr = 1.0 / (norm*PetscAbsScalar(rho));

    491:       /* approximate eigenvalue is the shift */
    492:       eps->eigr[eps->nconv] = rho;

    494:       /* compute new shift */
    495:       if (relerr<eps->tol) {
    496:         rho = sigma;  /* if converged, restore original shift */
    497:         PetscCall(STSetShift(eps->st,rho));
    498:       } else {
    499:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
    500:         if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
    501:           /* beta1 is the norm of the residual associated with R(v) */
    502:           PetscCall(BVGetColumn(eps->V,k,&v));
    503:           PetscCall(VecAXPY(v,-PetscConj(theta)/(delta*delta),y));
    504:           PetscCall(BVRestoreColumn(eps->V,k,&v));
    505:           PetscCall(BVScaleColumn(eps->V,k,1.0/delta));
    506:           PetscCall(BVNormColumn(eps->V,k,NORM_2,&norm1));
    507:           beta1 = norm1;

    509:           /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
    510:           PetscCall(STGetMatrix(eps->st,0,&A));
    511:           PetscCall(BVGetColumn(eps->V,k,&v));
    512:           PetscCall(MatMult(A,v,e));
    513:           PetscCall(VecDot(v,e,&alpha2));
    514:           PetscCall(BVRestoreColumn(eps->V,k,&v));
    515:           alpha2 = alpha2 / (beta1 * beta1);

    517:           /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
    518:           PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    519:           PetscCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
    520:           PetscCall(PetscFPTrapPop());
    521:           if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
    522:           else rho = rt2;
    523:         }
    524:         /* update operator according to new shift */
    525:         PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
    526:         PetscCall(STSetShift(eps->st,rho));
    527:         PetscCall(KSPGetConvergedReason(ksp,&reason));
    528:         if (reason) {
    529:           PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
    530:           rho *= 1+10*PETSC_MACHINE_EPSILON;
    531:           PetscCall(STSetShift(eps->st,rho));
    532:           PetscCall(KSPGetConvergedReason(ksp,&reason));
    533:           PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
    534:         }
    535:         PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
    536:       }
    537:     }
    538:     eps->errest[eps->nconv] = relerr;

    540:     /* normalize */
    541:     if (!power->nonlinear) PetscCall(Normalize(y,norm,power->idx,power->p,power->sign_normalization,NULL));
    542:     PetscCall(BVInsertVec(eps->V,k,y));

    544:     if (PetscUnlikely(power->update)) {
    545:       PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
    546:       /* For Newton eigensolver, we are ready to return once SNES converged. */
    547:       if (snesreason>0) eps->nconv = 1;
    548:     } else if (PetscUnlikely(relerr<eps->tol)) {   /* accept eigenpair */
    549:       eps->nconv = eps->nconv + 1;
    550:       if (eps->nconv<eps->nev) {
    551:         PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
    552:         if (breakdown) {
    553:           eps->reason = EPS_DIVERGED_BREAKDOWN;
    554:           PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
    555:           break;
    556:         }
    557:       }
    558:     }
    559:     /* For Newton eigensolver, monitor will be called from SNES monitor */
    560:     if (!power->update) PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
    561:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));

    563:     /**
    564:      * When a customized stopping test is used, and reason can be set to be converged (EPS_CONVERGED_USER).
    565:      * In this case, we need to increase eps->nconv to "1" so users can retrieve the solution.
    566:      */
    567:     if (PetscUnlikely(power->nonlinear && eps->reason>0)) eps->nconv = 1;
    568:   }

    570:   if (power->nonlinear) PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",NULL));
    571:   else {
    572:     PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
    573:     PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
    574:   }
    575:   PetscFunctionReturn(PETSC_SUCCESS);
    576: }

    578: static PetscErrorCode EPSSolve_TS_Power(EPS eps)
    579: {
    580:   EPS_POWER          *power = (EPS_POWER*)eps->data;
    581:   PetscInt           k,ld;
    582:   Vec                v,w,y,e,z;
    583:   KSP                ksp;
    584:   PetscReal          relerr=1.0,relerrl,delta;
    585:   PetscScalar        theta,rho,alpha,sigma;
    586:   PetscBool          breakdown,breakdownl;
    587:   KSPConvergedReason reason;

    589:   PetscFunctionBegin;
    590:   e = eps->work[0];
    591:   v = eps->work[1];
    592:   w = eps->work[2];

    594:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
    595:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
    596:   PetscCall(EPSGetStartVector(eps,0,NULL));
    597:   PetscCall(EPSGetLeftStartVector(eps,0,NULL));
    598:   PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL));
    599:   PetscCall(BVCopyVec(eps->V,0,v));
    600:   PetscCall(BVCopyVec(eps->W,0,w));
    601:   PetscCall(STGetShift(eps->st,&sigma));    /* original shift */
    602:   rho = sigma;

    604:   while (eps->reason == EPS_CONVERGED_ITERATING) {
    605:     eps->its++;
    606:     k = eps->nconv;

    608:     /* y = OP v, z = OP' w */
    609:     PetscCall(BVGetColumn(eps->V,k,&y));
    610:     PetscCall(STApply(eps->st,v,y));
    611:     PetscCall(BVRestoreColumn(eps->V,k,&y));
    612:     PetscCall(BVGetColumn(eps->W,k,&z));
    613:     PetscCall(STApplyHermitianTranspose(eps->st,w,z));
    614:     PetscCall(BVRestoreColumn(eps->W,k,&z));

    616:     /* purge previously converged eigenvectors */
    617:     PetscCall(BVBiorthogonalizeColumn(eps->V,eps->W,k));

    619:     /* theta = (w,y)_B */
    620:     PetscCall(BVSetActiveColumns(eps->V,k,k+1));
    621:     PetscCall(BVDotVec(eps->V,w,&theta));
    622:     theta = PetscConj(theta);

    624:     if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */

    626:       /* approximate eigenvalue is the Rayleigh quotient */
    627:       eps->eigr[eps->nconv] = theta;

    629:       /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
    630:       PetscCall(BVCopyVec(eps->V,k,e));
    631:       PetscCall(VecAXPY(e,-theta,v));
    632:       PetscCall(VecNorm(e,NORM_2,&relerr));
    633:       PetscCall(BVCopyVec(eps->W,k,e));
    634:       PetscCall(VecAXPY(e,-PetscConj(theta),w));
    635:       PetscCall(VecNorm(e,NORM_2,&relerrl));
    636:       relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
    637:     }

    639:     /* normalize */
    640:     PetscCall(BVSetActiveColumns(eps->V,k,k+1));
    641:     PetscCall(BVGetColumn(eps->W,k,&z));
    642:     PetscCall(BVDotVec(eps->V,z,&alpha));
    643:     PetscCall(BVRestoreColumn(eps->W,k,&z));
    644:     delta = PetscSqrtReal(PetscAbsScalar(alpha));
    645:     PetscCheck(delta!=0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in two-sided Power/RQI");
    646:     PetscCall(BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta)));
    647:     PetscCall(BVScaleColumn(eps->W,k,1.0/delta));
    648:     PetscCall(BVCopyVec(eps->V,k,v));
    649:     PetscCall(BVCopyVec(eps->W,k,w));

    651:     if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */

    653:       /* compute relative error */
    654:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
    655:       else relerr = 1.0 / PetscAbsScalar(delta*rho);

    657:       /* approximate eigenvalue is the shift */
    658:       eps->eigr[eps->nconv] = rho;

    660:       /* compute new shift */
    661:       if (relerr<eps->tol) {
    662:         rho = sigma;  /* if converged, restore original shift */
    663:         PetscCall(STSetShift(eps->st,rho));
    664:       } else {
    665:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
    666:         /* update operator according to new shift */
    667:         PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
    668:         PetscCall(STSetShift(eps->st,rho));
    669:         PetscCall(KSPGetConvergedReason(ksp,&reason));
    670:         if (reason) {
    671:           PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
    672:           rho *= 1+10*PETSC_MACHINE_EPSILON;
    673:           PetscCall(STSetShift(eps->st,rho));
    674:           PetscCall(KSPGetConvergedReason(ksp,&reason));
    675:           PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
    676:         }
    677:         PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
    678:       }
    679:     }
    680:     eps->errest[eps->nconv] = relerr;

    682:     /* if relerr<tol, accept eigenpair */
    683:     if (relerr<eps->tol) {
    684:       eps->nconv = eps->nconv + 1;
    685:       if (eps->nconv<eps->nev) {
    686:         PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
    687:         PetscCall(EPSGetLeftStartVector(eps,eps->nconv,&breakdownl));
    688:         if (breakdown || breakdownl) {
    689:           eps->reason = EPS_DIVERGED_BREAKDOWN;
    690:           PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
    691:           break;
    692:         }
    693:         PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL));
    694:       }
    695:     }
    696:     PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
    697:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
    698:   }

    700:   PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
    701:   PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

    705: static PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
    706: {
    707:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    708:   SNESConvergedReason snesreason;

    710:   PetscFunctionBegin;
    711:   if (PetscUnlikely(power->update)) {
    712:     PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
    713:     if (snesreason < 0) {
    714:       *reason = EPS_DIVERGED_BREAKDOWN;
    715:       PetscFunctionReturn(PETSC_SUCCESS);
    716:     }
    717:   }
    718:   PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx));
    719:   PetscFunctionReturn(PETSC_SUCCESS);
    720: }

    722: static PetscErrorCode EPSBackTransform_Power(EPS eps)
    723: {
    724:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    726:   PetscFunctionBegin;
    727:   if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
    728:   else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) PetscCall(EPSBackTransform_Default(eps));
    729:   PetscFunctionReturn(PETSC_SUCCESS);
    730: }

    732: static PetscErrorCode EPSSetFromOptions_Power(EPS eps,PetscOptionItems PetscOptionsObject)
    733: {
    734:   EPS_POWER         *power = (EPS_POWER*)eps->data;
    735:   PetscBool         flg,val;
    736:   EPSPowerShiftType shift;

    738:   PetscFunctionBegin;
    739:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Power Options");

    741:     PetscCall(PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg));
    742:     if (flg) PetscCall(EPSPowerSetShiftType(eps,shift));

    744:     PetscCall(PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg));
    745:     if (flg) PetscCall(EPSPowerSetNonlinear(eps,val));

    747:     PetscCall(PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg));
    748:     if (flg) PetscCall(EPSPowerSetUpdate(eps,val));

    750:     PetscCall(PetscOptionsBool("-eps_power_sign_normalization","Normalize Bx with sign of first nonzero entry","EPSPowerSetSignNormalization",power->sign_normalization,&power->sign_normalization,&flg));

    752:   PetscOptionsHeadEnd();
    753:   PetscFunctionReturn(PETSC_SUCCESS);
    754: }

    756: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
    757: {
    758:   EPS_POWER *power = (EPS_POWER*)eps->data;

    760:   PetscFunctionBegin;
    761:   switch (shift) {
    762:     case EPS_POWER_SHIFT_CONSTANT:
    763:     case EPS_POWER_SHIFT_RAYLEIGH:
    764:     case EPS_POWER_SHIFT_WILKINSON:
    765:       if (power->shift_type != shift) {
    766:         power->shift_type = shift;
    767:         eps->state = EPS_STATE_INITIAL;
    768:       }
    769:       break;
    770:     default:
    771:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
    772:   }
    773:   PetscFunctionReturn(PETSC_SUCCESS);
    774: }

    776: /*@
    777:    EPSPowerSetShiftType - Sets the type of shifts used during the power
    778:    iteration. This can be used to emulate the Rayleigh Quotient Iteration
    779:    (RQI) method.

    781:    Logically Collective

    783:    Input Parameters:
    784: +  eps - the linear eigensolver context
    785: -  shift - the type of shift, see `EPSPowerShiftType` for possible values

    787:    Options Database Key:
    788: .  -eps_power_shift_type \<shift\> - sets the shift type, either `constant`, `rayleigh` or `wilkinson`

    790:    Notes:
    791:    By default, shifts are constant (`EPS_POWER_SHIFT_CONSTANT`) and the iteration
    792:    is the simple power method (or inverse iteration if a shift-and-invert
    793:    transformation is being used).

    795:    A variable shift can be specified (`EPS_POWER_SHIFT_RAYLEIGH` or
    796:    `EPS_POWER_SHIFT_WILKINSON`). In this case, the iteration behaves rather like
    797:    a cubic converging method such as RQI.

    799:    Details of the three variants can be found in {cite:p}`Her05`.

    801:    Level: advanced

    803: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetShiftType()`, `STSetShift()`, `EPSPowerShiftType`
    804: @*/
    805: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
    806: {
    807:   PetscFunctionBegin;
    810:   PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
    811:   PetscFunctionReturn(PETSC_SUCCESS);
    812: }

    814: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
    815: {
    816:   EPS_POWER *power = (EPS_POWER*)eps->data;

    818:   PetscFunctionBegin;
    819:   *shift = power->shift_type;
    820:   PetscFunctionReturn(PETSC_SUCCESS);
    821: }

    823: /*@
    824:    EPSPowerGetShiftType - Gets the type of shifts used during the power
    825:    iteration.

    827:    Not Collective

    829:    Input Parameter:
    830: .  eps - the linear eigensolver context

    832:    Output Parameter:
    833: .  shift - the type of shift

    835:    Level: advanced

    837: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetShiftType()`, `EPSPowerShiftType`
    838: @*/
    839: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
    840: {
    841:   PetscFunctionBegin;
    843:   PetscAssertPointer(shift,2);
    844:   PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
    845:   PetscFunctionReturn(PETSC_SUCCESS);
    846: }

    848: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
    849: {
    850:   EPS_POWER *power = (EPS_POWER*)eps->data;

    852:   PetscFunctionBegin;
    853:   if (power->nonlinear != nonlinear) {
    854:     power->nonlinear = nonlinear;
    855:     eps->useds = PetscNot(nonlinear);
    856:     eps->ops->setupsort = nonlinear? NULL: EPSSetUpSort_Default;
    857:     eps->state = EPS_STATE_INITIAL;
    858:   }
    859:   PetscFunctionReturn(PETSC_SUCCESS);
    860: }

    862: /*@
    863:    EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.

    865:    Logically Collective

    867:    Input Parameters:
    868: +  eps - the linear eigensolver context
    869: -  nonlinear - whether the problem is nonlinear or not

    871:    Options Database Key:
    872: .  -eps_power_nonlinear - sets the nonlinear flag

    874:    Notes:
    875:    If this flag is set, the solver assumes that the problem is nonlinear,
    876:    that is, the operators that define the eigenproblem are not constant
    877:    matrices, but depend on the eigenvector, $A(x)x=\lambda B(x)x$. This is
    878:    different from the case of nonlinearity with respect to the eigenvalue
    879:    (use the `NEP` solver class for this kind of problems).

    881:    The way in which nonlinear operators are specified is very similar to
    882:    the case of PETSc's `SNES` solver. The difference is that the callback
    883:    functions are provided via composed functions `formFunction` and
    884:    `formJacobian` in each of the matrix objects passed as arguments of
    885:    `EPSSetOperators()`. The application context required for these functions
    886:    can be attached via a composed `PetscContainer`.

    888:    Level: advanced

    890: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetNonlinear()`, `EPSSetOperators()`
    891: @*/
    892: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
    893: {
    894:   PetscFunctionBegin;
    897:   PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
    898:   PetscFunctionReturn(PETSC_SUCCESS);
    899: }

    901: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
    902: {
    903:   EPS_POWER *power = (EPS_POWER*)eps->data;

    905:   PetscFunctionBegin;
    906:   *nonlinear = power->nonlinear;
    907:   PetscFunctionReturn(PETSC_SUCCESS);
    908: }

    910: /*@
    911:    EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.

    913:    Not Collective

    915:    Input Parameter:
    916: .  eps - the linear eigensolver context

    918:    Output Parameter:
    919: .  nonlinear - the nonlinear flag

    921:    Level: advanced

    923: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetUpdate()`, `EPSPowerSetNonlinear()`
    924: @*/
    925: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
    926: {
    927:   PetscFunctionBegin;
    929:   PetscAssertPointer(nonlinear,2);
    930:   PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
    931:   PetscFunctionReturn(PETSC_SUCCESS);
    932: }

    934: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
    935: {
    936:   EPS_POWER *power = (EPS_POWER*)eps->data;

    938:   PetscFunctionBegin;
    939:   PetscCheck(power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
    940:   power->update = update;
    941:   eps->state = EPS_STATE_INITIAL;
    942:   PetscFunctionReturn(PETSC_SUCCESS);
    943: }

    945: /*@
    946:    EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
    947:    for nonlinear problems. This potentially has a better convergence.

    949:    Logically Collective

    951:    Input Parameters:
    952: +  eps - the linear eigensolver context
    953: -  update - whether the residual is updated monolithically or not

    955:    Options Database Key:
    956: .  -eps_power_update - sets the update flag

    958:    Note:
    959:    This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.

    961:    Level: advanced

    963: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetUpdate()`, `EPSPowerGetNonlinear()`, `EPSSetOperators()`
    964: @*/
    965: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
    966: {
    967:   PetscFunctionBegin;
    970:   PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
    971:   PetscFunctionReturn(PETSC_SUCCESS);
    972: }

    974: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
    975: {
    976:   EPS_POWER *power = (EPS_POWER*)eps->data;

    978:   PetscFunctionBegin;
    979:   *update = power->update;
    980:   PetscFunctionReturn(PETSC_SUCCESS);
    981: }

    983: /*@
    984:    EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
    985:    for nonlinear problems.

    987:    Not Collective

    989:    Input Parameter:
    990: .  eps - the linear eigensolver context

    992:    Output Parameter:
    993: .  update - the update flag

    995:    Level: advanced

    997: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetUpdate()`, `EPSPowerSetNonlinear()`
    998: @*/
    999: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
    1000: {
    1001:   PetscFunctionBegin;
    1003:   PetscAssertPointer(update,2);
    1004:   PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
    1005:   PetscFunctionReturn(PETSC_SUCCESS);
    1006: }

    1008: static PetscErrorCode EPSPowerSetSignNormalization_Power(EPS eps,PetscBool sign_normalization)
    1009: {
    1010:   EPS_POWER *power = (EPS_POWER*)eps->data;

    1012:   PetscFunctionBegin;
    1013:   power->sign_normalization = sign_normalization;
    1014:   PetscFunctionReturn(PETSC_SUCCESS);
    1015: }

    1017: /*@
    1018:    EPSPowerSetSignNormalization - Sets a flag to indicate whether the $Bx$ vector
    1019:    should be normalized by the sign of the first non-zero element in the vector.
    1020:    E.g., if this is true, the post-normalization value of the first non-zero element
    1021:    in the vector is guaranteed to be positive.

    1023:    Logically Collective

    1025:    Input Parameters:
    1026: +  eps                - the linear eigensolver context
    1027: -  sign_normalization - whether $Bx$ should be multiplied by the sign of the first non-zero
    1028:                         element when performing normalization steps

    1030:    Options Database Key:
    1031: .  -eps_power_sign_normalization - sets the sign normalization flag

    1033:    Note:
    1034:    This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.

    1036:    Level: advanced

    1038: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetNonlinear()`, `EPSPowerGetSignNormalization()`
    1039: @*/
    1040: PetscErrorCode EPSPowerSetSignNormalization(EPS eps,PetscBool sign_normalization)
    1041: {
    1042:   PetscFunctionBegin;
    1045:   PetscTryMethod(eps,"EPSPowerSetSignNormalization_C",(EPS,PetscBool),(eps,sign_normalization));
    1046:   PetscFunctionReturn(PETSC_SUCCESS);
    1047: }

    1049: static PetscErrorCode EPSPowerGetSignNormalization_Power(EPS eps,PetscBool *sign_normalization)
    1050: {
    1051:   EPS_POWER *power = (EPS_POWER*)eps->data;

    1053:   PetscFunctionBegin;
    1054:   *sign_normalization = power->sign_normalization;
    1055:   PetscFunctionReturn(PETSC_SUCCESS);
    1056: }

    1058: /*@
    1059:    EPSPowerGetSignNormalization - Returns a flag indicating whether the $Bx$ vector
    1060:    is normalized by the sign of the first non-zero element in the vector. E.g.,
    1061:    if this is true, the post-normalization value of the first non-zero element in
    1062:    the vector is guaranteed to be positive.

    1064:    Not Collective

    1066:    Input Parameter:
    1067: .  eps - the linear eigensolver context

    1069:    Output Parameter:
    1070: .  sign_normalization - the sign normalization flag

    1072:    Level: advanced

    1074: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetSignNormalization()`
    1075: @*/
    1076: PetscErrorCode EPSPowerGetSignNormalization(EPS eps,PetscBool *sign_normalization)
    1077: {
    1078:   PetscFunctionBegin;
    1080:   PetscAssertPointer(sign_normalization,2);
    1081:   PetscUseMethod(eps,"EPSPowerGetSignNormalization_C",(EPS,PetscBool*),(eps,sign_normalization));
    1082:   PetscFunctionReturn(PETSC_SUCCESS);
    1083: }

    1085: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
    1086: {
    1087:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    1089:   PetscFunctionBegin;
    1090:   PetscCall(PetscObjectReference((PetscObject)snes));
    1091:   PetscCall(SNESDestroy(&power->snes));
    1092:   power->snes = snes;
    1093:   eps->state = EPS_STATE_INITIAL;
    1094:   PetscFunctionReturn(PETSC_SUCCESS);
    1095: }

    1097: /*@
    1098:    EPSPowerSetSNES - Associate a nonlinear solver object (`SNES`) to the
    1099:    eigenvalue solver (to be used in nonlinear inverse iteration).

    1101:    Collective

    1103:    Input Parameters:
    1104: +  eps  - the linear eigensolver context
    1105: -  snes - the nonlinear solver object

    1107:    Note:
    1108:    This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.

    1110:    Level: advanced

    1112: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetNonlinear()`, `EPSPowerGetSNES()`
    1113: @*/
    1114: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
    1115: {
    1116:   PetscFunctionBegin;
    1119:   PetscCheckSameComm(eps,1,snes,2);
    1120:   PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
    1121:   PetscFunctionReturn(PETSC_SUCCESS);
    1122: }

    1124: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
    1125: {
    1126:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    1128:   PetscFunctionBegin;
    1129:   if (!power->snes) {
    1130:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes));
    1131:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1));
    1132:     PetscCall(SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix));
    1133:     PetscCall(SNESAppendOptionsPrefix(power->snes,"eps_power_"));
    1134:     PetscCall(PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options));
    1135:   }
    1136:   *snes = power->snes;
    1137:   PetscFunctionReturn(PETSC_SUCCESS);
    1138: }

    1140: /*@
    1141:    EPSPowerGetSNES - Retrieve the nonlinear solver object (`SNES`) associated
    1142:    with the eigenvalue solver.

    1144:    Not Collective

    1146:    Input Parameter:
    1147: .  eps - the linear eigensolver context

    1149:    Output Parameter:
    1150: .  snes - the nonlinear solver object

    1152:    Level: advanced

    1154: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetSNES()`
    1155: @*/
    1156: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
    1157: {
    1158:   PetscFunctionBegin;
    1160:   PetscAssertPointer(snes,2);
    1161:   PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
    1162:   PetscFunctionReturn(PETSC_SUCCESS);
    1163: }

    1165: static PetscErrorCode EPSReset_Power(EPS eps)
    1166: {
    1167:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    1169:   PetscFunctionBegin;
    1170:   if (power->snes) PetscCall(SNESReset(power->snes));
    1171:   PetscFunctionReturn(PETSC_SUCCESS);
    1172: }

    1174: static PetscErrorCode EPSDestroy_Power(EPS eps)
    1175: {
    1176:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    1178:   PetscFunctionBegin;
    1179:   if (power->nonlinear) PetscCall(SNESDestroy(&power->snes));
    1180:   PetscCall(PetscFree(eps->data));
    1181:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL));
    1182:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL));
    1183:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL));
    1184:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL));
    1185:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL));
    1186:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL));
    1187:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",NULL));
    1188:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",NULL));
    1189:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL));
    1190:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL));
    1191:   PetscFunctionReturn(PETSC_SUCCESS);
    1192: }

    1194: static PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
    1195: {
    1196:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    1197:   PetscBool      isascii;

    1199:   PetscFunctionBegin;
    1200:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1201:   if (isascii) {
    1202:     if (power->nonlinear) {
    1203:       if (power->sign_normalization) PetscCall(PetscViewerASCIIPrintf(viewer,"  normalizing Bx by the sign of the first nonzero element\n"));
    1204:       else PetscCall(PetscViewerASCIIPrintf(viewer,"  not normalizing Bx by the sign of the first nonzero element\n"));
    1205:       PetscCall(PetscViewerASCIIPrintf(viewer,"  using nonlinear inverse iteration\n"));
    1206:       if (power->update) PetscCall(PetscViewerASCIIPrintf(viewer,"  updating the residual monolithically\n"));
    1207:       if (!power->snes) PetscCall(EPSPowerGetSNES(eps,&power->snes));
    1208:       PetscCall(PetscViewerASCIIPushTab(viewer));
    1209:       PetscCall(SNESView(power->snes,viewer));
    1210:       PetscCall(PetscViewerASCIIPopTab(viewer));
    1211:     } else PetscCall(PetscViewerASCIIPrintf(viewer,"  %s shifts\n",EPSPowerShiftTypes[power->shift_type]));
    1212:   }
    1213:   PetscFunctionReturn(PETSC_SUCCESS);
    1214: }

    1216: static PetscErrorCode EPSComputeVectors_Power(EPS eps)
    1217: {
    1218:   EPS_POWER      *power = (EPS_POWER*)eps->data;

    1220:   PetscFunctionBegin;
    1221:   if (eps->twosided) {
    1222:     PetscCall(EPSComputeVectors_Twosided(eps));
    1223:     PetscCall(BVNormalize(eps->V,NULL));
    1224:     PetscCall(BVNormalize(eps->W,NULL));
    1225:   } else if (!power->nonlinear) PetscCall(EPSComputeVectors_Schur(eps));
    1226:   PetscFunctionReturn(PETSC_SUCCESS);
    1227: }

    1229: static PetscErrorCode EPSSetDefaultST_Power(EPS eps)
    1230: {
    1231:   EPS_POWER      *power = (EPS_POWER*)eps->data;
    1232:   KSP            ksp;
    1233:   PC             pc;

    1235:   PetscFunctionBegin;
    1236:   if (power->nonlinear) {
    1237:     eps->categ=EPS_CATEGORY_PRECOND;
    1238:     PetscCall(STGetKSP(eps->st,&ksp));
    1239:     /* Set ST as STPRECOND so it can carry one preconditioning matrix
    1240:      * It is useful when A and B are shell matrices
    1241:      */
    1242:     PetscCall(STSetType(eps->st,STPRECOND));
    1243:     PetscCall(KSPGetPC(ksp,&pc));
    1244:     PetscCall(PCSetType(pc,PCNONE));
    1245:   }
    1246:   PetscFunctionReturn(PETSC_SUCCESS);
    1247: }

    1249: /*MC
    1250:    EPSPOWER - EPSPOWER = "power" - The simple power iteration and inverse
    1251:    iteration.

    1253:    Notes:
    1254:    This solver is very basic and is not recommended in general, since it
    1255:    will not be competitive with respect to other solvers.

    1257:    The implemented method is the power iteration, or inverse iteration in
    1258:    case the selected spectral transformation is `STSINVERT`. In this latter
    1259:    case it is possible to use a dynamic shift, as in the RQI method, see
    1260:    `EPSPowerSetShiftType()`.

    1262:    The solver incorporates deflation so that several eigenpairs can be
    1263:    computed. There is also a two-sided implementation that also computes
    1264:    left eigenvectors.

    1266:    This solver can also be used for nonlinear inverse iteration on the problem
    1267:    $A(x)x=\lambda B(x)x$, where $A$ and $B$ are not constant matrices but
    1268:    depend on the eigenvector $x$. This mode is enabled with `EPSPowerSetNonlinear()`.
    1269:    Note that this is a nonlinear eigenvector problem, as opposed to problems
    1270:    addressed in `NEP` that are nonlinear with respect to the eigenvalue.

    1272:    Level: beginner

    1274: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSSetTwoSided()`, `EPSPowerSetShiftType()`, `EPSPowerSetNonlinear()`
    1275: M*/
    1276: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
    1277: {
    1278:   EPS_POWER      *ctx;

    1280:   PetscFunctionBegin;
    1281:   PetscCall(PetscNew(&ctx));
    1282:   eps->data = (void*)ctx;

    1284:   eps->useds = PETSC_TRUE;
    1285:   eps->categ = EPS_CATEGORY_OTHER;

    1287:   eps->ops->setup          = EPSSetUp_Power;
    1288:   eps->ops->setupsort      = EPSSetUpSort_Default;
    1289:   eps->ops->setfromoptions = EPSSetFromOptions_Power;
    1290:   eps->ops->reset          = EPSReset_Power;
    1291:   eps->ops->destroy        = EPSDestroy_Power;
    1292:   eps->ops->view           = EPSView_Power;
    1293:   eps->ops->backtransform  = EPSBackTransform_Power;
    1294:   eps->ops->computevectors = EPSComputeVectors_Power;
    1295:   eps->ops->setdefaultst   = EPSSetDefaultST_Power;
    1296:   eps->stopping            = EPSStopping_Power;
    1297:   ctx->sign_normalization  = PETSC_TRUE;

    1299:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power));
    1300:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power));
    1301:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power));
    1302:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power));
    1303:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power));
    1304:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power));
    1305:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",EPSPowerSetSignNormalization_Power));
    1306:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",EPSPowerGetSignNormalization_Power));
    1307:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power));
    1308:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power));
    1309:   PetscFunctionReturn(PETSC_SUCCESS);
    1310: }