Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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"
12 :
13 : Method: Power Iteration
14 :
15 : Algorithm:
16 :
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.
23 :
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.
26 :
27 : References:
28 :
29 : [1] "Single Vector Iteration Methods in SLEPc", SLEPc Technical Report
30 : STR-2, available at https://slepc.upv.es.
31 : */
32 :
33 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
34 : #include <slepcblaslapack.h>
35 : /* petsc headers */
36 : #include <petscdm.h>
37 : #include <petscsnes.h>
38 :
39 : static PetscErrorCode EPSPowerFormFunction_Update(SNES,Vec,Vec,void*);
40 : static PetscErrorCode EPSSolve_Power(EPS);
41 : static PetscErrorCode EPSSolve_TS_Power(EPS);
42 :
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;
60 :
61 104 : static PetscErrorCode SNESMonitor_PowerUpdate(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
62 : {
63 104 : EPS eps;
64 :
65 104 : PetscFunctionBegin;
66 104 : PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
67 104 : PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
68 : /* Call EPS monitor on each SNES iteration */
69 104 : PetscCall(EPSMonitor(eps,its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
70 104 : PetscFunctionReturn(PETSC_SUCCESS);
71 : }
72 :
73 48 : static PetscErrorCode EPSSetUp_Power(EPS eps)
74 : {
75 48 : EPS_POWER *power = (EPS_POWER*)eps->data;
76 48 : STMatMode mode;
77 48 : Mat A,B,P;
78 48 : Vec res;
79 48 : PetscContainer container;
80 48 : PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
81 48 : PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
82 48 : PetscErrorCode (*formNorm)(SNES,Vec,PetscReal*,void*);
83 48 : void *ctx;
84 :
85 48 : PetscFunctionBegin;
86 48 : if (eps->ncv!=PETSC_DEFAULT) {
87 5 : PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
88 43 : } else eps->ncv = eps->nev;
89 48 : if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
90 48 : if (eps->max_it==PETSC_DEFAULT) {
91 : /* SNES will directly return the solution for us, and we need to do only one iteration */
92 28 : if (power->nonlinear && power->update) eps->max_it = 1;
93 16 : else eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
94 : }
95 48 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
96 48 : 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 48 : if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
98 4 : PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
99 4 : EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," (with variable shifts)");
100 4 : PetscCall(STGetMatMode(eps->st,&mode));
101 4 : PetscCheck(mode!=ST_MATMODE_INPLACE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
102 : }
103 48 : EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE);
104 48 : EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
105 48 : PetscCall(EPSAllocateSolution(eps,0));
106 48 : PetscCall(EPS_SetInnerProduct(eps));
107 :
108 48 : if (power->nonlinear) {
109 31 : PetscCheck(eps->nev==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
110 31 : PetscCall(EPSSetWorkVecs(eps,3));
111 31 : PetscCheck(!power->update || eps->max_it==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"More than one iteration is not allowed for Newton eigensolver (SNES)");
112 :
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 31 : if (power->snes) PetscCall(SNESReset(power->snes));
116 0 : else PetscCall(EPSPowerGetSNES(eps,&power->snes));
117 :
118 : /* For nonlinear solver (Newton), we should scale the initial vector back.
119 : The initial vector will be scaled in EPSSetUp. */
120 31 : if (eps->IS) PetscCall(VecNorm(eps->IS[0],NORM_2,&power->norm0));
121 :
122 31 : PetscCall(EPSGetOperators(eps,&A,&B));
123 :
124 : /* form function */
125 31 : PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA));
126 31 : PetscCheck(formFunctionA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
127 31 : PetscCall(PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container));
128 31 : if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
129 0 : else ctx = NULL;
130 31 : PetscCall(MatCreateVecs(A,&res,NULL));
131 31 : power->formFunctionA = formFunctionA;
132 31 : power->formFunctionActx = ctx;
133 31 : if (power->update) {
134 14 : PetscCall(SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx));
135 14 : PetscCall(PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB));
136 14 : PetscCall(SNESMonitorSet(power->snes,SNESMonitor_PowerUpdate,NULL,NULL));
137 : }
138 17 : else PetscCall(SNESSetFunction(power->snes,res,formFunctionA,ctx));
139 31 : PetscCall(VecDestroy(&res));
140 :
141 : /* form Jacobian */
142 31 : PetscCall(PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA));
143 31 : PetscCheck(formJacobianA,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
144 31 : PetscCall(PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container));
145 31 : if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
146 0 : 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 31 : PetscCall(STGetPreconditionerMat(eps->st,&P));
151 : /* If users do not provide a matrix, we simply use A */
152 31 : PetscCall(SNESSetJacobian(power->snes,A,P? P:A,formJacobianA,ctx));
153 31 : PetscCall(SNESSetFromOptions(power->snes));
154 31 : PetscCall(SNESSetUp(power->snes));
155 :
156 31 : ctx = NULL;
157 31 : formNorm = NULL;
158 31 : if (B) {
159 31 : PetscCall(PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB));
160 31 : PetscCall(PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container));
161 31 : if (power->formFunctionB && container) PetscCall(PetscContainerGetPointer(container,&power->formFunctionBctx));
162 0 : else power->formFunctionBctx = NULL;
163 :
164 : /* form norm */
165 31 : PetscCall(PetscObjectQueryFunction((PetscObject)B,"formNorm",&formNorm));
166 31 : if (formNorm) {
167 8 : PetscCall(PetscObjectQuery((PetscObject)B,"formNormCtx",(PetscObject*)&container));
168 8 : if (container) PetscCall(PetscContainerGetPointer(container,&ctx));
169 : }
170 : }
171 31 : power->formNorm = formNorm;
172 31 : power->formNormCtx = ctx;
173 : } else {
174 17 : if (eps->twosided) PetscCall(EPSSetWorkVecs(eps,3));
175 14 : else PetscCall(EPSSetWorkVecs(eps,2));
176 17 : PetscCall(DSSetType(eps->ds,DSNHEP));
177 17 : PetscCall(DSAllocate(eps->ds,eps->nev));
178 : }
179 : /* dispatch solve method */
180 48 : if (eps->twosided) {
181 3 : PetscCheck(!power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
182 3 : PetscCheck(power->shift_type!=EPS_POWER_SHIFT_WILKINSON,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
183 3 : eps->ops->solve = EPSSolve_TS_Power;
184 45 : } else eps->ops->solve = EPSSolve_Power;
185 48 : PetscFunctionReturn(PETSC_SUCCESS);
186 : }
187 :
188 : /*
189 : Find the index of the first nonzero entry of x, and the process that owns it.
190 : */
191 31 : static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscMPIInt *p)
192 : {
193 31 : PetscInt i,first,last,N;
194 31 : PetscLayout map;
195 31 : const PetscScalar *xx;
196 :
197 31 : PetscFunctionBegin;
198 31 : PetscCall(VecGetSize(x,&N));
199 31 : PetscCall(VecGetOwnershipRange(x,&first,&last));
200 31 : PetscCall(VecGetArrayRead(x,&xx));
201 341 : for (i=first;i<last;i++) {
202 341 : if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
203 : }
204 31 : PetscCall(VecRestoreArrayRead(x,&xx));
205 31 : if (i==last) i=N;
206 124 : PetscCall(MPIU_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x)));
207 31 : PetscCheck(*idx!=N,PetscObjectComm((PetscObject)x),PETSC_ERR_PLIB,"Zero vector found");
208 31 : PetscCall(VecGetLayout(x,&map));
209 31 : PetscCall(PetscLayoutFindOwner(map,*idx,p));
210 31 : PetscFunctionReturn(PETSC_SUCCESS);
211 : }
212 :
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 29308 : static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscMPIInt p,PetscBool sign_normalization,PetscScalar *sign)
218 : {
219 29308 : PetscScalar alpha=1.0;
220 29308 : PetscInt first,last;
221 29308 : PetscReal absv;
222 29308 : const PetscScalar *xx;
223 :
224 29308 : PetscFunctionBegin;
225 29308 : if (sign_normalization) {
226 29022 : PetscCall(VecGetOwnershipRange(x,&first,&last));
227 29022 : if (idx>=first && idx<last) {
228 29022 : PetscCall(VecGetArrayRead(x,&xx));
229 29022 : absv = PetscAbsScalar(xx[idx-first]);
230 29022 : if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
231 29022 : PetscCall(VecRestoreArrayRead(x,&xx));
232 : }
233 58044 : PetscCallMPI(MPI_Bcast(&alpha,1,MPIU_SCALAR,p,PetscObjectComm((PetscObject)x)));
234 : }
235 29308 : if (sign) *sign = alpha;
236 29308 : alpha *= norm;
237 29308 : PetscCall(VecScale(x,1.0/alpha));
238 29308 : PetscFunctionReturn(PETSC_SUCCESS);
239 : }
240 :
241 661 : static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
242 : {
243 661 : EPS_POWER *power = (EPS_POWER*)eps->data;
244 661 : Mat B;
245 :
246 661 : PetscFunctionBegin;
247 661 : PetscCall(STResetMatrixState(eps->st));
248 661 : PetscCall(EPSGetOperators(eps,NULL,&B));
249 661 : if (B) {
250 661 : if (power->formFunctionB) PetscCall((*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx));
251 0 : else PetscCall(MatMult(B,x,Bx));
252 0 : } else PetscCall(VecCopy(x,Bx));
253 661 : PetscFunctionReturn(PETSC_SUCCESS);
254 : }
255 :
256 445 : static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
257 : {
258 445 : EPS_POWER *power = (EPS_POWER*)eps->data;
259 445 : Mat A;
260 :
261 445 : PetscFunctionBegin;
262 445 : PetscCall(STResetMatrixState(eps->st));
263 445 : PetscCall(EPSGetOperators(eps,&A,NULL));
264 445 : PetscCheck(A,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
265 445 : if (power->formFunctionA) PetscCall((*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx));
266 0 : else PetscCall(MatMult(A,x,Ax));
267 445 : PetscFunctionReturn(PETSC_SUCCESS);
268 : }
269 :
270 682 : static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
271 : {
272 682 : EPS eps;
273 682 : PetscReal bx;
274 682 : Vec Bx;
275 682 : PetscScalar sign;
276 682 : EPS_POWER *power;
277 :
278 682 : PetscFunctionBegin;
279 682 : PetscCall(PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps));
280 682 : PetscCheck(eps,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
281 682 : PetscCall(STResetMatrixState(eps->st));
282 682 : power = (EPS_POWER*)eps->data;
283 682 : Bx = eps->work[2];
284 682 : if (power->formFunctionAB) PetscCall((*power->formFunctionAB)(snes,x,y,Bx,ctx));
285 : else {
286 445 : PetscCall(EPSPowerUpdateFunctionA(eps,x,y));
287 445 : PetscCall(EPSPowerUpdateFunctionB(eps,x,Bx));
288 : }
289 682 : if (power->formNorm) PetscCall((*power->formNorm)(snes,Bx,&bx,power->formNormCtx));
290 472 : else PetscCall(VecNorm(Bx,NORM_2,&bx));
291 682 : PetscCall(Normalize(Bx,bx,power->idx,power->p,power->sign_normalization,&sign));
292 682 : 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 682 : eps->eigr[(eps->nconv < eps->nev)? eps->nconv:(eps->nconv-1)] = 1.0/(bx*sign);
295 682 : PetscFunctionReturn(PETSC_SUCCESS);
296 : }
297 :
298 : /*
299 : Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
300 : */
301 199 : static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
302 : {
303 199 : EPS_POWER *power = (EPS_POWER*)eps->data;
304 199 : Vec Bx;
305 :
306 199 : PetscFunctionBegin;
307 199 : PetscCall(VecCopy(x,y));
308 199 : if (power->update) PetscCall(SNESSolve(power->snes,NULL,y));
309 : else {
310 185 : Bx = eps->work[2];
311 185 : PetscCall(SNESSolve(power->snes,Bx,y));
312 : }
313 199 : PetscFunctionReturn(PETSC_SUCCESS);
314 : }
315 :
316 : /*
317 : Use nonlinear inverse power to compute an initial guess
318 : */
319 12 : static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
320 : {
321 12 : EPS powereps;
322 12 : Mat A,B,P;
323 12 : Vec v1,v2;
324 12 : SNES snes;
325 12 : DM dm,newdm;
326 12 : PetscBool sign_normalization;
327 :
328 12 : PetscFunctionBegin;
329 12 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)eps),&powereps));
330 12 : PetscCall(EPSGetOperators(eps,&A,&B));
331 12 : PetscCall(EPSSetType(powereps,EPSPOWER));
332 12 : PetscCall(EPSSetOperators(powereps,A,B));
333 12 : PetscCall(EPSSetTolerances(powereps,1e-6,4));
334 12 : PetscCall(EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix));
335 12 : PetscCall(EPSAppendOptionsPrefix(powereps,"init_"));
336 12 : PetscCall(EPSSetProblemType(powereps,EPS_GNHEP));
337 12 : PetscCall(EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE));
338 12 : PetscCall(EPSPowerSetNonlinear(powereps,PETSC_TRUE));
339 12 : PetscCall(STGetPreconditionerMat(eps->st,&P));
340 : /* attach dm to initial solve */
341 12 : PetscCall(EPSPowerGetSNES(eps,&snes));
342 12 : PetscCall(SNESGetDM(snes,&dm));
343 : /* use dmshell to temporarily store snes context */
344 12 : PetscCall(DMCreate(PetscObjectComm((PetscObject)eps),&newdm));
345 12 : PetscCall(DMSetType(newdm,DMSHELL));
346 12 : PetscCall(DMSetUp(newdm));
347 12 : PetscCall(DMCopyDMSNES(dm,newdm));
348 12 : PetscCall(EPSPowerGetSNES(powereps,&snes));
349 12 : PetscCall(SNESSetDM(snes,dm));
350 12 : PetscCall(EPSPowerGetSignNormalization(eps,&sign_normalization));
351 12 : PetscCall(EPSPowerSetSignNormalization(powereps,sign_normalization));
352 12 : PetscCall(EPSSetFromOptions(powereps));
353 12 : if (P) PetscCall(STSetPreconditionerMat(powereps->st,P));
354 12 : PetscCall(EPSSolve(powereps));
355 12 : PetscCall(BVGetColumn(eps->V,0,&v2));
356 12 : PetscCall(BVGetColumn(powereps->V,0,&v1));
357 12 : PetscCall(VecCopy(v1,v2));
358 12 : PetscCall(BVRestoreColumn(powereps->V,0,&v1));
359 12 : PetscCall(BVRestoreColumn(eps->V,0,&v2));
360 12 : PetscCall(EPSDestroy(&powereps));
361 : /* restore context back to the old nonlinear solver */
362 12 : PetscCall(DMCopyDMSNES(newdm,dm));
363 12 : PetscCall(DMDestroy(&newdm));
364 12 : PetscFunctionReturn(PETSC_SUCCESS);
365 : }
366 :
367 45 : static PetscErrorCode EPSSolve_Power(EPS eps)
368 : {
369 45 : EPS_POWER *power = (EPS_POWER*)eps->data;
370 45 : PetscInt k,ld;
371 45 : Vec v,y,e,Bx;
372 45 : Mat A;
373 45 : KSP ksp;
374 45 : PetscReal relerr,norm,norm1,rt1,rt2,cs1;
375 45 : PetscScalar theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
376 45 : PetscBool breakdown;
377 45 : KSPConvergedReason reason;
378 45 : SNESConvergedReason snesreason;
379 :
380 45 : PetscFunctionBegin;
381 45 : e = eps->work[0];
382 45 : y = eps->work[1];
383 45 : if (power->nonlinear) Bx = eps->work[2];
384 : else Bx = NULL;
385 :
386 45 : if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
387 45 : if (power->nonlinear) {
388 31 : PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps));
389 : /* Compute an initial guess only when users do not provide one */
390 31 : if (power->update && !eps->nini) PetscCall(EPSPowerComputeInitialGuess_Update(eps));
391 14 : } else PetscCall(DSGetLeadingDimension(eps->ds,&ld));
392 45 : if (!power->update) PetscCall(EPSGetStartVector(eps,0,NULL));
393 45 : if (power->nonlinear) {
394 31 : PetscCall(BVGetColumn(eps->V,0,&v));
395 31 : if (eps->nini) {
396 : /* We scale the initial vector back if the initial vector was provided by users */
397 2 : 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 31 : PetscCall(EPSPowerUpdateFunctionB(eps,v,Bx));
408 31 : PetscCall(BVRestoreColumn(eps->V,0,&v));
409 31 : if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
410 23 : else PetscCall(VecNorm(Bx,NORM_2,&norm));
411 31 : PetscCall(FirstNonzeroIdx(Bx,&power->idx,&power->p));
412 31 : PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,NULL));
413 : }
414 :
415 45 : PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
416 45 : rho = sigma;
417 :
418 28654 : while (eps->reason == EPS_CONVERGED_ITERATING) {
419 28609 : eps->its++;
420 28609 : k = eps->nconv;
421 :
422 : /* y = OP v */
423 28609 : PetscCall(BVGetColumn(eps->V,k,&v));
424 28609 : if (power->nonlinear) PetscCall(EPSPowerApply_SNES(eps,v,y));
425 28410 : else PetscCall(STApply(eps->st,v,y));
426 28609 : PetscCall(BVRestoreColumn(eps->V,k,&v));
427 :
428 : /* purge previously converged eigenvectors */
429 28609 : 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 199 : if (!power->update) {
434 185 : PetscCall(EPSPowerUpdateFunctionB(eps,y,Bx));
435 185 : if (power->formNorm) PetscCall((*power->formNorm)(power->snes,Bx,&norm,power->formNormCtx));
436 115 : else PetscCall(VecNorm(Bx,NORM_2,&norm));
437 185 : PetscCall(Normalize(Bx,norm,power->idx,power->p,power->sign_normalization,&sign));
438 : }
439 : } else {
440 28410 : PetscCall(DSGetArray(eps->ds,DS_MAT_A,&T));
441 28410 : PetscCall(BVSetActiveColumns(eps->V,0,k));
442 28410 : PetscCall(BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL));
443 : }
444 :
445 : /* theta = (v,y)_B */
446 28609 : PetscCall(BVSetActiveColumns(eps->V,k,k+1));
447 28609 : PetscCall(BVDotVec(eps->V,y,&theta));
448 28609 : if (!power->nonlinear) {
449 28410 : T[k+k*ld] = theta;
450 28410 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&T));
451 : }
452 :
453 : /* Eigenvalue is already stored in function evaluations.
454 : * Assign eigenvalue to theta to make the rest of the code consistent
455 : */
456 28609 : if (power->update) theta = eps->eigr[eps->nconv];
457 28595 : else if (power->nonlinear) theta = 1.0/(norm*sign); /* Eigenvalue: 1/|Bx| */
458 :
459 28609 : if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
460 :
461 : /* approximate eigenvalue is the Rayleigh quotient */
462 28575 : eps->eigr[eps->nconv] = theta;
463 :
464 : /**
465 : * If the Newton method (update, SNES) is used, we do not compute "relerr"
466 : * since SNES determines the convergence.
467 : */
468 28575 : if (PetscUnlikely(power->update)) relerr = 0.;
469 : else {
470 : /* compute relative error as ||y-theta v||_2/|theta| */
471 28561 : PetscCall(VecCopy(y,e));
472 28561 : PetscCall(BVGetColumn(eps->V,k,&v));
473 28561 : PetscCall(VecAXPY(e,power->nonlinear?-1.0:-theta,v));
474 28561 : PetscCall(BVRestoreColumn(eps->V,k,&v));
475 28561 : PetscCall(VecNorm(e,NORM_2,&relerr));
476 28561 : if (PetscUnlikely(power->nonlinear)) relerr *= PetscAbsScalar(theta);
477 28376 : else relerr /= PetscAbsScalar(theta);
478 : }
479 :
480 : } else { /* RQI */
481 :
482 : /* delta = ||y||_B */
483 34 : delta = norm;
484 :
485 : /* compute relative error */
486 34 : if (rho == 0.0) relerr = PETSC_MAX_REAL;
487 34 : else relerr = 1.0 / (norm*PetscAbsScalar(rho));
488 :
489 : /* approximate eigenvalue is the shift */
490 34 : eps->eigr[eps->nconv] = rho;
491 :
492 : /* compute new shift */
493 34 : if (relerr<eps->tol) {
494 8 : rho = sigma; /* if converged, restore original shift */
495 8 : PetscCall(STSetShift(eps->st,rho));
496 : } else {
497 26 : rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
498 26 : if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
499 : /* beta1 is the norm of the residual associated with R(v) */
500 9 : PetscCall(BVGetColumn(eps->V,k,&v));
501 9 : PetscCall(VecAXPY(v,-PetscConj(theta)/(delta*delta),y));
502 9 : PetscCall(BVRestoreColumn(eps->V,k,&v));
503 9 : PetscCall(BVScaleColumn(eps->V,k,1.0/delta));
504 9 : PetscCall(BVNormColumn(eps->V,k,NORM_2,&norm1));
505 9 : beta1 = norm1;
506 :
507 : /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
508 9 : PetscCall(STGetMatrix(eps->st,0,&A));
509 9 : PetscCall(BVGetColumn(eps->V,k,&v));
510 9 : PetscCall(MatMult(A,v,e));
511 9 : PetscCall(VecDot(v,e,&alpha2));
512 9 : PetscCall(BVRestoreColumn(eps->V,k,&v));
513 9 : alpha2 = alpha2 / (beta1 * beta1);
514 :
515 : /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
516 9 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
517 9 : PetscCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
518 9 : PetscCall(PetscFPTrapPop());
519 9 : if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
520 3 : else rho = rt2;
521 : }
522 : /* update operator according to new shift */
523 26 : PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
524 26 : PetscCall(STSetShift(eps->st,rho));
525 26 : PetscCall(KSPGetConvergedReason(ksp,&reason));
526 26 : if (reason) {
527 0 : PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
528 0 : rho *= 1+10*PETSC_MACHINE_EPSILON;
529 0 : PetscCall(STSetShift(eps->st,rho));
530 0 : PetscCall(KSPGetConvergedReason(ksp,&reason));
531 0 : PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
532 : }
533 26 : PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
534 : }
535 : }
536 28609 : eps->errest[eps->nconv] = relerr;
537 :
538 : /* normalize */
539 28609 : if (!power->nonlinear) PetscCall(Normalize(y,norm,power->idx,power->p,power->sign_normalization,NULL));
540 28609 : PetscCall(BVInsertVec(eps->V,k,y));
541 :
542 28609 : if (PetscUnlikely(power->update)) {
543 14 : PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
544 : /* For Newton eigensolver, we are ready to return once SNES converged. */
545 14 : if (snesreason>0) eps->nconv = 1;
546 28595 : } else if (PetscUnlikely(relerr<eps->tol)) { /* accept eigenpair */
547 50 : eps->nconv = eps->nconv + 1;
548 50 : if (eps->nconv<eps->nev) {
549 31 : PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
550 31 : if (breakdown) {
551 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
552 0 : 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 28609 : if (!power->update) PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
559 28609 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
560 :
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 28609 : if (PetscUnlikely(power->nonlinear && eps->reason>0)) eps->nconv = 1;
566 : }
567 :
568 45 : if (power->nonlinear) PetscCall(PetscObjectCompose((PetscObject)power->snes,"eps",NULL));
569 : else {
570 14 : PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
571 14 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
572 : }
573 45 : PetscFunctionReturn(PETSC_SUCCESS);
574 : }
575 :
576 3 : static PetscErrorCode EPSSolve_TS_Power(EPS eps)
577 : {
578 3 : EPS_POWER *power = (EPS_POWER*)eps->data;
579 3 : PetscInt k,ld;
580 3 : Vec v,w,y,e,z;
581 3 : KSP ksp;
582 3 : PetscReal relerr=1.0,relerrl,delta;
583 3 : PetscScalar theta,rho,alpha,sigma;
584 3 : PetscBool breakdown,breakdownl;
585 3 : KSPConvergedReason reason;
586 :
587 3 : PetscFunctionBegin;
588 3 : e = eps->work[0];
589 3 : v = eps->work[1];
590 3 : w = eps->work[2];
591 :
592 3 : if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) PetscCall(STGetKSP(eps->st,&ksp));
593 3 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
594 3 : PetscCall(EPSGetStartVector(eps,0,NULL));
595 3 : PetscCall(EPSGetLeftStartVector(eps,0,NULL));
596 3 : PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL));
597 3 : PetscCall(BVCopyVec(eps->V,0,v));
598 3 : PetscCall(BVCopyVec(eps->W,0,w));
599 3 : PetscCall(STGetShift(eps->st,&sigma)); /* original shift */
600 3 : rho = sigma;
601 :
602 3 : while (eps->reason == EPS_CONVERGED_ITERATING) {
603 423 : eps->its++;
604 423 : k = eps->nconv;
605 :
606 : /* y = OP v, z = OP' w */
607 423 : PetscCall(BVGetColumn(eps->V,k,&y));
608 423 : PetscCall(STApply(eps->st,v,y));
609 423 : PetscCall(BVRestoreColumn(eps->V,k,&y));
610 423 : PetscCall(BVGetColumn(eps->W,k,&z));
611 423 : PetscCall(STApplyHermitianTranspose(eps->st,w,z));
612 423 : PetscCall(BVRestoreColumn(eps->W,k,&z));
613 :
614 : /* purge previously converged eigenvectors */
615 423 : PetscCall(BVBiorthogonalizeColumn(eps->V,eps->W,k));
616 :
617 : /* theta = (w,y)_B */
618 423 : PetscCall(BVSetActiveColumns(eps->V,k,k+1));
619 423 : PetscCall(BVDotVec(eps->V,w,&theta));
620 423 : theta = PetscConj(theta);
621 :
622 423 : if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
623 :
624 : /* approximate eigenvalue is the Rayleigh quotient */
625 411 : eps->eigr[eps->nconv] = theta;
626 :
627 : /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
628 411 : PetscCall(BVCopyVec(eps->V,k,e));
629 411 : PetscCall(VecAXPY(e,-theta,v));
630 411 : PetscCall(VecNorm(e,NORM_2,&relerr));
631 411 : PetscCall(BVCopyVec(eps->W,k,e));
632 411 : PetscCall(VecAXPY(e,-PetscConj(theta),w));
633 411 : PetscCall(VecNorm(e,NORM_2,&relerrl));
634 792 : relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
635 : }
636 :
637 : /* normalize */
638 423 : PetscCall(BVSetActiveColumns(eps->V,k,k+1));
639 423 : PetscCall(BVGetColumn(eps->W,k,&z));
640 423 : PetscCall(BVDotVec(eps->V,z,&alpha));
641 423 : PetscCall(BVRestoreColumn(eps->W,k,&z));
642 423 : delta = PetscSqrtReal(PetscAbsScalar(alpha));
643 423 : PetscCheck(delta!=0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in two-sided Power/RQI");
644 423 : PetscCall(BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta)));
645 423 : PetscCall(BVScaleColumn(eps->W,k,1.0/delta));
646 423 : PetscCall(BVCopyVec(eps->V,k,v));
647 423 : PetscCall(BVCopyVec(eps->W,k,w));
648 :
649 423 : if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */
650 :
651 : /* compute relative error */
652 12 : if (rho == 0.0) relerr = PETSC_MAX_REAL;
653 12 : else relerr = 1.0 / PetscAbsScalar(delta*rho);
654 :
655 : /* approximate eigenvalue is the shift */
656 12 : eps->eigr[eps->nconv] = rho;
657 :
658 : /* compute new shift */
659 12 : if (relerr<eps->tol) {
660 2 : rho = sigma; /* if converged, restore original shift */
661 2 : PetscCall(STSetShift(eps->st,rho));
662 : } else {
663 10 : rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
664 : /* update operator according to new shift */
665 10 : PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_FALSE));
666 10 : PetscCall(STSetShift(eps->st,rho));
667 10 : PetscCall(KSPGetConvergedReason(ksp,&reason));
668 10 : if (reason) {
669 0 : PetscCall(PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n"));
670 0 : rho *= 1+10*PETSC_MACHINE_EPSILON;
671 0 : PetscCall(STSetShift(eps->st,rho));
672 0 : PetscCall(KSPGetConvergedReason(ksp,&reason));
673 0 : PetscCheck(!reason,PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
674 : }
675 10 : PetscCall(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE));
676 : }
677 : }
678 423 : eps->errest[eps->nconv] = relerr;
679 :
680 : /* if relerr<tol, accept eigenpair */
681 423 : if (relerr<eps->tol) {
682 8 : eps->nconv = eps->nconv + 1;
683 8 : if (eps->nconv<eps->nev) {
684 5 : PetscCall(EPSGetStartVector(eps,eps->nconv,&breakdown));
685 5 : PetscCall(EPSGetLeftStartVector(eps,eps->nconv,&breakdownl));
686 5 : if (breakdown || breakdownl) {
687 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
688 0 : PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
689 : break;
690 : }
691 5 : PetscCall(BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL));
692 : }
693 : }
694 423 : PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev)));
695 426 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
696 : }
697 :
698 3 : PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
699 3 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
700 3 : PetscFunctionReturn(PETSC_SUCCESS);
701 : }
702 :
703 29032 : static PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
704 : {
705 29032 : EPS_POWER *power = (EPS_POWER*)eps->data;
706 29032 : SNESConvergedReason snesreason;
707 :
708 29032 : PetscFunctionBegin;
709 29032 : if (PetscUnlikely(power->update)) {
710 14 : PetscCall(SNESGetConvergedReason(power->snes,&snesreason));
711 14 : if (snesreason < 0) {
712 0 : *reason = EPS_DIVERGED_BREAKDOWN;
713 0 : PetscFunctionReturn(PETSC_SUCCESS);
714 : }
715 : }
716 29032 : PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx));
717 29032 : PetscFunctionReturn(PETSC_SUCCESS);
718 : }
719 :
720 17 : static PetscErrorCode EPSBackTransform_Power(EPS eps)
721 : {
722 17 : EPS_POWER *power = (EPS_POWER*)eps->data;
723 :
724 17 : PetscFunctionBegin;
725 17 : if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
726 17 : else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) PetscCall(EPSBackTransform_Default(eps));
727 17 : PetscFunctionReturn(PETSC_SUCCESS);
728 : }
729 :
730 44 : static PetscErrorCode EPSSetFromOptions_Power(EPS eps,PetscOptionItems *PetscOptionsObject)
731 : {
732 44 : EPS_POWER *power = (EPS_POWER*)eps->data;
733 44 : PetscBool flg,val;
734 44 : EPSPowerShiftType shift;
735 :
736 44 : PetscFunctionBegin;
737 44 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS Power Options");
738 :
739 44 : PetscCall(PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg));
740 44 : if (flg) PetscCall(EPSPowerSetShiftType(eps,shift));
741 :
742 44 : PetscCall(PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg));
743 44 : if (flg) PetscCall(EPSPowerSetNonlinear(eps,val));
744 :
745 44 : PetscCall(PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg));
746 44 : if (flg) PetscCall(EPSPowerSetUpdate(eps,val));
747 :
748 44 : PetscCall(PetscOptionsBool("-eps_power_sign_normalization","Normalize Bx with sign of first nonzero entry","EPSPowerSetSignNormalization",power->sign_normalization,&power->sign_normalization,&flg));
749 :
750 44 : PetscOptionsHeadEnd();
751 44 : PetscFunctionReturn(PETSC_SUCCESS);
752 : }
753 :
754 7 : static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
755 : {
756 7 : EPS_POWER *power = (EPS_POWER*)eps->data;
757 :
758 7 : PetscFunctionBegin;
759 7 : switch (shift) {
760 7 : case EPS_POWER_SHIFT_CONSTANT:
761 : case EPS_POWER_SHIFT_RAYLEIGH:
762 : case EPS_POWER_SHIFT_WILKINSON:
763 7 : if (power->shift_type != shift) {
764 4 : power->shift_type = shift;
765 4 : eps->state = EPS_STATE_INITIAL;
766 : }
767 7 : break;
768 0 : default:
769 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
770 : }
771 7 : PetscFunctionReturn(PETSC_SUCCESS);
772 : }
773 :
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.
778 :
779 : Logically Collective
780 :
781 : Input Parameters:
782 : + eps - the eigenproblem solver context
783 : - shift - the type of shift
784 :
785 : Options Database Key:
786 : . -eps_power_shift_type - Sets the shift type (either 'constant' or
787 : 'rayleigh' or 'wilkinson')
788 :
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).
793 :
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.
797 :
798 : Level: advanced
799 :
800 : .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
801 : @*/
802 7 : PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
803 : {
804 7 : PetscFunctionBegin;
805 7 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
806 28 : PetscValidLogicalCollectiveEnum(eps,shift,2);
807 7 : PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
808 7 : PetscFunctionReturn(PETSC_SUCCESS);
809 : }
810 :
811 1 : static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
812 : {
813 1 : EPS_POWER *power = (EPS_POWER*)eps->data;
814 :
815 1 : PetscFunctionBegin;
816 1 : *shift = power->shift_type;
817 1 : PetscFunctionReturn(PETSC_SUCCESS);
818 : }
819 :
820 : /*@
821 : EPSPowerGetShiftType - Gets the type of shifts used during the power
822 : iteration.
823 :
824 : Not Collective
825 :
826 : Input Parameter:
827 : . eps - the eigenproblem solver context
828 :
829 : Output Parameter:
830 : . shift - the type of shift
831 :
832 : Level: advanced
833 :
834 : .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
835 : @*/
836 1 : PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
837 : {
838 1 : PetscFunctionBegin;
839 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
840 1 : PetscAssertPointer(shift,2);
841 1 : PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
842 1 : PetscFunctionReturn(PETSC_SUCCESS);
843 : }
844 :
845 29 : static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
846 : {
847 29 : EPS_POWER *power = (EPS_POWER*)eps->data;
848 :
849 29 : PetscFunctionBegin;
850 29 : if (power->nonlinear != nonlinear) {
851 29 : power->nonlinear = nonlinear;
852 29 : eps->useds = PetscNot(nonlinear);
853 29 : eps->ops->setupsort = nonlinear? NULL: EPSSetUpSort_Default;
854 29 : eps->state = EPS_STATE_INITIAL;
855 : }
856 29 : PetscFunctionReturn(PETSC_SUCCESS);
857 : }
858 :
859 : /*@
860 : EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.
861 :
862 : Logically Collective
863 :
864 : Input Parameters:
865 : + eps - the eigenproblem solver context
866 : - nonlinear - whether the problem is nonlinear or not
867 :
868 : Options Database Key:
869 : . -eps_power_nonlinear - Sets the nonlinear flag
870 :
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).
877 :
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.
884 :
885 : Level: advanced
886 :
887 : .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
888 : @*/
889 29 : PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
890 : {
891 29 : PetscFunctionBegin;
892 29 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
893 116 : PetscValidLogicalCollectiveBool(eps,nonlinear,2);
894 29 : PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
895 29 : PetscFunctionReturn(PETSC_SUCCESS);
896 : }
897 :
898 17 : static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
899 : {
900 17 : EPS_POWER *power = (EPS_POWER*)eps->data;
901 :
902 17 : PetscFunctionBegin;
903 17 : *nonlinear = power->nonlinear;
904 17 : PetscFunctionReturn(PETSC_SUCCESS);
905 : }
906 :
907 : /*@
908 : EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.
909 :
910 : Not Collective
911 :
912 : Input Parameter:
913 : . eps - the eigenproblem solver context
914 :
915 : Output Parameter:
916 : . nonlinear - the nonlinear flag
917 :
918 : Level: advanced
919 :
920 : .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
921 : @*/
922 17 : PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
923 : {
924 17 : PetscFunctionBegin;
925 17 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
926 17 : PetscAssertPointer(nonlinear,2);
927 17 : PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
928 17 : PetscFunctionReturn(PETSC_SUCCESS);
929 : }
930 :
931 12 : static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
932 : {
933 12 : EPS_POWER *power = (EPS_POWER*)eps->data;
934 :
935 12 : PetscFunctionBegin;
936 12 : PetscCheck(power->nonlinear,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
937 12 : power->update = update;
938 12 : eps->state = EPS_STATE_INITIAL;
939 12 : PetscFunctionReturn(PETSC_SUCCESS);
940 : }
941 :
942 : /*@
943 : EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
944 : for nonlinear problems. This potentially has a better convergence.
945 :
946 : Logically Collective
947 :
948 : Input Parameters:
949 : + eps - the eigenproblem solver context
950 : - update - whether the residual is updated monolithically or not
951 :
952 : Options Database Key:
953 : . -eps_power_update - Sets the update flag
954 :
955 : Level: advanced
956 :
957 : .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
958 : @*/
959 12 : PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
960 : {
961 12 : PetscFunctionBegin;
962 12 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
963 48 : PetscValidLogicalCollectiveBool(eps,update,2);
964 12 : PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
965 12 : PetscFunctionReturn(PETSC_SUCCESS);
966 : }
967 :
968 17 : static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
969 : {
970 17 : EPS_POWER *power = (EPS_POWER*)eps->data;
971 :
972 17 : PetscFunctionBegin;
973 17 : *update = power->update;
974 17 : PetscFunctionReturn(PETSC_SUCCESS);
975 : }
976 :
977 : /*@
978 : EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
979 : for nonlinear problems.
980 :
981 : Not Collective
982 :
983 : Input Parameter:
984 : . eps - the eigenproblem solver context
985 :
986 : Output Parameter:
987 : . update - the update flag
988 :
989 : Level: advanced
990 :
991 : .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
992 : @*/
993 17 : PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
994 : {
995 17 : PetscFunctionBegin;
996 17 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
997 17 : PetscAssertPointer(update,2);
998 17 : PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
999 17 : PetscFunctionReturn(PETSC_SUCCESS);
1000 : }
1001 :
1002 29 : static PetscErrorCode EPSPowerSetSignNormalization_Power(EPS eps,PetscBool sign_normalization)
1003 : {
1004 29 : EPS_POWER *power = (EPS_POWER*)eps->data;
1005 :
1006 29 : PetscFunctionBegin;
1007 29 : power->sign_normalization = sign_normalization;
1008 29 : PetscFunctionReturn(PETSC_SUCCESS);
1009 : }
1010 :
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.
1016 :
1017 : Logically Collective
1018 :
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
1023 :
1024 : Options Database Key:
1025 : . -eps_power_sign_normalization - Sets the sign normalization flag
1026 :
1027 : Level: advanced
1028 :
1029 : .seealso: EPSPowerGetSignNormalization()
1030 : @*/
1031 29 : PetscErrorCode EPSPowerSetSignNormalization(EPS eps,PetscBool sign_normalization)
1032 : {
1033 29 : PetscFunctionBegin;
1034 29 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1035 116 : PetscValidLogicalCollectiveBool(eps,sign_normalization,2);
1036 29 : PetscTryMethod(eps,"EPSPowerSetSignNormalization_C",(EPS,PetscBool),(eps,sign_normalization));
1037 29 : PetscFunctionReturn(PETSC_SUCCESS);
1038 : }
1039 :
1040 12 : static PetscErrorCode EPSPowerGetSignNormalization_Power(EPS eps,PetscBool *sign_normalization)
1041 : {
1042 12 : EPS_POWER *power = (EPS_POWER*)eps->data;
1043 :
1044 12 : PetscFunctionBegin;
1045 12 : *sign_normalization = power->sign_normalization;
1046 12 : PetscFunctionReturn(PETSC_SUCCESS);
1047 : }
1048 :
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.
1054 :
1055 : Not Collective
1056 :
1057 : Input Parameter:
1058 : . eps - the eigenproblem solver context
1059 :
1060 : Output Parameter:
1061 : . sign_normalization - the sign normalization flag
1062 :
1063 : Level: advanced
1064 :
1065 : .seealso: EPSPowerSetSignNormalization()
1066 : @*/
1067 12 : PetscErrorCode EPSPowerGetSignNormalization(EPS eps,PetscBool *sign_normalization)
1068 : {
1069 12 : PetscFunctionBegin;
1070 12 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1071 12 : PetscAssertPointer(sign_normalization,2);
1072 12 : PetscUseMethod(eps,"EPSPowerGetSignNormalization_C",(EPS,PetscBool*),(eps,sign_normalization));
1073 12 : PetscFunctionReturn(PETSC_SUCCESS);
1074 : }
1075 :
1076 0 : static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
1077 : {
1078 0 : EPS_POWER *power = (EPS_POWER*)eps->data;
1079 :
1080 0 : PetscFunctionBegin;
1081 0 : PetscCall(PetscObjectReference((PetscObject)snes));
1082 0 : PetscCall(SNESDestroy(&power->snes));
1083 0 : power->snes = snes;
1084 0 : eps->state = EPS_STATE_INITIAL;
1085 0 : PetscFunctionReturn(PETSC_SUCCESS);
1086 : }
1087 :
1088 : /*@
1089 : EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
1090 : eigenvalue solver (to be used in nonlinear inverse iteration).
1091 :
1092 : Collective
1093 :
1094 : Input Parameters:
1095 : + eps - the eigenvalue solver
1096 : - snes - the nonlinear solver object
1097 :
1098 : Level: advanced
1099 :
1100 : .seealso: EPSPowerGetSNES()
1101 : @*/
1102 0 : PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1103 : {
1104 0 : PetscFunctionBegin;
1105 0 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1106 0 : PetscValidHeaderSpecific(snes,SNES_CLASSID,2);
1107 0 : PetscCheckSameComm(eps,1,snes,2);
1108 0 : PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1109 0 : PetscFunctionReturn(PETSC_SUCCESS);
1110 : }
1111 :
1112 41 : static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1113 : {
1114 41 : EPS_POWER *power = (EPS_POWER*)eps->data;
1115 :
1116 41 : PetscFunctionBegin;
1117 41 : if (!power->snes) {
1118 29 : PetscCall(SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes));
1119 29 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1));
1120 29 : PetscCall(SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix));
1121 29 : PetscCall(SNESAppendOptionsPrefix(power->snes,"eps_power_"));
1122 29 : PetscCall(PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options));
1123 : }
1124 41 : *snes = power->snes;
1125 41 : PetscFunctionReturn(PETSC_SUCCESS);
1126 : }
1127 :
1128 : /*@
1129 : EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1130 : with the eigenvalue solver.
1131 :
1132 : Not Collective
1133 :
1134 : Input Parameter:
1135 : . eps - the eigenvalue solver
1136 :
1137 : Output Parameter:
1138 : . snes - the nonlinear solver object
1139 :
1140 : Level: advanced
1141 :
1142 : .seealso: EPSPowerSetSNES()
1143 : @*/
1144 41 : PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1145 : {
1146 41 : PetscFunctionBegin;
1147 41 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1148 41 : PetscAssertPointer(snes,2);
1149 41 : PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1150 41 : PetscFunctionReturn(PETSC_SUCCESS);
1151 : }
1152 :
1153 45 : static PetscErrorCode EPSReset_Power(EPS eps)
1154 : {
1155 45 : EPS_POWER *power = (EPS_POWER*)eps->data;
1156 :
1157 45 : PetscFunctionBegin;
1158 45 : if (power->snes) PetscCall(SNESReset(power->snes));
1159 45 : PetscFunctionReturn(PETSC_SUCCESS);
1160 : }
1161 :
1162 44 : static PetscErrorCode EPSDestroy_Power(EPS eps)
1163 : {
1164 44 : EPS_POWER *power = (EPS_POWER*)eps->data;
1165 :
1166 44 : PetscFunctionBegin;
1167 44 : if (power->nonlinear) PetscCall(SNESDestroy(&power->snes));
1168 44 : PetscCall(PetscFree(eps->data));
1169 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL));
1170 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL));
1171 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL));
1172 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL));
1173 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL));
1174 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL));
1175 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",NULL));
1176 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",NULL));
1177 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL));
1178 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL));
1179 44 : PetscFunctionReturn(PETSC_SUCCESS);
1180 : }
1181 :
1182 0 : static PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1183 : {
1184 0 : EPS_POWER *power = (EPS_POWER*)eps->data;
1185 0 : PetscBool isascii;
1186 :
1187 0 : PetscFunctionBegin;
1188 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1189 0 : if (isascii) {
1190 0 : if (power->nonlinear) {
1191 0 : if (power->sign_normalization) PetscCall(PetscViewerASCIIPrintf(viewer," normalizing Bx by the sign of the first nonzero element\n"));
1192 0 : else PetscCall(PetscViewerASCIIPrintf(viewer," not normalizing Bx by the sign of the first nonzero element\n"));
1193 0 : PetscCall(PetscViewerASCIIPrintf(viewer," using nonlinear inverse iteration\n"));
1194 0 : if (power->update) PetscCall(PetscViewerASCIIPrintf(viewer," updating the residual monolithically\n"));
1195 0 : if (!power->snes) PetscCall(EPSPowerGetSNES(eps,&power->snes));
1196 0 : PetscCall(PetscViewerASCIIPushTab(viewer));
1197 0 : PetscCall(SNESView(power->snes,viewer));
1198 0 : PetscCall(PetscViewerASCIIPopTab(viewer));
1199 0 : } else PetscCall(PetscViewerASCIIPrintf(viewer," %s shifts\n",EPSPowerShiftTypes[power->shift_type]));
1200 : }
1201 0 : PetscFunctionReturn(PETSC_SUCCESS);
1202 : }
1203 :
1204 36 : static PetscErrorCode EPSComputeVectors_Power(EPS eps)
1205 : {
1206 36 : EPS_POWER *power = (EPS_POWER*)eps->data;
1207 :
1208 36 : PetscFunctionBegin;
1209 36 : if (eps->twosided) {
1210 3 : PetscCall(EPSComputeVectors_Twosided(eps));
1211 3 : PetscCall(BVNormalize(eps->V,NULL));
1212 3 : PetscCall(BVNormalize(eps->W,NULL));
1213 33 : } else if (!power->nonlinear) PetscCall(EPSComputeVectors_Schur(eps));
1214 36 : PetscFunctionReturn(PETSC_SUCCESS);
1215 : }
1216 :
1217 92 : static PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1218 : {
1219 92 : EPS_POWER *power = (EPS_POWER*)eps->data;
1220 92 : KSP ksp;
1221 92 : PC pc;
1222 :
1223 92 : PetscFunctionBegin;
1224 92 : if (power->nonlinear) {
1225 60 : eps->categ=EPS_CATEGORY_PRECOND;
1226 60 : 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 60 : PetscCall(STSetType(eps->st,STPRECOND));
1231 60 : PetscCall(KSPGetPC(ksp,&pc));
1232 60 : PetscCall(PCSetType(pc,PCNONE));
1233 : }
1234 92 : PetscFunctionReturn(PETSC_SUCCESS);
1235 : }
1236 :
1237 44 : SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1238 : {
1239 44 : EPS_POWER *ctx;
1240 :
1241 44 : PetscFunctionBegin;
1242 44 : PetscCall(PetscNew(&ctx));
1243 44 : eps->data = (void*)ctx;
1244 :
1245 44 : eps->useds = PETSC_TRUE;
1246 44 : eps->categ = EPS_CATEGORY_OTHER;
1247 :
1248 44 : eps->ops->setup = EPSSetUp_Power;
1249 44 : eps->ops->setupsort = EPSSetUpSort_Default;
1250 44 : eps->ops->setfromoptions = EPSSetFromOptions_Power;
1251 44 : eps->ops->reset = EPSReset_Power;
1252 44 : eps->ops->destroy = EPSDestroy_Power;
1253 44 : eps->ops->view = EPSView_Power;
1254 44 : eps->ops->backtransform = EPSBackTransform_Power;
1255 44 : eps->ops->computevectors = EPSComputeVectors_Power;
1256 44 : eps->ops->setdefaultst = EPSSetDefaultST_Power;
1257 44 : eps->stopping = EPSStopping_Power;
1258 44 : ctx->sign_normalization = PETSC_TRUE;
1259 :
1260 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power));
1261 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power));
1262 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power));
1263 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power));
1264 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power));
1265 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power));
1266 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSignNormalization_C",EPSPowerSetSignNormalization_Power));
1267 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSignNormalization_C",EPSPowerGetSignNormalization_Power));
1268 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power));
1269 44 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power));
1270 44 : PetscFunctionReturn(PETSC_SUCCESS);
1271 : }
|
