Actual source code: power.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

480:     } else {  /* RQI */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

736:   PetscFunctionBegin;

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

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

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

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

751:   PetscFunctionReturn(PETSC_SUCCESS);
752: }

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

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

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

779:    Logically Collective

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

785:    Options Database Key:
786: .  -eps_power_shift_type - Sets the shift type (either 'constant' or
787:                            'rayleigh' or 'wilkinson')

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

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

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

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

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

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

824:    Not Collective

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

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

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

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

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

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

862:    Logically Collective

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

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

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

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

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

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

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

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

910:    Not Collective

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

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

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

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

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

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

946:    Logically Collective

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

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

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

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

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

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

981:    Not Collective

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

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

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

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

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

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

1017:    Logically Collective

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

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

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

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

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

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

1055:    Not Collective

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

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

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

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

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

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

1092:    Collective

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

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

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

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

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

1132:    Not Collective

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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