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