Actual source code: power.c
slepc-3.21.0 2024-03-30
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: if (eps->ncv!=PETSC_DEFAULT) {
87: PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
88: } else eps->ncv = eps->nev;
89: if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
90: if (eps->max_it==PETSC_DEFAULT) {
91: /* SNES will directly return the solution for us, and we need to do only one iteration */
92: if (power->nonlinear && power->update) eps->max_it = 1;
93: else eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
94: }
95: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
96: PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
97: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
98: PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
99: EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," (with variable shifts)");
100: PetscCall(STGetMatMode(eps->st,&mode));
101: PetscCheck(mode!=ST_MATMODE_INPLACE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
102: }
103: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE);
104: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
105: PetscCall(EPSAllocateSolution(eps,0));
106: PetscCall(EPS_SetInnerProduct(eps));
108: if (power->nonlinear) {
109: PetscCheck(eps->nev==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
110: PetscCall(EPSSetWorkVecs(eps,3));
111: PetscCheck(!power->update || eps->max_it==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"More than one iteration is not allowed for Newton eigensolver (SNES)");
113: /* Set up SNES solver */
114: /* If SNES was setup, we need to reset it so that it will be consistent with the current EPS */
115: if (power->snes) PetscCall(SNESReset(power->snes));
116: else PetscCall(EPSPowerGetSNES(eps,&power->snes));
118: /* For nonlinear solver (Newton), we should scale the initial vector back.
119: The initial vector will be scaled in EPSSetUp. */
120: if (eps->IS) PetscCall(VecNorm(eps->IS[0],NORM_2,&power->norm0));
122: PetscCall(EPSGetOperators(eps,&A,&B));
124: /* form function */
125: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA));
126: PetscCheck(formFunctionA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
127: PetscCall(PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container));
128: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
129: else ctx = NULL;
130: PetscCall(MatCreateVecs(A,&res,NULL));
131: power->formFunctionA = formFunctionA;
132: power->formFunctionActx = ctx;
133: if (power->update) {
134: PetscCall(SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx));
135: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB));
136: PetscCall(SNESMonitorSet(power->snes,SNESMonitor_PowerUpdate,NULL,NULL));
137: }
138: else PetscCall(SNESSetFunction(power->snes,res,formFunctionA,ctx));
139: PetscCall(VecDestroy(&res));
141: /* form Jacobian */
142: PetscCall(PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA));
143: PetscCheck(formJacobianA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
144: PetscCall(PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container));
145: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
146: else ctx = NULL;
147: /* This allows users to compute a different preconditioning matrix than A.
148: * It is useful when A and B are shell matrices.
149: */
150: PetscCall(STGetPreconditionerMat(eps->st,&P));
151: /* If users do not provide a matrix, we simply use A */
152: PetscCall(SNESSetJacobian(power->snes,A,P? P:A,formJacobianA,ctx));
153: PetscCall(SNESSetFromOptions(power->snes));
154: PetscCall(SNESSetUp(power->snes));
156: ctx = NULL;
157: formNorm = NULL;
158: if (B) {
159: PetscCall(PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB));
160: PetscCall(PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container));
161: if (power->formFunctionB && container) PetscCall(PetscContainerGetPointer(container,&power->formFunctionBctx));
162: else power->formFunctionBctx = NULL;
164: /* form norm */
165: PetscCall(PetscObjectQueryFunction((PetscObject)B,"formNorm",&formNorm));
166: if (formNorm) {
167: PetscCall(PetscObjectQuery((PetscObject)B,"formNormCtx",(PetscObject*)&container));
168: if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
169: }
170: }
171: power->formNorm = formNorm;
172: power->formNormCtx = ctx;
173: } else {
174: if (eps->twosided) PetscCall(EPSSetWorkVecs(eps,3));
175: else PetscCall(EPSSetWorkVecs(eps,2));
176: PetscCall(DSSetType(eps->ds,DSNHEP));
177: PetscCall(DSAllocate(eps->ds,eps->nev));
178: }
179: /* dispatch solve method */
180: if (eps->twosided) {
181: PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
182: PetscCheck(power->shift_type!=EPS_POWER_SHIFT_WILKINSON,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
183: eps->ops->solve = EPSSolve_TS_Power;
184: } else eps->ops->solve = EPSSolve_Power;
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*
189: Find the index of the first nonzero entry of x, and the process that owns it.
190: */
191: static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscMPIInt *p)
192: {
193: PetscInt i,first,last,N;
194: PetscLayout map;
195: const PetscScalar *xx;
197: PetscFunctionBegin;
198: PetscCall(VecGetSize(x,&N));
199: PetscCall(VecGetOwnershipRange(x,&first,&last));
200: PetscCall(VecGetArrayRead(x,&xx));
201: for (i=first;i<last;i++) {
202: if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
203: }
204: PetscCall(VecRestoreArrayRead(x,&xx));
205: if (i==last) i=N;
206: PetscCall(MPIU_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x)));
207: PetscCheck(*idx!=N,PetscObjectComm((PetscObject)x),PETSC_ERR_PLIB,"Zero vector found");
208: PetscCall(VecGetLayout(x,&map));
209: PetscCall(PetscLayoutFindOwner(map,*idx,p));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*
214: Normalize a vector x with respect to a given norm as well as, optionally, the
215: sign of the first nonzero entry (assumed to be idx in process p).
216: */
217: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscMPIInt p,PetscBool sign_normalization,PetscScalar *sign)
218: {
219: PetscScalar alpha=1.0;
220: PetscInt first,last;
221: PetscReal absv;
222: const PetscScalar *xx;
224: PetscFunctionBegin;
225: if (sign_normalization) {
226: PetscCall(VecGetOwnershipRange(x,&first,&last));
227: if (idx>=first && idx<last) {
228: PetscCall(VecGetArrayRead(x,&xx));
229: absv = PetscAbsScalar(xx[idx-first]);
230: if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
231: PetscCall(VecRestoreArrayRead(x,&xx));
232: }
233: PetscCallMPI(MPI_Bcast(&alpha,1,MPIU_SCALAR,p,PetscObjectComm((PetscObject)x)));
234: }
235: if (sign) *sign = alpha;
236: alpha *= norm;
237: PetscCall(VecScale(x,1.0/alpha));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
242: {
243: EPS_POWER *power = (EPS_POWER*)eps->data;
244: Mat B;
246: PetscFunctionBegin;
247: PetscCall(STResetMatrixState(eps->st));
248: PetscCall(EPSGetOperators(eps,NULL,&B));
249: if (B) {
250: if (power->formFunctionB) PetscCall((*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx));
251: else PetscCall(MatMult(B,x,Bx));
252: } else PetscCall(VecCopy(x,Bx));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
257: {
258: EPS_POWER *power = (EPS_POWER*)eps->data;
259: Mat A;
261: PetscFunctionBegin;
262: PetscCall(STResetMatrixState(eps->st));
263: PetscCall(EPSGetOperators(eps,&A,NULL));
264: PetscCheck(A,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
265: if (power->formFunctionA) PetscCall((*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx));
266: else PetscCall(MatMult(A,x,Ax));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
271: {
272: EPS eps;
273: PetscReal bx;
274: Vec Bx;
275: PetscScalar sign;
276: EPS_POWER *power;
278: PetscFunctionBegin;
279: PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
280: PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
281: PetscCall(STResetMatrixState(eps->st));
282: power = (EPS_POWER*)eps->data;
283: Bx = eps->work[2];
284: if (power->formFunctionAB) PetscCall((*power->formFunctionAB)(snes,x,y,Bx,ctx));
285: else {
286: PetscCall(EPSPowerUpdateFunctionA(eps,x,y));
287: PetscCall(EPSPowerUpdateFunctionB(eps,x,Bx));
288: }
289: if (power->formNorm) PetscCall((*power->formNorm)(snes,Bx,&bx,power->formNormCtx));
290: else PetscCall(VecNorm(Bx,NORM_2,&bx));
291: PetscCall(Normalize(Bx,bx,power->idx,power->p,power->sign_normalization,&sign));
292: PetscCall(VecAXPY(y,-1.0,Bx));
293: /* Keep tracking eigenvalue update. It would be useful when we want to monitor solver progress via snes monitor. */
294: eps->eigr[(eps->nconv < eps->nev)? eps->nconv:(eps->nconv-1)] = 1.0/(bx*sign);
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*
299: Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
300: */
301: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
302: {
303: EPS_POWER *power = (EPS_POWER*)eps->data;
304: Vec Bx;
306: PetscFunctionBegin;
307: PetscCall(VecCopy(x,y));
308: if (power->update) PetscCall(SNESSolve(power->snes,NULL,y));
309: else {
310: Bx = eps->work[2];
311: PetscCall(SNESSolve(power->snes,Bx,y));
312: }
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*
317: Use nonlinear inverse power to compute an initial guess
318: */
319: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
320: {
321: EPS powereps;
322: Mat A,B,P;
323: Vec v1,v2;
324: SNES snes;
325: DM dm,newdm;
326: PetscBool sign_normalization;
328: PetscFunctionBegin;
329: PetscCall(EPSCreate(PetscObjectComm((PetscObject)eps),&powereps));
330: PetscCall(EPSGetOperators(eps,&A,&B));
331: PetscCall(EPSSetType(powereps,EPSPOWER));
332: PetscCall(EPSSetOperators(powereps,A,B));
333: PetscCall(EPSSetTolerances(powereps,1e-6,4));
334: PetscCall(EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix));
335: PetscCall(EPSAppendOptionsPrefix(powereps,"init_"));
336: PetscCall(EPSSetProblemType(powereps,EPS_GNHEP));
337: PetscCall(EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE));
338: PetscCall(EPSPowerSetNonlinear(powereps,PETSC_TRUE));
339: PetscCall(STGetPreconditionerMat(eps->st,&P));
340: /* attach dm to initial solve */
341: PetscCall(EPSPowerGetSNES(eps,&snes));
342: PetscCall(SNESGetDM(snes,&dm));
343: /* use dmshell to temporarily store snes context */
344: PetscCall(DMCreate(PetscObjectComm((PetscObject)eps),&newdm));
345: PetscCall(DMSetType(newdm,DMSHELL));
346: PetscCall(DMSetUp(newdm));
347: PetscCall(DMCopyDMSNES(dm,newdm));
348: PetscCall(EPSPowerGetSNES(powereps,&snes));
349: PetscCall(SNESSetDM(snes,dm));
350: PetscCall(EPSPowerGetSignNormalization(eps,&sign_normalization));
351: PetscCall(EPSPowerSetSignNormalization(powereps,sign_normalization));
352: PetscCall(EPSSetFromOptions(powereps));
353: if (P) PetscCall(STSetPreconditionerMat(powereps->st,P));
354: PetscCall(EPSSolve(powereps));
355: PetscCall(BVGetColumn(eps->V,0,&v2));
356: PetscCall(BVGetColumn(powereps->V,0,&v1));
357: PetscCall(VecCopy(v1,v2));
358: PetscCall(BVRestoreColumn(powereps->V,0,&v1));
359: PetscCall(BVRestoreColumn(eps->V,0,&v2));
360: PetscCall(EPSDestroy(&powereps));
361: /* restore context back to the old nonlinear solver */
362: PetscCall(DMCopyDMSNES(newdm,dm));
363: PetscCall(DMDestroy(&newdm));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: static PetscErrorCode EPSSolve_Power(EPS eps)
368: {
369: EPS_POWER *power = (EPS_POWER*)eps->data;
370: PetscInt k,ld;
371: Vec v,y,e,Bx;
372: Mat A;
373: KSP ksp;
374: PetscReal relerr,norm,norm1,rt1,rt2,cs1;
375: PetscScalar theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
376: PetscBool breakdown;
377: KSPConvergedReason reason;
378: SNESConvergedReason snesreason;
380: PetscFunctionBegin;
381: e = eps->work[0];
382: y = eps->work[1];
383: if (power->nonlinear) Bx = eps->work[2];
384: else Bx = NULL;
386: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
387: if (power->nonlinear) {
388: PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps));
389: /* Compute an initial guess only when users do not provide one */
390: if (power->update && !eps->nini) PetscCall(EPSPowerComputeInitialGuess_Update(eps));
391: } else PetscCall(DSGetLeadingDimension(eps->ds,&ld));
392: if (!power->update) PetscCall(EPSGetStartVector(eps,0,NULL));
393: if (power->nonlinear) {
394: PetscCall(BVGetColumn(eps->V,0,&v));
395: if (eps->nini) {
396: /* We scale the initial vector back if the initial vector was provided by users */
397: PetscCall(VecScale(v,power->norm0));
398: }
399: /* Do a couple things:
400: * 1) Determine the first non-zero index for Bx and the proc that owns it. This will be
401: * used if performing normalization by the sign of the first nonzero element of Bx.
402: * 2) Scale Bx by its norm. For non-update power iterations, Bx (in code terms) is used
403: * as the RHS argument to SNESSolve. And recall that the formula for generalized
404: * inverse power iterations in this case is: (Ax)_n = (Bx)_{n-1}/|(Bx)_{n-1}| (in
405: * math terms)
406: */
407: PetscCall(EPSPowerUpdateFunctionB(eps,v,Bx));
408: PetscCall(BVRestoreColumn(eps->V,0,&v));
409: if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
410: else PetscCall(VecNorm(Bx,NORM_2,&norm));
411: PetscCall(FirstNonzeroIdx(Bx,&power->idx,&power->p));
412: PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,NULL));
413: }
415: PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
416: rho = sigma;
418: while (eps->reason == EPS_CONVERGED_ITERATING) {
419: eps->its++;
420: k = eps->nconv;
422: /* y = OP v */
423: PetscCall(BVGetColumn(eps->V,k,&v));
424: if (power->nonlinear) PetscCall(EPSPowerApply_SNES(eps,v,y));
425: else PetscCall(STApply(eps->st,v,y));
426: PetscCall(BVRestoreColumn(eps->V,k,&v));
428: /* purge previously converged eigenvectors */
429: if (PetscUnlikely(power->nonlinear)) {
430: /* We do not need to call this for Newton eigensolver since eigenvalue is
431: * updated in function evaluations.
432: */
433: if (!power->update) {
434: PetscCall(EPSPowerUpdateFunctionB(eps,y,Bx));
435: if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
436: else PetscCall(VecNorm(Bx,NORM_2,&norm));
437: PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,&sign));
438: }
439: } else {
440: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&T));
441: PetscCall(BVSetActiveColumns(eps->V,0,k));
442: PetscCall(BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL));
443: }
445: /* theta = (v,y)_B */
446: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
447: PetscCall(BVDotVec(eps->V,y,&theta));
448: if (!power->nonlinear) {
449: T[k+k*ld] = theta;
450: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&T));
451: }
453: /* Eigenvalue is already stored in function evaluations.
454: * Assign eigenvalue to theta to make the rest of the code consistent
455: */
456: if (power->update) theta = eps->eigr[eps->nconv];
457: else if (power->nonlinear) theta = 1.0/(norm*sign); /* Eigenvalue: 1/|Bx| */
459: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
461: /* approximate eigenvalue is the Rayleigh quotient */
462: eps->eigr[eps->nconv] = theta;
464: /**
465: * If the Newton method (update, SNES) is used, we do not compute "relerr"
466: * since SNES determines the convergence.
467: */
468: if (PetscUnlikely(power->update)) relerr = 0.;
469: else {
470: /* compute relative error as ||y-theta v||_2/|theta| */
471: PetscCall(VecCopy(y,e));
472: PetscCall(BVGetColumn(eps->V,k,&v));
473: PetscCall(VecAXPY(e,power->nonlinear?-1.0:-theta,v));
474: PetscCall(BVRestoreColumn(eps->V,k,&v));
475: PetscCall(VecNorm(e,NORM_2,&relerr));
476: if (PetscUnlikely(power->nonlinear)) relerr *= PetscAbsScalar(theta);
477: else relerr /= PetscAbsScalar(theta);
478: }
480: } else { /* RQI */
482: /* delta = ||y||_B */
483: delta = norm;
485: /* compute relative error */
486: if (rho == 0.0) relerr = PETSC_MAX_REAL;
487: else relerr = 1.0 / (norm*PetscAbsScalar(rho));
489: /* approximate eigenvalue is the shift */
490: eps->eigr[eps->nconv] = rho;
492: /* compute new shift */
493: if (relerr<eps->tol) {
494: rho = sigma; /* if converged, restore original shift */
495: PetscCall(STSetShift(eps->st,rho));
496: } else {
497: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
498: if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
499: /* beta1 is the norm of the residual associated with R(v) */
500: PetscCall(BVGetColumn(eps->V,k,&v));
501: PetscCall(VecAXPY(v,-PetscConj(theta)/(delta*delta),y));
502: PetscCall(BVRestoreColumn(eps->V,k,&v));
503: PetscCall(BVScaleColumn(eps->V,k,1.0/delta));
504: PetscCall(BVNormColumn(eps->V,k,NORM_2,&norm1));
505: beta1 = norm1;
507: /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
508: PetscCall(STGetMatrix(eps->st,0,&A));
509: PetscCall(BVGetColumn(eps->V,k,&v));
510: PetscCall(MatMult(A,v,e));
511: PetscCall(VecDot(v,e,&alpha2));
512: PetscCall(BVRestoreColumn(eps->V,k,&v));
513: alpha2 = alpha2 / (beta1 * beta1);
515: /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
516: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
517: PetscCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
518: PetscCall(PetscFPTrapPop());
519: if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
520: else rho = rt2;
521: }
522: /* update operator according to new shift */
523: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
524: PetscCall(STSetShift(eps->st,rho));
525: PetscCall(KSPGetConvergedReason(ksp,&reason));
526: if (reason) {
527: PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
528: rho *= 1+10*PETSC_MACHINE_EPSILON;
529: PetscCall(STSetShift(eps->st,rho));
530: PetscCall(KSPGetConvergedReason(ksp,&reason));
531: PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
532: }
533: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
534: }
535: }
536: eps->errest[eps->nconv] = relerr;
538: /* normalize */
539: if (!power->nonlinear) PetscCall(Normalize(y,norm,power->idx,power->p,power->sign_normalization,NULL));
540: PetscCall(BVInsertVec(eps->V,k,y));
542: if (PetscUnlikely(power->update)) {
543: PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
544: /* For Newton eigensolver, we are ready to return once SNES converged. */
545: if (snesreason>0) eps->nconv = 1;
546: } else if (PetscUnlikely(relerr<eps->tol)) { /* accept eigenpair */
547: eps->nconv = eps->nconv + 1;
548: if (eps->nconv<eps->nev) {
549: PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
550: if (breakdown) {
551: eps->reason = EPS_DIVERGED_BREAKDOWN;
552: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
553: break;
554: }
555: }
556: }
557: /* For Newton eigensolver, monitor will be called from SNES monitor */
558: if (!power->update) PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
559: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
561: /**
562: * When a customized stopping test is used, and reason can be set to be converged (EPS_CONVERGED_USER).
563: * In this case, we need to increase eps->nconv to "1" so users can retrieve the solution.
564: */
565: if (PetscUnlikely(power->nonlinear && eps->reason>0)) eps->nconv = 1;
566: }
568: if (power->nonlinear) PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",NULL));
569: else {
570: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
571: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
572: }
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: static PetscErrorCode EPSSolve_TS_Power(EPS eps)
577: {
578: EPS_POWER *power = (EPS_POWER*)eps->data;
579: PetscInt k,ld;
580: Vec v,w,y,e,z;
581: KSP ksp;
582: PetscReal relerr=1.0,relerrl,delta;
583: PetscScalar theta,rho,alpha,sigma;
584: PetscBool breakdown,breakdownl;
585: KSPConvergedReason reason;
587: PetscFunctionBegin;
588: e = eps->work[0];
589: v = eps->work[1];
590: w = eps->work[2];
592: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
593: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
594: PetscCall(EPSGetStartVector(eps,0,NULL));
595: PetscCall(EPSGetLeftStartVector(eps,0,NULL));
596: PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL));
597: PetscCall(BVCopyVec(eps->V,0,v));
598: PetscCall(BVCopyVec(eps->W,0,w));
599: PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
600: rho = sigma;
602: while (eps->reason == EPS_CONVERGED_ITERATING) {
603: eps->its++;
604: k = eps->nconv;
606: /* y = OP v, z = OP' w */
607: PetscCall(BVGetColumn(eps->V,k,&y));
608: PetscCall(STApply(eps->st,v,y));
609: PetscCall(BVRestoreColumn(eps->V,k,&y));
610: PetscCall(BVGetColumn(eps->W,k,&z));
611: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
612: PetscCall(BVRestoreColumn(eps->W,k,&z));
614: /* purge previously converged eigenvectors */
615: PetscCall(BVBiorthogonalizeColumn(eps->V,eps->W,k));
617: /* theta = (w,y)_B */
618: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
619: PetscCall(BVDotVec(eps->V,w,&theta));
620: theta = PetscConj(theta);
622: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
624: /* approximate eigenvalue is the Rayleigh quotient */
625: eps->eigr[eps->nconv] = theta;
627: /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
628: PetscCall(BVCopyVec(eps->V,k,e));
629: PetscCall(VecAXPY(e,-theta,v));
630: PetscCall(VecNorm(e,NORM_2,&relerr));
631: PetscCall(BVCopyVec(eps->W,k,e));
632: PetscCall(VecAXPY(e,-PetscConj(theta),w));
633: PetscCall(VecNorm(e,NORM_2,&relerrl));
634: relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
635: }
637: /* normalize */
638: PetscCall(BVSetActiveColumns(eps->V,k,k+1));
639: PetscCall(BVGetColumn(eps->W,k,&z));
640: PetscCall(BVDotVec(eps->V,z,&alpha));
641: PetscCall(BVRestoreColumn(eps->W,k,&z));
642: delta = PetscSqrtReal(PetscAbsScalar(alpha));
643: PetscCheck(delta!=0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in two-sided Power/RQI");
644: PetscCall(BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta)));
645: PetscCall(BVScaleColumn(eps->W,k,1.0/delta));
646: PetscCall(BVCopyVec(eps->V,k,v));
647: PetscCall(BVCopyVec(eps->W,k,w));
649: if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */
651: /* compute relative error */
652: if (rho == 0.0) relerr = PETSC_MAX_REAL;
653: else relerr = 1.0 / PetscAbsScalar(delta*rho);
655: /* approximate eigenvalue is the shift */
656: eps->eigr[eps->nconv] = rho;
658: /* compute new shift */
659: if (relerr<eps->tol) {
660: rho = sigma; /* if converged, restore original shift */
661: PetscCall(STSetShift(eps->st,rho));
662: } else {
663: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
664: /* update operator according to new shift */
665: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
666: PetscCall(STSetShift(eps->st,rho));
667: PetscCall(KSPGetConvergedReason(ksp,&reason));
668: if (reason) {
669: PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
670: rho *= 1+10*PETSC_MACHINE_EPSILON;
671: PetscCall(STSetShift(eps->st,rho));
672: PetscCall(KSPGetConvergedReason(ksp,&reason));
673: PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
674: }
675: PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
676: }
677: }
678: eps->errest[eps->nconv] = relerr;
680: /* if relerr<tol, accept eigenpair */
681: if (relerr<eps->tol) {
682: eps->nconv = eps->nconv + 1;
683: if (eps->nconv<eps->nev) {
684: PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
685: PetscCall(EPSGetLeftStartVector(eps,eps->nconv,&breakdownl));
686: if (breakdown || breakdownl) {
687: eps->reason = EPS_DIVERGED_BREAKDOWN;
688: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
689: break;
690: }
691: PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL));
692: }
693: }
694: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
695: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
696: }
698: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
699: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
704: {
705: EPS_POWER *power = (EPS_POWER*)eps->data;
706: SNESConvergedReason snesreason;
708: PetscFunctionBegin;
709: if (PetscUnlikely(power->update)) {
710: PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
711: if (snesreason < 0) {
712: *reason = EPS_DIVERGED_BREAKDOWN;
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
715: }
716: PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode EPSBackTransform_Power(EPS eps)
721: {
722: EPS_POWER *power = (EPS_POWER*)eps->data;
724: PetscFunctionBegin;
725: if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
726: else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) PetscCall(EPSBackTransform_Default(eps));
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: static PetscErrorCode EPSSetFromOptions_Power(EPS eps,PetscOptionItems *PetscOptionsObject)
731: {
732: EPS_POWER *power = (EPS_POWER*)eps->data;
733: PetscBool flg,val;
734: EPSPowerShiftType shift;
736: PetscFunctionBegin;
737: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Power Options");
739: PetscCall(PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg));
740: if (flg) PetscCall(EPSPowerSetShiftType(eps,shift));
742: PetscCall(PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg));
743: if (flg) PetscCall(EPSPowerSetNonlinear(eps,val));
745: PetscCall(PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg));
746: if (flg) PetscCall(EPSPowerSetUpdate(eps,val));
748: PetscCall(PetscOptionsBool("-eps_power_sign_normalization","Normalize Bx with sign of first nonzero entry","EPSPowerSetSignNormalization",power->sign_normalization,&power->sign_normalization,&flg));
750: PetscOptionsHeadEnd();
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
755: {
756: EPS_POWER *power = (EPS_POWER*)eps->data;
758: PetscFunctionBegin;
759: switch (shift) {
760: case EPS_POWER_SHIFT_CONSTANT:
761: case EPS_POWER_SHIFT_RAYLEIGH:
762: case EPS_POWER_SHIFT_WILKINSON:
763: if (power->shift_type != shift) {
764: power->shift_type = shift;
765: eps->state = EPS_STATE_INITIAL;
766: }
767: break;
768: default:
769: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
770: }
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: EPSPowerSetShiftType - Sets the type of shifts used during the power
776: iteration. This can be used to emulate the Rayleigh Quotient Iteration
777: (RQI) method.
779: Logically Collective
781: Input Parameters:
782: + eps - the eigenproblem solver context
783: - shift - the type of shift
785: Options Database Key:
786: . -eps_power_shift_type - Sets the shift type (either 'constant' or
787: 'rayleigh' or 'wilkinson')
789: Notes:
790: By default, shifts are constant (EPS_POWER_SHIFT_CONSTANT) and the iteration
791: is the simple power method (or inverse iteration if a shift-and-invert
792: transformation is being used).
794: A variable shift can be specified (EPS_POWER_SHIFT_RAYLEIGH or
795: EPS_POWER_SHIFT_WILKINSON). In this case, the iteration behaves rather like
796: a cubic converging method such as RQI.
798: Level: advanced
800: .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
801: @*/
802: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
803: {
804: PetscFunctionBegin;
807: PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
812: {
813: EPS_POWER *power = (EPS_POWER*)eps->data;
815: PetscFunctionBegin;
816: *shift = power->shift_type;
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: /*@
821: EPSPowerGetShiftType - Gets the type of shifts used during the power
822: iteration.
824: Not Collective
826: Input Parameter:
827: . eps - the eigenproblem solver context
829: Output Parameter:
830: . shift - the type of shift
832: Level: advanced
834: .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
835: @*/
836: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
837: {
838: PetscFunctionBegin;
840: PetscAssertPointer(shift,2);
841: PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
845: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
846: {
847: EPS_POWER *power = (EPS_POWER*)eps->data;
849: PetscFunctionBegin;
850: if (power->nonlinear != nonlinear) {
851: power->nonlinear = nonlinear;
852: eps->useds = PetscNot(nonlinear);
853: eps->ops->setupsort = nonlinear? NULL: EPSSetUpSort_Default;
854: eps->state = EPS_STATE_INITIAL;
855: }
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: /*@
860: EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.
862: Logically Collective
864: Input Parameters:
865: + eps - the eigenproblem solver context
866: - nonlinear - whether the problem is nonlinear or not
868: Options Database Key:
869: . -eps_power_nonlinear - Sets the nonlinear flag
871: Notes:
872: If this flag is set, the solver assumes that the problem is nonlinear,
873: that is, the operators that define the eigenproblem are not constant
874: matrices, but depend on the eigenvector, A(x)*x=lambda*B(x)*x. This is
875: different from the case of nonlinearity with respect to the eigenvalue
876: (use the NEP solver class for this kind of problems).
878: The way in which nonlinear operators are specified is very similar to
879: the case of PETSc's SNES solver. The difference is that the callback
880: functions are provided via composed functions "formFunction" and
881: "formJacobian" in each of the matrix objects passed as arguments of
882: EPSSetOperators(). The application context required for these functions
883: can be attached via a composed PetscContainer.
885: Level: advanced
887: .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
888: @*/
889: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
890: {
891: PetscFunctionBegin;
894: PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
899: {
900: EPS_POWER *power = (EPS_POWER*)eps->data;
902: PetscFunctionBegin;
903: *nonlinear = power->nonlinear;
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: /*@
908: EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.
910: Not Collective
912: Input Parameter:
913: . eps - the eigenproblem solver context
915: Output Parameter:
916: . nonlinear - the nonlinear flag
918: Level: advanced
920: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
921: @*/
922: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
923: {
924: PetscFunctionBegin;
926: PetscAssertPointer(nonlinear,2);
927: PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
932: {
933: EPS_POWER *power = (EPS_POWER*)eps->data;
935: PetscFunctionBegin;
936: PetscCheck(power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
937: power->update = update;
938: eps->state = EPS_STATE_INITIAL;
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: /*@
943: EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
944: for nonlinear problems. This potentially has a better convergence.
946: Logically Collective
948: Input Parameters:
949: + eps - the eigenproblem solver context
950: - update - whether the residual is updated monolithically or not
952: Options Database Key:
953: . -eps_power_update - Sets the update flag
955: Level: advanced
957: .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
958: @*/
959: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
960: {
961: PetscFunctionBegin;
964: PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
969: {
970: EPS_POWER *power = (EPS_POWER*)eps->data;
972: PetscFunctionBegin;
973: *update = power->update;
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: /*@
978: EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
979: for nonlinear problems.
981: Not Collective
983: Input Parameter:
984: . eps - the eigenproblem solver context
986: Output Parameter:
987: . update - the update flag
989: Level: advanced
991: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
992: @*/
993: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
994: {
995: PetscFunctionBegin;
997: PetscAssertPointer(update,2);
998: PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: static PetscErrorCode EPSPowerSetSignNormalization_Power(EPS eps,PetscBool sign_normalization)
1003: {
1004: EPS_POWER *power = (EPS_POWER*)eps->data;
1006: PetscFunctionBegin;
1007: power->sign_normalization = sign_normalization;
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@
1012: EPSPowerSetSignNormalization - Sets a flag to indicate whether the Bx vector
1013: should be normalized by the sign of the first non-zero element in the vector.
1014: E.g., if this is true, the post-normalization value of the first non-zero element
1015: in the vector is guaranteed to be positive.
1017: Logically Collective
1019: Input Parameters:
1020: + eps - the eigenproblem solver context
1021: - sign_normalization - whether Bx should be multiplied by the sign of the first non-zero
1022: element when performing normalization steps
1024: Options Database Key:
1025: . -eps_power_sign_normalization - Sets the sign normalization flag
1027: Level: advanced
1029: .seealso: EPSPowerGetSignNormalization()
1030: @*/
1031: PetscErrorCode EPSPowerSetSignNormalization(EPS eps,PetscBool sign_normalization)
1032: {
1033: PetscFunctionBegin;
1036: PetscTryMethod(eps,"EPSPowerSetSignNormalization_C",(EPS,PetscBool),(eps,sign_normalization));
1037: PetscFunctionReturn(PETSC_SUCCESS);
1038: }
1040: static PetscErrorCode EPSPowerGetSignNormalization_Power(EPS eps,PetscBool *sign_normalization)
1041: {
1042: EPS_POWER *power = (EPS_POWER*)eps->data;
1044: PetscFunctionBegin;
1045: *sign_normalization = power->sign_normalization;
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: /*@
1050: EPSPowerGetSignNormalization - Returns a flag indicating whether the Bx vector
1051: is normalized by the sign of the first non-zero element in the vector. E.g.,
1052: if this is true, the post-normalization value of the first non-zero element in
1053: the vector is guaranteed to be positive.
1055: Not Collective
1057: Input Parameter:
1058: . eps - the eigenproblem solver context
1060: Output Parameter:
1061: . sign_normalization - the sign normalization flag
1063: Level: advanced
1065: .seealso: EPSPowerSetSignNormalization()
1066: @*/
1067: PetscErrorCode EPSPowerGetSignNormalization(EPS eps,PetscBool *sign_normalization)
1068: {
1069: PetscFunctionBegin;
1071: PetscAssertPointer(sign_normalization,2);
1072: PetscUseMethod(eps,"EPSPowerGetSignNormalization_C",(EPS,PetscBool*),(eps,sign_normalization));
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }
1076: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
1077: {
1078: EPS_POWER *power = (EPS_POWER*)eps->data;
1080: PetscFunctionBegin;
1081: PetscCall(PetscObjectReference((PetscObject)snes));
1082: PetscCall(SNESDestroy(&power->snes));
1083: power->snes = snes;
1084: eps->state = EPS_STATE_INITIAL;
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: /*@
1089: EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
1090: eigenvalue solver (to be used in nonlinear inverse iteration).
1092: Collective
1094: Input Parameters:
1095: + eps - the eigenvalue solver
1096: - snes - the nonlinear solver object
1098: Level: advanced
1100: .seealso: EPSPowerGetSNES()
1101: @*/
1102: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1103: {
1104: PetscFunctionBegin;
1107: PetscCheckSameComm(eps,1,snes,2);
1108: PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1113: {
1114: EPS_POWER *power = (EPS_POWER*)eps->data;
1116: PetscFunctionBegin;
1117: if (!power->snes) {
1118: PetscCall(SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes));
1119: PetscCall(PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1));
1120: PetscCall(SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix));
1121: PetscCall(SNESAppendOptionsPrefix(power->snes,"eps_power_"));
1122: PetscCall(PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options));
1123: }
1124: *snes = power->snes;
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: /*@
1129: EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1130: with the eigenvalue solver.
1132: Not Collective
1134: Input Parameter:
1135: . eps - the eigenvalue solver
1137: Output Parameter:
1138: . snes - the nonlinear solver object
1140: Level: advanced
1142: .seealso: EPSPowerSetSNES()
1143: @*/
1144: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1145: {
1146: PetscFunctionBegin;
1148: PetscAssertPointer(snes,2);
1149: PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: static PetscErrorCode EPSReset_Power(EPS eps)
1154: {
1155: EPS_POWER *power = (EPS_POWER*)eps->data;
1157: PetscFunctionBegin;
1158: if (power->snes) PetscCall(SNESReset(power->snes));
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: static PetscErrorCode EPSDestroy_Power(EPS eps)
1163: {
1164: EPS_POWER *power = (EPS_POWER*)eps->data;
1166: PetscFunctionBegin;
1167: if (power->nonlinear) PetscCall(SNESDestroy(&power->snes));
1168: PetscCall(PetscFree(eps->data));
1169: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL));
1170: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL));
1171: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL));
1172: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL));
1173: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL));
1174: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL));
1175: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",NULL));
1176: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",NULL));
1177: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL));
1178: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL));
1179: PetscFunctionReturn(PETSC_SUCCESS);
1180: }
1182: static PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1183: {
1184: EPS_POWER *power = (EPS_POWER*)eps->data;
1185: PetscBool isascii;
1187: PetscFunctionBegin;
1188: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1189: if (isascii) {
1190: if (power->nonlinear) {
1191: if (power->sign_normalization) PetscCall(PetscViewerASCIIPrintf(viewer," normalizing Bx by the sign of the first nonzero element\n"));
1192: else PetscCall(PetscViewerASCIIPrintf(viewer," not normalizing Bx by the sign of the first nonzero element\n"));
1193: PetscCall(PetscViewerASCIIPrintf(viewer," using nonlinear inverse iteration\n"));
1194: if (power->update) PetscCall(PetscViewerASCIIPrintf(viewer," updating the residual monolithically\n"));
1195: if (!power->snes) PetscCall(EPSPowerGetSNES(eps,&power->snes));
1196: PetscCall(PetscViewerASCIIPushTab(viewer));
1197: PetscCall(SNESView(power->snes,viewer));
1198: PetscCall(PetscViewerASCIIPopTab(viewer));
1199: } else PetscCall(PetscViewerASCIIPrintf(viewer," %s shifts\n",EPSPowerShiftTypes[power->shift_type]));
1200: }
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: static PetscErrorCode EPSComputeVectors_Power(EPS eps)
1205: {
1206: EPS_POWER *power = (EPS_POWER*)eps->data;
1208: PetscFunctionBegin;
1209: if (eps->twosided) {
1210: PetscCall(EPSComputeVectors_Twosided(eps));
1211: PetscCall(BVNormalize(eps->V,NULL));
1212: PetscCall(BVNormalize(eps->W,NULL));
1213: } else if (!power->nonlinear) PetscCall(EPSComputeVectors_Schur(eps));
1214: PetscFunctionReturn(PETSC_SUCCESS);
1215: }
1217: static PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1218: {
1219: EPS_POWER *power = (EPS_POWER*)eps->data;
1220: KSP ksp;
1221: PC pc;
1223: PetscFunctionBegin;
1224: if (power->nonlinear) {
1225: eps->categ=EPS_CATEGORY_PRECOND;
1226: PetscCall(STGetKSP(eps->st,&ksp));
1227: /* Set ST as STPRECOND so it can carry one preconditioning matrix
1228: * It is useful when A and B are shell matrices
1229: */
1230: PetscCall(STSetType(eps->st,STPRECOND));
1231: PetscCall(KSPGetPC(ksp,&pc));
1232: PetscCall(PCSetType(pc,PCNONE));
1233: }
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1238: {
1239: EPS_POWER *ctx;
1241: PetscFunctionBegin;
1242: PetscCall(PetscNew(&ctx));
1243: eps->data = (void*)ctx;
1245: eps->useds = PETSC_TRUE;
1246: eps->categ = EPS_CATEGORY_OTHER;
1248: eps->ops->setup = EPSSetUp_Power;
1249: eps->ops->setupsort = EPSSetUpSort_Default;
1250: eps->ops->setfromoptions = EPSSetFromOptions_Power;
1251: eps->ops->reset = EPSReset_Power;
1252: eps->ops->destroy = EPSDestroy_Power;
1253: eps->ops->view = EPSView_Power;
1254: eps->ops->backtransform = EPSBackTransform_Power;
1255: eps->ops->computevectors = EPSComputeVectors_Power;
1256: eps->ops->setdefaultst = EPSSetDefaultST_Power;
1257: eps->stopping = EPSStopping_Power;
1258: ctx->sign_normalization = PETSC_TRUE;
1260: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power));
1261: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power));
1262: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power));
1263: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power));
1264: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power));
1265: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power));
1266: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",EPSPowerSetSignNormalization_Power));
1267: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",EPSPowerGetSignNormalization_Power));
1268: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power));
1269: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power));
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }