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