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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

481:     } else {  /* RQI */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

780:    Logically Collective

782:    Input Parameters:
783: +  eps - the eigenproblem solver context
784: -  shift - the type of shift

786:    Options Database Key:
787: .  -eps_power_shift_type - Sets the shift type (either 'constant' or
788:                            '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:    Level: advanced

801: .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
802: @*/
803: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
804: {
805:   PetscFunctionBegin;
808:   PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
809:   PetscFunctionReturn(PETSC_SUCCESS);
810: }

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

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

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

825:    Not Collective

827:    Input Parameter:
828: .  eps - the eigenproblem solver context

830:    Output Parameter:
831: .  shift - the type of shift

833:    Level: advanced

835: .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
836: @*/
837: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
838: {
839:   PetscFunctionBegin;
841:   PetscAssertPointer(shift,2);
842:   PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
843:   PetscFunctionReturn(PETSC_SUCCESS);
844: }

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

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

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

863:    Logically Collective

865:    Input Parameters:
866: +  eps - the eigenproblem solver context
867: -  nonlinear - whether the problem is nonlinear or not

869:    Options Database Key:
870: .  -eps_power_nonlinear - Sets the nonlinear flag

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

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

886:    Level: advanced

888: .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
889: @*/
890: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
891: {
892:   PetscFunctionBegin;
895:   PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
896:   PetscFunctionReturn(PETSC_SUCCESS);
897: }

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

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

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

911:    Not Collective

913:    Input Parameter:
914: .  eps - the eigenproblem solver context

916:    Output Parameter:
917: .  nonlinear - the nonlinear flag

919:    Level: advanced

921: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
922: @*/
923: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
924: {
925:   PetscFunctionBegin;
927:   PetscAssertPointer(nonlinear,2);
928:   PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
929:   PetscFunctionReturn(PETSC_SUCCESS);
930: }

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

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

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

947:    Logically Collective

949:    Input Parameters:
950: +  eps - the eigenproblem solver context
951: -  update - whether the residual is updated monolithically or not

953:    Options Database Key:
954: .  -eps_power_update - Sets the update flag

956:    Level: advanced

958: .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
959: @*/
960: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
961: {
962:   PetscFunctionBegin;
965:   PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
970: {
971:   EPS_POWER *power = (EPS_POWER*)eps->data;

973:   PetscFunctionBegin;
974:   *update = power->update;
975:   PetscFunctionReturn(PETSC_SUCCESS);
976: }

978: /*@
979:    EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
980:    for nonlinear problems.

982:    Not Collective

984:    Input Parameter:
985: .  eps - the eigenproblem solver context

987:    Output Parameter:
988: .  update - the update flag

990:    Level: advanced

992: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
993: @*/
994: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
995: {
996:   PetscFunctionBegin;
998:   PetscAssertPointer(update,2);
999:   PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
1000:   PetscFunctionReturn(PETSC_SUCCESS);
1001: }

1003: static PetscErrorCode EPSPowerSetSignNormalization_Power(EPS eps,PetscBool sign_normalization)
1004: {
1005:   EPS_POWER *power = (EPS_POWER*)eps->data;

1007:   PetscFunctionBegin;
1008:   power->sign_normalization = sign_normalization;
1009:   PetscFunctionReturn(PETSC_SUCCESS);
1010: }

1012: /*@
1013:    EPSPowerSetSignNormalization - Sets a flag to indicate whether the Bx vector
1014:    should be normalized by the sign of the first non-zero element in the vector.
1015:    E.g., if this is true, the post-normalization value of the first non-zero element
1016:    in the vector is guaranteed to be positive.

1018:    Logically Collective

1020:    Input Parameters:
1021: +  eps - the eigenproblem solver context
1022: -  sign_normalization - whether Bx should be multiplied by the sign of the first non-zero
1023:                         element when performing normalization steps

1025:    Options Database Key:
1026: .  -eps_power_sign_normalization - Sets the sign normalization flag

1028:    Level: advanced

1030: .seealso: EPSPowerGetSignNormalization()
1031: @*/
1032: PetscErrorCode EPSPowerSetSignNormalization(EPS eps,PetscBool sign_normalization)
1033: {
1034:   PetscFunctionBegin;
1037:   PetscTryMethod(eps,"EPSPowerSetSignNormalization_C",(EPS,PetscBool),(eps,sign_normalization));
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

1041: static PetscErrorCode EPSPowerGetSignNormalization_Power(EPS eps,PetscBool *sign_normalization)
1042: {
1043:   EPS_POWER *power = (EPS_POWER*)eps->data;

1045:   PetscFunctionBegin;
1046:   *sign_normalization = power->sign_normalization;
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: /*@
1051:    EPSPowerGetSignNormalization - Returns a flag indicating whether the Bx vector
1052:    is normalized by the sign of the first non-zero element in the vector. E.g.,
1053:    if this is true, the post-normalization value of the first non-zero element in
1054:    the vector is guaranteed to be positive.

1056:    Not Collective

1058:    Input Parameter:
1059: .  eps - the eigenproblem solver context

1061:    Output Parameter:
1062: .  sign_normalization - the sign normalization flag

1064:    Level: advanced

1066: .seealso: EPSPowerSetSignNormalization()
1067: @*/
1068: PetscErrorCode EPSPowerGetSignNormalization(EPS eps,PetscBool *sign_normalization)
1069: {
1070:   PetscFunctionBegin;
1072:   PetscAssertPointer(sign_normalization,2);
1073:   PetscUseMethod(eps,"EPSPowerGetSignNormalization_C",(EPS,PetscBool*),(eps,sign_normalization));
1074:   PetscFunctionReturn(PETSC_SUCCESS);
1075: }

1077: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
1078: {
1079:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1081:   PetscFunctionBegin;
1082:   PetscCall(PetscObjectReference((PetscObject)snes));
1083:   PetscCall(SNESDestroy(&power->snes));
1084:   power->snes = snes;
1085:   eps->state = EPS_STATE_INITIAL;
1086:   PetscFunctionReturn(PETSC_SUCCESS);
1087: }

1089: /*@
1090:    EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
1091:    eigenvalue solver (to be used in nonlinear inverse iteration).

1093:    Collective

1095:    Input Parameters:
1096: +  eps  - the eigenvalue solver
1097: -  snes - the nonlinear solver object

1099:    Level: advanced

1101: .seealso: EPSPowerGetSNES()
1102: @*/
1103: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1104: {
1105:   PetscFunctionBegin;
1108:   PetscCheckSameComm(eps,1,snes,2);
1109:   PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1110:   PetscFunctionReturn(PETSC_SUCCESS);
1111: }

1113: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1114: {
1115:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1117:   PetscFunctionBegin;
1118:   if (!power->snes) {
1119:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes));
1120:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1));
1121:     PetscCall(SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix));
1122:     PetscCall(SNESAppendOptionsPrefix(power->snes,"eps_power_"));
1123:     PetscCall(PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options));
1124:   }
1125:   *snes = power->snes;
1126:   PetscFunctionReturn(PETSC_SUCCESS);
1127: }

1129: /*@
1130:    EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1131:    with the eigenvalue solver.

1133:    Not Collective

1135:    Input Parameter:
1136: .  eps - the eigenvalue solver

1138:    Output Parameter:
1139: .  snes - the nonlinear solver object

1141:    Level: advanced

1143: .seealso: EPSPowerSetSNES()
1144: @*/
1145: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1146: {
1147:   PetscFunctionBegin;
1149:   PetscAssertPointer(snes,2);
1150:   PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1151:   PetscFunctionReturn(PETSC_SUCCESS);
1152: }

1154: static PetscErrorCode EPSReset_Power(EPS eps)
1155: {
1156:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1158:   PetscFunctionBegin;
1159:   if (power->snes) PetscCall(SNESReset(power->snes));
1160:   PetscFunctionReturn(PETSC_SUCCESS);
1161: }

1163: static PetscErrorCode EPSDestroy_Power(EPS eps)
1164: {
1165:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1167:   PetscFunctionBegin;
1168:   if (power->nonlinear) PetscCall(SNESDestroy(&power->snes));
1169:   PetscCall(PetscFree(eps->data));
1170:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL));
1171:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL));
1172:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL));
1173:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL));
1174:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL));
1175:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL));
1176:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",NULL));
1177:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",NULL));
1178:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL));
1179:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL));
1180:   PetscFunctionReturn(PETSC_SUCCESS);
1181: }

1183: static PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1184: {
1185:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1186:   PetscBool      isascii;

1188:   PetscFunctionBegin;
1189:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1190:   if (isascii) {
1191:     if (power->nonlinear) {
1192:       if (power->sign_normalization) PetscCall(PetscViewerASCIIPrintf(viewer,"  normalizing Bx by the sign of the first nonzero element\n"));
1193:       else PetscCall(PetscViewerASCIIPrintf(viewer,"  not normalizing Bx by the sign of the first nonzero element\n"));
1194:       PetscCall(PetscViewerASCIIPrintf(viewer,"  using nonlinear inverse iteration\n"));
1195:       if (power->update) PetscCall(PetscViewerASCIIPrintf(viewer,"  updating the residual monolithically\n"));
1196:       if (!power->snes) PetscCall(EPSPowerGetSNES(eps,&power->snes));
1197:       PetscCall(PetscViewerASCIIPushTab(viewer));
1198:       PetscCall(SNESView(power->snes,viewer));
1199:       PetscCall(PetscViewerASCIIPopTab(viewer));
1200:     } else PetscCall(PetscViewerASCIIPrintf(viewer,"  %s shifts\n",EPSPowerShiftTypes[power->shift_type]));
1201:   }
1202:   PetscFunctionReturn(PETSC_SUCCESS);
1203: }

1205: static PetscErrorCode EPSComputeVectors_Power(EPS eps)
1206: {
1207:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1209:   PetscFunctionBegin;
1210:   if (eps->twosided) {
1211:     PetscCall(EPSComputeVectors_Twosided(eps));
1212:     PetscCall(BVNormalize(eps->V,NULL));
1213:     PetscCall(BVNormalize(eps->W,NULL));
1214:   } else if (!power->nonlinear) PetscCall(EPSComputeVectors_Schur(eps));
1215:   PetscFunctionReturn(PETSC_SUCCESS);
1216: }

1218: static PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1219: {
1220:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1221:   KSP            ksp;
1222:   PC             pc;

1224:   PetscFunctionBegin;
1225:   if (power->nonlinear) {
1226:     eps->categ=EPS_CATEGORY_PRECOND;
1227:     PetscCall(STGetKSP(eps->st,&ksp));
1228:     /* Set ST as STPRECOND so it can carry one preconditioning matrix
1229:      * It is useful when A and B are shell matrices
1230:      */
1231:     PetscCall(STSetType(eps->st,STPRECOND));
1232:     PetscCall(KSPGetPC(ksp,&pc));
1233:     PetscCall(PCSetType(pc,PCNONE));
1234:   }
1235:   PetscFunctionReturn(PETSC_SUCCESS);
1236: }

1238: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1239: {
1240:   EPS_POWER      *ctx;

1242:   PetscFunctionBegin;
1243:   PetscCall(PetscNew(&ctx));
1244:   eps->data = (void*)ctx;

1246:   eps->useds = PETSC_TRUE;
1247:   eps->categ = EPS_CATEGORY_OTHER;

1249:   eps->ops->setup          = EPSSetUp_Power;
1250:   eps->ops->setupsort      = EPSSetUpSort_Default;
1251:   eps->ops->setfromoptions = EPSSetFromOptions_Power;
1252:   eps->ops->reset          = EPSReset_Power;
1253:   eps->ops->destroy        = EPSDestroy_Power;
1254:   eps->ops->view           = EPSView_Power;
1255:   eps->ops->backtransform  = EPSBackTransform_Power;
1256:   eps->ops->computevectors = EPSComputeVectors_Power;
1257:   eps->ops->setdefaultst   = EPSSetDefaultST_Power;
1258:   eps->stopping            = EPSStopping_Power;
1259:   ctx->sign_normalization  = PETSC_TRUE;

1261:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power));
1262:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power));
1263:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power));
1264:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power));
1265:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power));
1266:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power));
1267:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",EPSPowerSetSignNormalization_Power));
1268:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",EPSPowerGetSignNormalization_Power));
1269:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power));
1270:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power));
1271:   PetscFunctionReturn(PETSC_SUCCESS);
1272: }