Actual source code: power.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "power"
13: Method: Power Iteration
15: Algorithm:
17: This solver implements the power iteration for finding dominant
18: eigenpairs. It also includes the following well-known methods:
19: - Inverse Iteration: when used in combination with shift-and-invert
20: spectral transformation.
21: - Rayleigh Quotient Iteration (RQI): also with shift-and-invert plus
22: a variable shift.
24: It can also be used for nonlinear inverse iteration on the problem
25: A(x)*x=lambda*B(x)*x, where A and B are not constant but depend on x.
27: References:
29: [1] "Single Vector Iteration Methods in SLEPc", SLEPc Technical Report
30: STR-2, available at https://slepc.upv.es.
31: */
33: #include <slepc/private/epsimpl.h>
34: #include <slepcblaslapack.h>
35: /* PETSc headers */
36: #include <petscdm.h>
37: #include <petscsnes.h>
39: static PetscErrorCode EPSPowerFormFunction_Update(SNES,Vec,Vec,void*);
40: static PetscErrorCode EPSSolve_Power(EPS);
41: static PetscErrorCode EPSSolve_TS_Power(EPS);
43: typedef struct {
44: EPSPowerShiftType shift_type;
45: PetscBool nonlinear;
46: PetscBool update;
47: PetscBool sign_normalization;
48: SNES snes;
49: PetscErrorCode (*formFunctionB)(SNES,Vec,Vec,void*);
50: void *formFunctionBctx;
51: PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
52: void *formFunctionActx;
53: PetscErrorCode (*formFunctionAB)(SNES,Vec,Vec,Vec,void*);
54: PetscErrorCode (*formNorm)(SNES,Vec,PetscReal*,void*);
55: void *formNormCtx;
56: PetscInt idx; /* index of the first nonzero entry in the iteration vector */
57: PetscMPIInt p; /* process id of the owner of idx */
58: PetscReal norm0; /* norm of initial vector */
59: } EPS_POWER;
61: static PetscErrorCode SNESMonitor_PowerUpdate(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
62: {
63: EPS eps;
65: PetscFunctionBegin;
66: PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
67: PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
68: /* Call EPS monitor on each SNES iteration */
69: PetscCall(EPSMonitor(eps,its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode EPSSetUp_Power(EPS eps)
74: {
75: EPS_POWER *power = (EPS_POWER*)eps->data;
76: STMatMode mode;
77: Mat A,B,P;
78: Vec res;
79: PetscContainer container;
80: PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
81: PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
82: PetscErrorCode (*formNorm)(SNES,Vec,PetscReal*,void*);
83: void *ctx;
85: PetscFunctionBegin;
86: EPSCheckNotStructured(eps);
87: if (eps->nev==0) eps->nev = 1;
88: if (eps->ncv!=PETSC_DETERMINE) {
89: PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
90: } else eps->ncv = eps->nev;
91: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
92: if (eps->max_it==PETSC_DETERMINE) {
93: /* SNES will directly return the solution for us, and we need to do only one iteration */
94: if (power->nonlinear && power->update) eps->max_it = 1;
95: else eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
96: }
97: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
98: PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
99: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
100: PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
101: EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," (with variable shifts)");
102: PetscCall(STGetMatMode(eps->st,&mode));
103: PetscCheck(mode!=ST_MATMODE_INPLACE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
104: }
105: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_THRESHOLD);
106: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
107: PetscCall(EPSAllocateSolution(eps,0));
108: PetscCall(EPS_SetInnerProduct(eps));
110: if (power->nonlinear) {
111: PetscCheck(eps->nev==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
112: PetscCall(EPSSetWorkVecs(eps,3));
113: PetscCheck(!power->update || eps->max_it==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"More than one iteration is not allowed for Newton eigensolver (SNES)");
115: /* Set up SNES solver */
116: /* If SNES was setup, we need to reset it so that it will be consistent with the current EPS */
117: if (power->snes) PetscCall(SNESReset(power->snes));
118: else PetscCall(EPSPowerGetSNES(eps,&power->snes));
120: /* For nonlinear solver (Newton), we should scale the initial vector back.
121: The initial vector will be scaled in EPSSetUp. */
122: if (eps->IS) PetscCall(VecNorm(eps->IS[0],NORM_2,&power->norm0));
124: PetscCall(EPSGetOperators(eps,&A,&B));
126: /* form function */
127: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA));
128: PetscCheck(formFunctionA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
129: PetscCall(PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container));
130: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
131: else ctx = NULL;
132: PetscCall(MatCreateVecs(A,&res,NULL));
133: power->formFunctionA = formFunctionA;
134: power->formFunctionActx = ctx;
135: if (power->update) {
136: PetscCall(SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx));
137: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB));
138: PetscCall(SNESMonitorSet(power->snes,SNESMonitor_PowerUpdate,NULL,NULL));
139: }
140: else PetscCall(SNESSetFunction(power->snes,res,formFunctionA,ctx));
141: PetscCall(VecDestroy(&res));
143: /* form Jacobian */
144: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA));
145: PetscCheck(formJacobianA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
146: PetscCall(PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container));
147: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
148: else ctx = NULL;
149: /* This allows users to compute a different preconditioning matrix than A.
150: * It is useful when A and B are shell matrices.
151: */
152: PetscCall(STGetPreconditionerMat(eps->st,&P));
153: /* If users do not provide a matrix, we simply use A */
154: PetscCall(SNESSetJacobian(power->snes,A,P? P:A,formJacobianA,ctx));
155: PetscCall(SNESSetFromOptions(power->snes));
156: PetscCall(SNESSetUp(power->snes));
158: ctx = NULL;
159: formNorm = NULL;
160: if (B) {
161: PetscCall(PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB));
162: PetscCall(PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container));
163: if (power->formFunctionB && container) PetscCall(PetscContainerGetPointer(container,&power->formFunctionBctx));
164: else power->formFunctionBctx = NULL;
166: /* form norm */
167: PetscCall(PetscObjectQueryFunction((PetscObject)B,"formNorm",&formNorm));
168: if (formNorm) {
169: PetscCall(PetscObjectQuery((PetscObject)B,"formNormCtx",(PetscObject*)&container));
170: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
171: }
172: }
173: power->formNorm = formNorm;
174: power->formNormCtx = ctx;
175: } else {
176: if (eps->twosided) PetscCall(EPSSetWorkVecs(eps,3));
177: else PetscCall(EPSSetWorkVecs(eps,2));
178: PetscCall(DSSetType(eps->ds,DSNHEP));
179: PetscCall(DSAllocate(eps->ds,eps->nev));
180: }
181: /* dispatch solve method */
182: if (eps->twosided) {
183: PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
184: PetscCheck(power->shift_type!=EPS_POWER_SHIFT_WILKINSON,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
185: eps->ops->solve = EPSSolve_TS_Power;
186: } else eps->ops->solve = EPSSolve_Power;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*
191: Find the index of the first nonzero entry of x, and the process that owns it.
192: */
193: static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscMPIInt *p)
194: {
195: PetscInt i,first,last,N;
196: PetscLayout map;
197: const PetscScalar *xx;
199: PetscFunctionBegin;
200: PetscCall(VecGetSize(x,&N));
201: PetscCall(VecGetOwnershipRange(x,&first,&last));
202: PetscCall(VecGetArrayRead(x,&xx));
203: for (i=first;i<last;i++) {
204: if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
205: }
206: PetscCall(VecRestoreArrayRead(x,&xx));
207: if (i==last) i=N;
208: PetscCallMPI(MPIU_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x)));
209: PetscCheck(*idx!=N,PetscObjectComm((PetscObject)x),PETSC_ERR_PLIB,"Zero vector found");
210: PetscCall(VecGetLayout(x,&map));
211: PetscCall(PetscLayoutFindOwner(map,*idx,p));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*
216: Normalize a vector x with respect to a given norm as well as, optionally, the
217: sign of the first nonzero entry (assumed to be idx in process p).
218: */
219: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscMPIInt p,PetscBool sign_normalization,PetscScalar *sign)
220: {
221: PetscScalar alpha=1.0;
222: PetscInt first,last;
223: PetscReal absv;
224: const PetscScalar *xx;
226: PetscFunctionBegin;
227: if (sign_normalization) {
228: PetscCall(VecGetOwnershipRange(x,&first,&last));
229: if (idx>=first && idx<last) {
230: PetscCall(VecGetArrayRead(x,&xx));
231: absv = PetscAbsScalar(xx[idx-first]);
232: if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
233: PetscCall(VecRestoreArrayRead(x,&xx));
234: }
235: PetscCallMPI(MPI_Bcast(&alpha,1,MPIU_SCALAR,p,PetscObjectComm((PetscObject)x)));
236: }
237: if (sign) *sign = alpha;
238: alpha *= norm;
239: PetscCall(VecScale(x,1.0/alpha));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
244: {
245: EPS_POWER *power = (EPS_POWER*)eps->data;
246: Mat B;
248: PetscFunctionBegin;
249: PetscCall(STResetMatrixState(eps->st));
250: PetscCall(EPSGetOperators(eps,NULL,&B));
251: if (B) {
252: if (power->formFunctionB) PetscCall((*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx));
253: else PetscCall(MatMult(B,x,Bx));
254: } else PetscCall(VecCopy(x,Bx));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
259: {
260: EPS_POWER *power = (EPS_POWER*)eps->data;
261: Mat A;
263: PetscFunctionBegin;
264: PetscCall(STResetMatrixState(eps->st));
265: PetscCall(EPSGetOperators(eps,&A,NULL));
266: PetscCheck(A,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
267: if (power->formFunctionA) PetscCall((*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx));
268: else PetscCall(MatMult(A,x,Ax));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
273: {
274: EPS eps;
275: PetscReal bx;
276: Vec Bx;
277: PetscScalar sign;
278: EPS_POWER *power;
280: PetscFunctionBegin;
281: PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
282: PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
283: PetscCall(STResetMatrixState(eps->st));
284: power = (EPS_POWER*)eps->data;
285: Bx = eps->work[2];
286: if (power->formFunctionAB) PetscCall((*power->formFunctionAB)(snes,x,y,Bx,ctx));
287: else {
288: PetscCall(EPSPowerUpdateFunctionA(eps,x,y));
289: PetscCall(EPSPowerUpdateFunctionB(eps,x,Bx));
290: }
291: if (power->formNorm) PetscCall((*power->formNorm)(snes,Bx,&bx,power->formNormCtx));
292: else PetscCall(VecNorm(Bx,NORM_2,&bx));
293: PetscCall(Normalize(Bx,bx,power->idx,power->p,power->sign_normalization,&sign));
294: PetscCall(VecAXPY(y,-1.0,Bx));
295: /* Keep tracking eigenvalue update. It would be useful when we want to monitor solver progress via snes monitor. */
296: eps->eigr[(eps->nconv < eps->nev)? eps->nconv:(eps->nconv-1)] = 1.0/(bx*sign);
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*
301: Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
302: */
303: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
304: {
305: EPS_POWER *power = (EPS_POWER*)eps->data;
306: Vec Bx;
308: PetscFunctionBegin;
309: PetscCall(VecCopy(x,y));
310: if (power->update) PetscCall(SNESSolve(power->snes,NULL,y));
311: else {
312: Bx = eps->work[2];
313: PetscCall(SNESSolve(power->snes,Bx,y));
314: }
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*
319: Use nonlinear inverse power to compute an initial guess
320: */
321: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
322: {
323: EPS powereps;
324: Mat A,B,P;
325: Vec v1,v2;
326: SNES snes;
327: DM dm,newdm;
328: PetscBool sign_normalization;
330: PetscFunctionBegin;
331: PetscCall(EPSCreate(PetscObjectComm((PetscObject)eps),&powereps));
332: PetscCall(EPSGetOperators(eps,&A,&B));
333: PetscCall(EPSSetType(powereps,EPSPOWER));
334: PetscCall(EPSSetOperators(powereps,A,B));
335: PetscCall(EPSSetTolerances(powereps,1e-6,4));
336: PetscCall(EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix));
337: PetscCall(EPSAppendOptionsPrefix(powereps,"init_"));
338: PetscCall(EPSSetProblemType(powereps,EPS_GNHEP));
339: PetscCall(EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE));
340: PetscCall(EPSPowerSetNonlinear(powereps,PETSC_TRUE));
341: PetscCall(STGetPreconditionerMat(eps->st,&P));
342: /* attach dm to initial solve */
343: PetscCall(EPSPowerGetSNES(eps,&snes));
344: PetscCall(SNESGetDM(snes,&dm));
345: /* use dmshell to temporarily store snes context */
346: PetscCall(DMCreate(PetscObjectComm((PetscObject)eps),&newdm));
347: PetscCall(DMSetType(newdm,DMSHELL));
348: PetscCall(DMSetUp(newdm));
349: PetscCall(DMCopyDMSNES(dm,newdm));
350: PetscCall(EPSPowerGetSNES(powereps,&snes));
351: PetscCall(SNESSetDM(snes,dm));
352: PetscCall(EPSPowerGetSignNormalization(eps,&sign_normalization));
353: PetscCall(EPSPowerSetSignNormalization(powereps,sign_normalization));
354: PetscCall(EPSSetFromOptions(powereps));
355: if (P) PetscCall(STSetPreconditionerMat(powereps->st,P));
356: PetscCall(EPSSolve(powereps));
357: PetscCall(BVGetColumn(eps->V,0,&v2));
358: PetscCall(BVGetColumn(powereps->V,0,&v1));
359: PetscCall(VecCopy(v1,v2));
360: PetscCall(BVRestoreColumn(powereps->V,0,&v1));
361: PetscCall(BVRestoreColumn(eps->V,0,&v2));
362: PetscCall(EPSDestroy(&powereps));
363: /* restore context back to the old nonlinear solver */
364: PetscCall(DMCopyDMSNES(newdm,dm));
365: PetscCall(DMDestroy(&newdm));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: static PetscErrorCode EPSSolve_Power(EPS eps)
370: {
371: EPS_POWER *power = (EPS_POWER*)eps->data;
372: PetscInt k,ld;
373: Vec v,y,e,Bx;
374: Mat A;
375: KSP ksp;
376: PetscReal relerr,norm,norm1,rt1,rt2,cs1;
377: PetscScalar theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
378: PetscBool breakdown;
379: KSPConvergedReason reason;
380: SNESConvergedReason snesreason;
382: PetscFunctionBegin;
383: e = eps->work[0];
384: y = eps->work[1];
385: if (power->nonlinear) Bx = eps->work[2];
386: else Bx = NULL;
388: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
389: if (power->nonlinear) {
390: PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps));
391: /* Compute an initial guess only when users do not provide one */
392: if (power->update && !eps->nini) PetscCall(EPSPowerComputeInitialGuess_Update(eps));
393: } else PetscCall(DSGetLeadingDimension(eps->ds,&ld));
394: if (!power->update) PetscCall(EPSGetStartVector(eps,0,NULL));
395: if (power->nonlinear) {
396: PetscCall(BVGetColumn(eps->V,0,&v));
397: if (eps->nini) {
398: /* We scale the initial vector back if the initial vector was provided by users */
399: PetscCall(VecScale(v,power->norm0));
400: }
401: /* Do a couple things:
402: * 1) Determine the first non-zero index for Bx and the proc that owns it. This will be
403: * used if performing normalization by the sign of the first nonzero element of Bx.
404: * 2) Scale Bx by its norm. For non-update power iterations, Bx (in code terms) is used
405: * as the RHS argument to SNESSolve. And recall that the formula for generalized
406: * inverse power iterations in this case is: (Ax)_n = (Bx)_{n-1}/|(Bx)_{n-1}| (in
407: * math terms)
408: */
409: PetscCall(EPSPowerUpdateFunctionB(eps,v,Bx));
410: PetscCall(BVRestoreColumn(eps->V,0,&v));
411: if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
412: else PetscCall(VecNorm(Bx,NORM_2,&norm));
413: PetscCall(FirstNonzeroIdx(Bx,&power->idx,&power->p));
414: PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,NULL));
415: }
417: PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
418: rho = sigma;
420: while (eps->reason == EPS_CONVERGED_ITERATING) {
421: eps->its++;
422: k = eps->nconv;
424: /* y = OP v */
425: PetscCall(BVGetColumn(eps->V,k,&v));
426: if (power->nonlinear) PetscCall(EPSPowerApply_SNES(eps,v,y));
427: else PetscCall(STApply(eps->st,v,y));
428: PetscCall(BVRestoreColumn(eps->V,k,&v));
430: /* purge previously converged eigenvectors */
431: if (PetscUnlikely(power->nonlinear)) {
432: /* We do not need to call this for Newton eigensolver since eigenvalue is
433: * updated in function evaluations.
434: */
435: if (!power->update) {
436: PetscCall(EPSPowerUpdateFunctionB(eps,y,Bx));
437: if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
438: else PetscCall(VecNorm(Bx,NORM_2,&norm));
439: PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,&sign));
440: }
441: } else {
442: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&T));
443: PetscCall(BVSetActiveColumns(eps->V,0,k));
444: PetscCall(BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL));
445: }
447: /* theta = (v,y)_B */
448: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
449: PetscCall(BVDotVec(eps->V,y,&theta));
450: if (!power->nonlinear) {
451: T[k+k*ld] = theta;
452: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&T));
453: }
455: /* Eigenvalue is already stored in function evaluations.
456: * Assign eigenvalue to theta to make the rest of the code consistent
457: */
458: if (power->update) theta = eps->eigr[eps->nconv];
459: else if (power->nonlinear) theta = 1.0/(norm*sign); /* Eigenvalue: 1/|Bx| */
461: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
463: /* approximate eigenvalue is the Rayleigh quotient */
464: eps->eigr[eps->nconv] = theta;
466: /**
467: * If the Newton method (update, SNES) is used, we do not compute "relerr"
468: * since SNES determines the convergence.
469: */
470: if (PetscUnlikely(power->update)) relerr = 0.;
471: else {
472: /* compute relative error as ||y-theta v||_2/|theta| */
473: PetscCall(VecCopy(y,e));
474: PetscCall(BVGetColumn(eps->V,k,&v));
475: PetscCall(VecAXPY(e,power->nonlinear?-1.0:-theta,v));
476: PetscCall(BVRestoreColumn(eps->V,k,&v));
477: PetscCall(VecNorm(e,NORM_2,&relerr));
478: if (PetscUnlikely(power->nonlinear)) relerr *= PetscAbsScalar(theta);
479: else relerr /= PetscAbsScalar(theta);
480: }
482: } else { /* RQI */
484: /* delta = ||y||_B */
485: delta = norm;
487: /* compute relative error */
488: if (rho == 0.0) relerr = PETSC_MAX_REAL;
489: else relerr = 1.0 / (norm*PetscAbsScalar(rho));
491: /* approximate eigenvalue is the shift */
492: eps->eigr[eps->nconv] = rho;
494: /* compute new shift */
495: if (relerr<eps->tol) {
496: rho = sigma; /* if converged, restore original shift */
497: PetscCall(STSetShift(eps->st,rho));
498: } else {
499: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
500: if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
501: /* beta1 is the norm of the residual associated with R(v) */
502: PetscCall(BVGetColumn(eps->V,k,&v));
503: PetscCall(VecAXPY(v,-PetscConj(theta)/(delta*delta),y));
504: PetscCall(BVRestoreColumn(eps->V,k,&v));
505: PetscCall(BVScaleColumn(eps->V,k,1.0/delta));
506: PetscCall(BVNormColumn(eps->V,k,NORM_2,&norm1));
507: beta1 = norm1;
509: /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
510: PetscCall(STGetMatrix(eps->st,0,&A));
511: PetscCall(BVGetColumn(eps->V,k,&v));
512: PetscCall(MatMult(A,v,e));
513: PetscCall(VecDot(v,e,&alpha2));
514: PetscCall(BVRestoreColumn(eps->V,k,&v));
515: alpha2 = alpha2 / (beta1 * beta1);
517: /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
518: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
519: PetscCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
520: PetscCall(PetscFPTrapPop());
521: if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
522: else rho = rt2;
523: }
524: /* update operator according to new shift */
525: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
526: PetscCall(STSetShift(eps->st,rho));
527: PetscCall(KSPGetConvergedReason(ksp,&reason));
528: if (reason) {
529: PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
530: rho *= 1+10*PETSC_MACHINE_EPSILON;
531: PetscCall(STSetShift(eps->st,rho));
532: PetscCall(KSPGetConvergedReason(ksp,&reason));
533: PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
534: }
535: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
536: }
537: }
538: eps->errest[eps->nconv] = relerr;
540: /* normalize */
541: if (!power->nonlinear) PetscCall(Normalize(y,norm,power->idx,power->p,power->sign_normalization,NULL));
542: PetscCall(BVInsertVec(eps->V,k,y));
544: if (PetscUnlikely(power->update)) {
545: PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
546: /* For Newton eigensolver, we are ready to return once SNES converged. */
547: if (snesreason>0) eps->nconv = 1;
548: } else if (PetscUnlikely(relerr<eps->tol)) { /* accept eigenpair */
549: eps->nconv = eps->nconv + 1;
550: if (eps->nconv<eps->nev) {
551: PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
552: if (breakdown) {
553: eps->reason = EPS_DIVERGED_BREAKDOWN;
554: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
555: break;
556: }
557: }
558: }
559: /* For Newton eigensolver, monitor will be called from SNES monitor */
560: if (!power->update) PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
561: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
563: /**
564: * When a customized stopping test is used, and reason can be set to be converged (EPS_CONVERGED_USER).
565: * In this case, we need to increase eps->nconv to "1" so users can retrieve the solution.
566: */
567: if (PetscUnlikely(power->nonlinear && eps->reason>0)) eps->nconv = 1;
568: }
570: if (power->nonlinear) PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",NULL));
571: else {
572: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
573: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
574: }
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode EPSSolve_TS_Power(EPS eps)
579: {
580: EPS_POWER *power = (EPS_POWER*)eps->data;
581: PetscInt k,ld;
582: Vec v,w,y,e,z;
583: KSP ksp;
584: PetscReal relerr=1.0,relerrl,delta;
585: PetscScalar theta,rho,alpha,sigma;
586: PetscBool breakdown,breakdownl;
587: KSPConvergedReason reason;
589: PetscFunctionBegin;
590: e = eps->work[0];
591: v = eps->work[1];
592: w = eps->work[2];
594: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
595: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
596: PetscCall(EPSGetStartVector(eps,0,NULL));
597: PetscCall(EPSGetLeftStartVector(eps,0,NULL));
598: PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL));
599: PetscCall(BVCopyVec(eps->V,0,v));
600: PetscCall(BVCopyVec(eps->W,0,w));
601: PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
602: rho = sigma;
604: while (eps->reason == EPS_CONVERGED_ITERATING) {
605: eps->its++;
606: k = eps->nconv;
608: /* y = OP v, z = OP' w */
609: PetscCall(BVGetColumn(eps->V,k,&y));
610: PetscCall(STApply(eps->st,v,y));
611: PetscCall(BVRestoreColumn(eps->V,k,&y));
612: PetscCall(BVGetColumn(eps->W,k,&z));
613: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
614: PetscCall(BVRestoreColumn(eps->W,k,&z));
616: /* purge previously converged eigenvectors */
617: PetscCall(BVBiorthogonalizeColumn(eps->V,eps->W,k));
619: /* theta = (w,y)_B */
620: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
621: PetscCall(BVDotVec(eps->V,w,&theta));
622: theta = PetscConj(theta);
624: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
626: /* approximate eigenvalue is the Rayleigh quotient */
627: eps->eigr[eps->nconv] = theta;
629: /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
630: PetscCall(BVCopyVec(eps->V,k,e));
631: PetscCall(VecAXPY(e,-theta,v));
632: PetscCall(VecNorm(e,NORM_2,&relerr));
633: PetscCall(BVCopyVec(eps->W,k,e));
634: PetscCall(VecAXPY(e,-PetscConj(theta),w));
635: PetscCall(VecNorm(e,NORM_2,&relerrl));
636: relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
637: }
639: /* normalize */
640: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
641: PetscCall(BVGetColumn(eps->W,k,&z));
642: PetscCall(BVDotVec(eps->V,z,&alpha));
643: PetscCall(BVRestoreColumn(eps->W,k,&z));
644: delta = PetscSqrtReal(PetscAbsScalar(alpha));
645: PetscCheck(delta!=0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in two-sided Power/RQI");
646: PetscCall(BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta)));
647: PetscCall(BVScaleColumn(eps->W,k,1.0/delta));
648: PetscCall(BVCopyVec(eps->V,k,v));
649: PetscCall(BVCopyVec(eps->W,k,w));
651: if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */
653: /* compute relative error */
654: if (rho == 0.0) relerr = PETSC_MAX_REAL;
655: else relerr = 1.0 / PetscAbsScalar(delta*rho);
657: /* approximate eigenvalue is the shift */
658: eps->eigr[eps->nconv] = rho;
660: /* compute new shift */
661: if (relerr<eps->tol) {
662: rho = sigma; /* if converged, restore original shift */
663: PetscCall(STSetShift(eps->st,rho));
664: } else {
665: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
666: /* update operator according to new shift */
667: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
668: PetscCall(STSetShift(eps->st,rho));
669: PetscCall(KSPGetConvergedReason(ksp,&reason));
670: if (reason) {
671: PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
672: rho *= 1+10*PETSC_MACHINE_EPSILON;
673: PetscCall(STSetShift(eps->st,rho));
674: PetscCall(KSPGetConvergedReason(ksp,&reason));
675: PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
676: }
677: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
678: }
679: }
680: eps->errest[eps->nconv] = relerr;
682: /* if relerr<tol, accept eigenpair */
683: if (relerr<eps->tol) {
684: eps->nconv = eps->nconv + 1;
685: if (eps->nconv<eps->nev) {
686: PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
687: PetscCall(EPSGetLeftStartVector(eps,eps->nconv,&breakdownl));
688: if (breakdown || breakdownl) {
689: eps->reason = EPS_DIVERGED_BREAKDOWN;
690: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
691: break;
692: }
693: PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL));
694: }
695: }
696: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
697: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
698: }
700: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
701: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: static PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
706: {
707: EPS_POWER *power = (EPS_POWER*)eps->data;
708: SNESConvergedReason snesreason;
710: PetscFunctionBegin;
711: if (PetscUnlikely(power->update)) {
712: PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
713: if (snesreason < 0) {
714: *reason = EPS_DIVERGED_BREAKDOWN;
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
717: }
718: PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx));
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: static PetscErrorCode EPSBackTransform_Power(EPS eps)
723: {
724: EPS_POWER *power = (EPS_POWER*)eps->data;
726: PetscFunctionBegin;
727: if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
728: else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) PetscCall(EPSBackTransform_Default(eps));
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: static PetscErrorCode EPSSetFromOptions_Power(EPS eps,PetscOptionItems PetscOptionsObject)
733: {
734: EPS_POWER *power = (EPS_POWER*)eps->data;
735: PetscBool flg,val;
736: EPSPowerShiftType shift;
738: PetscFunctionBegin;
739: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Power Options");
741: PetscCall(PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg));
742: if (flg) PetscCall(EPSPowerSetShiftType(eps,shift));
744: PetscCall(PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg));
745: if (flg) PetscCall(EPSPowerSetNonlinear(eps,val));
747: PetscCall(PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg));
748: if (flg) PetscCall(EPSPowerSetUpdate(eps,val));
750: PetscCall(PetscOptionsBool("-eps_power_sign_normalization","Normalize Bx with sign of first nonzero entry","EPSPowerSetSignNormalization",power->sign_normalization,&power->sign_normalization,&flg));
752: PetscOptionsHeadEnd();
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
757: {
758: EPS_POWER *power = (EPS_POWER*)eps->data;
760: PetscFunctionBegin;
761: switch (shift) {
762: case EPS_POWER_SHIFT_CONSTANT:
763: case EPS_POWER_SHIFT_RAYLEIGH:
764: case EPS_POWER_SHIFT_WILKINSON:
765: if (power->shift_type != shift) {
766: power->shift_type = shift;
767: eps->state = EPS_STATE_INITIAL;
768: }
769: break;
770: default:
771: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
772: }
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /*@
777: EPSPowerSetShiftType - Sets the type of shifts used during the power
778: iteration. This can be used to emulate the Rayleigh Quotient Iteration
779: (RQI) method.
781: Logically Collective
783: Input Parameters:
784: + eps - the linear eigensolver context
785: - shift - the type of shift, see `EPSPowerShiftType` for possible values
787: Options Database Key:
788: . -eps_power_shift_type \<shift\> - sets the shift type, either `constant`, `rayleigh` or `wilkinson`
790: Notes:
791: By default, shifts are constant (`EPS_POWER_SHIFT_CONSTANT`) and the iteration
792: is the simple power method (or inverse iteration if a shift-and-invert
793: transformation is being used).
795: A variable shift can be specified (`EPS_POWER_SHIFT_RAYLEIGH` or
796: `EPS_POWER_SHIFT_WILKINSON`). In this case, the iteration behaves rather like
797: a cubic converging method such as RQI.
799: Details of the three variants can be found in {cite:p}`Her05`.
801: Level: advanced
803: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetShiftType()`, `STSetShift()`, `EPSPowerShiftType`
804: @*/
805: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
806: {
807: PetscFunctionBegin;
810: PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
815: {
816: EPS_POWER *power = (EPS_POWER*)eps->data;
818: PetscFunctionBegin;
819: *shift = power->shift_type;
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: /*@
824: EPSPowerGetShiftType - Gets the type of shifts used during the power
825: iteration.
827: Not Collective
829: Input Parameter:
830: . eps - the linear eigensolver context
832: Output Parameter:
833: . shift - the type of shift
835: Level: advanced
837: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetShiftType()`, `EPSPowerShiftType`
838: @*/
839: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
840: {
841: PetscFunctionBegin;
843: PetscAssertPointer(shift,2);
844: PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
849: {
850: EPS_POWER *power = (EPS_POWER*)eps->data;
852: PetscFunctionBegin;
853: if (power->nonlinear != nonlinear) {
854: power->nonlinear = nonlinear;
855: eps->useds = PetscNot(nonlinear);
856: eps->ops->setupsort = nonlinear? NULL: EPSSetUpSort_Default;
857: eps->state = EPS_STATE_INITIAL;
858: }
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: /*@
863: EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.
865: Logically Collective
867: Input Parameters:
868: + eps - the linear eigensolver context
869: - nonlinear - whether the problem is nonlinear or not
871: Options Database Key:
872: . -eps_power_nonlinear - sets the nonlinear flag
874: Notes:
875: If this flag is set, the solver assumes that the problem is nonlinear,
876: that is, the operators that define the eigenproblem are not constant
877: matrices, but depend on the eigenvector, $A(x)x=\lambda B(x)x$. This is
878: different from the case of nonlinearity with respect to the eigenvalue
879: (use the `NEP` solver class for this kind of problems).
881: The way in which nonlinear operators are specified is very similar to
882: the case of PETSc's `SNES` solver. The difference is that the callback
883: functions are provided via composed functions `formFunction` and
884: `formJacobian` in each of the matrix objects passed as arguments of
885: `EPSSetOperators()`. The application context required for these functions
886: can be attached via a composed `PetscContainer`.
888: Level: advanced
890: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetNonlinear()`, `EPSSetOperators()`
891: @*/
892: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
893: {
894: PetscFunctionBegin;
897: PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
902: {
903: EPS_POWER *power = (EPS_POWER*)eps->data;
905: PetscFunctionBegin;
906: *nonlinear = power->nonlinear;
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.
913: Not Collective
915: Input Parameter:
916: . eps - the linear eigensolver context
918: Output Parameter:
919: . nonlinear - the nonlinear flag
921: Level: advanced
923: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetUpdate()`, `EPSPowerSetNonlinear()`
924: @*/
925: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
926: {
927: PetscFunctionBegin;
929: PetscAssertPointer(nonlinear,2);
930: PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
935: {
936: EPS_POWER *power = (EPS_POWER*)eps->data;
938: PetscFunctionBegin;
939: PetscCheck(power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
940: power->update = update;
941: eps->state = EPS_STATE_INITIAL;
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
947: for nonlinear problems. This potentially has a better convergence.
949: Logically Collective
951: Input Parameters:
952: + eps - the linear eigensolver context
953: - update - whether the residual is updated monolithically or not
955: Options Database Key:
956: . -eps_power_update - sets the update flag
958: Note:
959: This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.
961: Level: advanced
963: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerGetUpdate()`, `EPSPowerGetNonlinear()`, `EPSSetOperators()`
964: @*/
965: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
966: {
967: PetscFunctionBegin;
970: PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
975: {
976: EPS_POWER *power = (EPS_POWER*)eps->data;
978: PetscFunctionBegin;
979: *update = power->update;
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
985: for nonlinear problems.
987: Not Collective
989: Input Parameter:
990: . eps - the linear eigensolver context
992: Output Parameter:
993: . update - the update flag
995: Level: advanced
997: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetUpdate()`, `EPSPowerSetNonlinear()`
998: @*/
999: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
1000: {
1001: PetscFunctionBegin;
1003: PetscAssertPointer(update,2);
1004: PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: static PetscErrorCode EPSPowerSetSignNormalization_Power(EPS eps,PetscBool sign_normalization)
1009: {
1010: EPS_POWER *power = (EPS_POWER*)eps->data;
1012: PetscFunctionBegin;
1013: power->sign_normalization = sign_normalization;
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: /*@
1018: EPSPowerSetSignNormalization - Sets a flag to indicate whether the $Bx$ vector
1019: should be normalized by the sign of the first non-zero element in the vector.
1020: E.g., if this is true, the post-normalization value of the first non-zero element
1021: in the vector is guaranteed to be positive.
1023: Logically Collective
1025: Input Parameters:
1026: + eps - the linear eigensolver context
1027: - sign_normalization - whether $Bx$ should be multiplied by the sign of the first non-zero
1028: element when performing normalization steps
1030: Options Database Key:
1031: . -eps_power_sign_normalization - sets the sign normalization flag
1033: Note:
1034: This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.
1036: Level: advanced
1038: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetNonlinear()`, `EPSPowerGetSignNormalization()`
1039: @*/
1040: PetscErrorCode EPSPowerSetSignNormalization(EPS eps,PetscBool sign_normalization)
1041: {
1042: PetscFunctionBegin;
1045: PetscTryMethod(eps,"EPSPowerSetSignNormalization_C",(EPS,PetscBool),(eps,sign_normalization));
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: static PetscErrorCode EPSPowerGetSignNormalization_Power(EPS eps,PetscBool *sign_normalization)
1050: {
1051: EPS_POWER *power = (EPS_POWER*)eps->data;
1053: PetscFunctionBegin;
1054: *sign_normalization = power->sign_normalization;
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: /*@
1059: EPSPowerGetSignNormalization - Returns a flag indicating whether the $Bx$ vector
1060: is normalized by the sign of the first non-zero element in the vector. E.g.,
1061: if this is true, the post-normalization value of the first non-zero element in
1062: the vector is guaranteed to be positive.
1064: Not Collective
1066: Input Parameter:
1067: . eps - the linear eigensolver context
1069: Output Parameter:
1070: . sign_normalization - the sign normalization flag
1072: Level: advanced
1074: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetSignNormalization()`
1075: @*/
1076: PetscErrorCode EPSPowerGetSignNormalization(EPS eps,PetscBool *sign_normalization)
1077: {
1078: PetscFunctionBegin;
1080: PetscAssertPointer(sign_normalization,2);
1081: PetscUseMethod(eps,"EPSPowerGetSignNormalization_C",(EPS,PetscBool*),(eps,sign_normalization));
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
1086: {
1087: EPS_POWER *power = (EPS_POWER*)eps->data;
1089: PetscFunctionBegin;
1090: PetscCall(PetscObjectReference((PetscObject)snes));
1091: PetscCall(SNESDestroy(&power->snes));
1092: power->snes = snes;
1093: eps->state = EPS_STATE_INITIAL;
1094: PetscFunctionReturn(PETSC_SUCCESS);
1095: }
1097: /*@
1098: EPSPowerSetSNES - Associate a nonlinear solver object (`SNES`) to the
1099: eigenvalue solver (to be used in nonlinear inverse iteration).
1101: Collective
1103: Input Parameters:
1104: + eps - the linear eigensolver context
1105: - snes - the nonlinear solver object
1107: Note:
1108: This flag is relevant only in nonlinear problems, see `EPSPowerSetNonlinear()`.
1110: Level: advanced
1112: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetNonlinear()`, `EPSPowerGetSNES()`
1113: @*/
1114: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1115: {
1116: PetscFunctionBegin;
1119: PetscCheckSameComm(eps,1,snes,2);
1120: PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1125: {
1126: EPS_POWER *power = (EPS_POWER*)eps->data;
1128: PetscFunctionBegin;
1129: if (!power->snes) {
1130: PetscCall(SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes));
1131: PetscCall(PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1));
1132: PetscCall(SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix));
1133: PetscCall(SNESAppendOptionsPrefix(power->snes,"eps_power_"));
1134: PetscCall(PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options));
1135: }
1136: *snes = power->snes;
1137: PetscFunctionReturn(PETSC_SUCCESS);
1138: }
1140: /*@
1141: EPSPowerGetSNES - Retrieve the nonlinear solver object (`SNES`) associated
1142: with the eigenvalue solver.
1144: Not Collective
1146: Input Parameter:
1147: . eps - the linear eigensolver context
1149: Output Parameter:
1150: . snes - the nonlinear solver object
1152: Level: advanced
1154: .seealso: [](ch:eps), `EPSPOWER`, `EPSPowerSetSNES()`
1155: @*/
1156: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1157: {
1158: PetscFunctionBegin;
1160: PetscAssertPointer(snes,2);
1161: PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: static PetscErrorCode EPSReset_Power(EPS eps)
1166: {
1167: EPS_POWER *power = (EPS_POWER*)eps->data;
1169: PetscFunctionBegin;
1170: if (power->snes) PetscCall(SNESReset(power->snes));
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: static PetscErrorCode EPSDestroy_Power(EPS eps)
1175: {
1176: EPS_POWER *power = (EPS_POWER*)eps->data;
1178: PetscFunctionBegin;
1179: if (power->nonlinear) PetscCall(SNESDestroy(&power->snes));
1180: PetscCall(PetscFree(eps->data));
1181: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL));
1182: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL));
1183: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL));
1184: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL));
1185: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL));
1186: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL));
1187: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",NULL));
1188: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",NULL));
1189: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL));
1190: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL));
1191: PetscFunctionReturn(PETSC_SUCCESS);
1192: }
1194: static PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1195: {
1196: EPS_POWER *power = (EPS_POWER*)eps->data;
1197: PetscBool isascii;
1199: PetscFunctionBegin;
1200: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1201: if (isascii) {
1202: if (power->nonlinear) {
1203: if (power->sign_normalization) PetscCall(PetscViewerASCIIPrintf(viewer," normalizing Bx by the sign of the first nonzero element\n"));
1204: else PetscCall(PetscViewerASCIIPrintf(viewer," not normalizing Bx by the sign of the first nonzero element\n"));
1205: PetscCall(PetscViewerASCIIPrintf(viewer," using nonlinear inverse iteration\n"));
1206: if (power->update) PetscCall(PetscViewerASCIIPrintf(viewer," updating the residual monolithically\n"));
1207: if (!power->snes) PetscCall(EPSPowerGetSNES(eps,&power->snes));
1208: PetscCall(PetscViewerASCIIPushTab(viewer));
1209: PetscCall(SNESView(power->snes,viewer));
1210: PetscCall(PetscViewerASCIIPopTab(viewer));
1211: } else PetscCall(PetscViewerASCIIPrintf(viewer," %s shifts\n",EPSPowerShiftTypes[power->shift_type]));
1212: }
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: static PetscErrorCode EPSComputeVectors_Power(EPS eps)
1217: {
1218: EPS_POWER *power = (EPS_POWER*)eps->data;
1220: PetscFunctionBegin;
1221: if (eps->twosided) {
1222: PetscCall(EPSComputeVectors_Twosided(eps));
1223: PetscCall(BVNormalize(eps->V,NULL));
1224: PetscCall(BVNormalize(eps->W,NULL));
1225: } else if (!power->nonlinear) PetscCall(EPSComputeVectors_Schur(eps));
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: static PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1230: {
1231: EPS_POWER *power = (EPS_POWER*)eps->data;
1232: KSP ksp;
1233: PC pc;
1235: PetscFunctionBegin;
1236: if (power->nonlinear) {
1237: eps->categ=EPS_CATEGORY_PRECOND;
1238: PetscCall(STGetKSP(eps->st,&ksp));
1239: /* Set ST as STPRECOND so it can carry one preconditioning matrix
1240: * It is useful when A and B are shell matrices
1241: */
1242: PetscCall(STSetType(eps->st,STPRECOND));
1243: PetscCall(KSPGetPC(ksp,&pc));
1244: PetscCall(PCSetType(pc,PCNONE));
1245: }
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: /*MC
1250: EPSPOWER - EPSPOWER = "power" - The simple power iteration and inverse
1251: iteration.
1253: Notes:
1254: This solver is very basic and is not recommended in general, since it
1255: will not be competitive with respect to other solvers.
1257: The implemented method is the power iteration, or inverse iteration in
1258: case the selected spectral transformation is `STSINVERT`. In this latter
1259: case it is possible to use a dynamic shift, as in the RQI method, see
1260: `EPSPowerSetShiftType()`.
1262: The solver incorporates deflation so that several eigenpairs can be
1263: computed. There is also a two-sided implementation that also computes
1264: left eigenvectors.
1266: This solver can also be used for nonlinear inverse iteration on the problem
1267: $A(x)x=\lambda B(x)x$, where $A$ and $B$ are not constant matrices but
1268: depend on the eigenvector $x$. This mode is enabled with `EPSPowerSetNonlinear()`.
1269: Note that this is a nonlinear eigenvector problem, as opposed to problems
1270: addressed in `NEP` that are nonlinear with respect to the eigenvalue.
1272: Level: beginner
1274: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSSetTwoSided()`, `EPSPowerSetShiftType()`, `EPSPowerSetNonlinear()`
1275: M*/
1276: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1277: {
1278: EPS_POWER *ctx;
1280: PetscFunctionBegin;
1281: PetscCall(PetscNew(&ctx));
1282: eps->data = (void*)ctx;
1284: eps->useds = PETSC_TRUE;
1285: eps->categ = EPS_CATEGORY_OTHER;
1287: eps->ops->setup = EPSSetUp_Power;
1288: eps->ops->setupsort = EPSSetUpSort_Default;
1289: eps->ops->setfromoptions = EPSSetFromOptions_Power;
1290: eps->ops->reset = EPSReset_Power;
1291: eps->ops->destroy = EPSDestroy_Power;
1292: eps->ops->view = EPSView_Power;
1293: eps->ops->backtransform = EPSBackTransform_Power;
1294: eps->ops->computevectors = EPSComputeVectors_Power;
1295: eps->ops->setdefaultst = EPSSetDefaultST_Power;
1296: eps->stopping = EPSStopping_Power;
1297: ctx->sign_normalization = PETSC_TRUE;
1299: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power));
1300: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power));
1301: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power));
1302: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power));
1303: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power));
1304: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power));
1305: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",EPSPowerSetSignNormalization_Power));
1306: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",EPSPowerGetSignNormalization_Power));
1307: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power));
1308: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power));
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }