Actual source code: power.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "power"
13: Method: Power Iteration
15: Algorithm:
17: This solver implements the power iteration for finding dominant
18: eigenpairs. It also includes the following well-known methods:
19: - Inverse Iteration: when used in combination with shift-and-invert
20: spectral transformation.
21: - Rayleigh Quotient Iteration (RQI): also with shift-and-invert plus
22: a variable shift.
24: It can also be used for nonlinear inverse iteration on the problem
25: A(x)*x=lambda*B(x)*x, where A and B are not constant but depend on x.
27: References:
29: [1] "Single Vector Iteration Methods in SLEPc", SLEPc Technical Report
30: STR-2, available at http://slepc.upv.es.
31: */
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>
39: static PetscErrorCode EPSPowerFormFunction_Update(SNES,Vec,Vec,void*);
40: static PetscErrorCode SNESLineSearchPostheckFunction(SNESLineSearch linesearch,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx);
41: PetscErrorCode EPSSolve_Power(EPS);
42: PetscErrorCode EPSSolve_TS_Power(EPS);
44: typedef struct {
45: EPSPowerShiftType shift_type;
46: PetscBool nonlinear;
47: PetscBool update;
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: } EPS_POWER;
56: PetscErrorCode EPSSetUp_Power(EPS eps)
57: {
59: EPS_POWER *power = (EPS_POWER*)eps->data;
60: PetscBool flg,istrivial;
61: STMatMode mode;
62: Mat A,B;
63: Vec res;
64: PetscContainer container;
65: PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
66: PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
67: void *ctx;
68: SNESLineSearch linesearch;
71: if (eps->ncv) {
72: if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
73: } else eps->ncv = eps->nev;
74: if (eps->mpd) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
75: if (!eps->max_it) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
76: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
77: if (eps->which!=EPS_LARGEST_MAGNITUDE && eps->which !=EPS_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
78: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
79: if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
80: PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STCAYLEY,"");
81: if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts only allowed in shift-and-invert or Cayley ST");
82: STGetMatMode(eps->st,&mode);
83: if (mode == ST_MATMODE_INPLACE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
84: }
85: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
86: if (eps->balance!=EPS_BALANCE_NONE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Balancing not supported in this solver");
87: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
88: RGIsTrivial(eps->rg,&istrivial);
89: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");
90: EPSAllocateSolution(eps,0);
91: EPS_SetInnerProduct(eps);
93: if (power->nonlinear) {
94: if (eps->nev>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
95: EPSSetWorkVecs(eps,4);
97: /* set up SNES solver */
98: if (!power->snes) { EPSPowerGetSNES(eps,&power->snes); }
99: EPSGetOperators(eps,&A,&B);
101: PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA);
102: if (!formFunctionA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
103: PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container);
104: if (container) {
105: PetscContainerGetPointer(container,&ctx);
106: } else ctx = NULL;
107: MatCreateVecs(A,&res,NULL);
108: power->formFunctionA = formFunctionA;
109: power->formFunctionActx = ctx;
110: if (power->update) {
111: SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx);
112: PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB);
113: }
114: else { SNESSetFunction(power->snes,res,formFunctionA,ctx); }
115: VecDestroy(&res);
117: PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA);
118: if (!formJacobianA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
119: PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container);
120: if (container) {
121: PetscContainerGetPointer(container,&ctx);
122: } else ctx = NULL;
123: SNESSetJacobian(power->snes,A,A,formJacobianA,ctx);
124: SNESSetFromOptions(power->snes);
125: SNESGetLineSearch(power->snes,&linesearch);
126: if (power->update) {
127: SNESLineSearchSetPostCheck(linesearch,SNESLineSearchPostheckFunction,ctx);
128: }
129: SNESSetUp(power->snes);
130: if (B) {
131: PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB);
132: PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container);
133: if (power->formFunctionB && container) {
134: PetscContainerGetPointer(container,&power->formFunctionBctx);
135: } else power->formFunctionBctx = NULL;
136: }
137: } else {
138: if (eps->twosided) {
139: EPSSetWorkVecs(eps,3);
140: } else {
141: EPSSetWorkVecs(eps,2);
142: }
143: DSSetType(eps->ds,DSNHEP);
144: DSAllocate(eps->ds,eps->nev);
145: }
146: /* dispatch solve method */
147: if (eps->twosided) {
148: if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
149: if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
150: eps->ops->solve = EPSSolve_TS_Power;
151: } else eps->ops->solve = EPSSolve_Power;
152: return(0);
153: }
155: /*
156: Normalize a vector x with respect to a given norm as well as the
157: sign of the first entry.
158: */
159: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscScalar *sign)
160: {
161: PetscErrorCode ierr;
162: PetscScalar alpha = 1.0;
163: PetscMPIInt rank;
164: PetscReal absv;
165: const PetscScalar *xx;
168: MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank);
169: if (!rank) {
170: VecGetArrayRead(x,&xx);
171: absv = PetscAbsScalar(*xx);
172: if (absv>10*PETSC_MACHINE_EPSILON) alpha = *xx/absv;
173: VecRestoreArrayRead(x,&xx);
174: }
175: MPI_Bcast(&alpha,1,MPIU_SCALAR,0,PetscObjectComm((PetscObject)x));
176: if (sign) *sign = alpha;
177: alpha *= norm;
178: VecScale(x,1.0/alpha);
179: return(0);
180: }
182: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
183: {
185: EPS_POWER *power = (EPS_POWER*)eps->data;
186: Mat B;
189: STResetMatrixState(eps->st);
190: EPSGetOperators(eps,NULL,&B);
191: if (B) {
192: if (power->formFunctionB) {
193: (*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx);
194: } else {
195: MatMult(B,x,Bx);
196: }
197: } else {
198: VecCopy(x,Bx);
199: }
200: return(0);
201: }
203: static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
204: {
206: EPS_POWER *power = (EPS_POWER*)eps->data;
207: Mat A;
210: STResetMatrixState(eps->st);
211: EPSGetOperators(eps,&A,NULL);
212: if (A) {
213: if (power->formFunctionA) {
214: (*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx);
215: } else {
216: MatMult(A,x,Ax);
217: }
218: } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem\n");
219: return(0);
220: }
222: static PetscErrorCode SNESLineSearchPostheckFunction(SNESLineSearch linesearch,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
223: {
225: SNES snes;
226: EPS eps;
227: Vec oldx;
230: SNESLineSearchGetSNES(linesearch,&snes);
231: PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
232: if (!eps) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No composed EPS");
233: oldx = eps->work[3];
234: VecCopy(x,oldx);
235: return(0);
236: }
238: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
239: {
241: EPS eps;
242: PetscReal bx;
243: Vec Bx;
244: EPS_POWER *power;
247: PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
248: if (!eps) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No composed EPS");
249: power = (EPS_POWER*)eps->data;
250: Bx = eps->work[2];
251: if (power->formFunctionAB) {
252: (*power->formFunctionAB)(snes,x,y,Bx,ctx);
253: } else {
254: EPSPowerUpdateFunctionA(eps,x,y);
255: EPSPowerUpdateFunctionB(eps,x,Bx);
256: }
257: VecNorm(Bx,NORM_2,&bx);
258: Normalize(Bx,bx,NULL);
259: VecAXPY(y,-1.0,Bx);
260: return(0);
261: }
263: /*
264: Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
265: */
266: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
267: {
269: EPS_POWER *power = (EPS_POWER*)eps->data;
270: Vec Bx;
273: VecCopy(x,y);
274: if (power->update) {
275: SNESSolve(power->snes,NULL,y);
276: } else {
277: Bx = eps->work[2];
278: SNESSolve(power->snes,Bx,y);
279: }
280: return(0);
281: }
283: /*
284: Use nonlinear inverse power to compute an initial guess
285: */
286: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
287: {
288: EPS powereps;
289: Mat A,B;
290: Vec v1,v2;
291: SNES snes;
292: DM dm,newdm;
296: EPSCreate(PetscObjectComm((PetscObject)eps),&powereps);
297: EPSGetOperators(eps,&A,&B);
298: EPSSetType(powereps,EPSPOWER);
299: EPSSetOperators(powereps,A,B);
300: EPSSetTolerances(powereps,1e-6,4);
301: EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix);
302: EPSAppendOptionsPrefix(powereps,"init_");
303: EPSSetProblemType(powereps,EPS_GNHEP);
304: EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE);
305: EPSPowerSetNonlinear(powereps,PETSC_TRUE);
306: STSetType(powereps->st,STSINVERT);
307: /* attach dm to initial solve */
308: EPSPowerGetSNES(eps,&snes);
309: SNESGetDM(snes,&dm);
310: /* use dmshell to temporarily store snes context */
311: DMCreate(PetscObjectComm((PetscObject)eps),&newdm);
312: DMSetType(newdm,DMSHELL);
313: DMSetUp(newdm);
314: DMCopyDMSNES(dm,newdm);
315: EPSPowerGetSNES(powereps,&snes);
316: SNESSetDM(snes,dm);
317: EPSSetFromOptions(powereps);
318: EPSSolve(powereps);
319: BVGetColumn(eps->V,0,&v2);
320: BVGetColumn(powereps->V,0,&v1);
321: VecCopy(v1,v2);
322: BVRestoreColumn(powereps->V,0,&v1);
323: BVRestoreColumn(eps->V,0,&v2);
324: EPSDestroy(&powereps);
325: /* restore context back to the old nonlinear solver */
326: DMCopyDMSNES(newdm,dm);
327: DMDestroy(&newdm);
328: return(0);
329: }
331: PetscErrorCode EPSSolve_Power(EPS eps)
332: {
333: #if defined(SLEPC_MISSING_LAPACK_LAEV2)
335: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LAEV2 - Lapack routine is unavailable");
336: #else
337: PetscErrorCode ierr;
338: EPS_POWER *power = (EPS_POWER*)eps->data;
339: PetscInt k,ld;
340: Vec v,y,e,Bx;
341: Mat A;
342: KSP ksp;
343: PetscReal relerr,norm,norm1,rt1,rt2,cs1;
344: PetscScalar theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
345: PetscBool breakdown;
346: KSPConvergedReason reason;
347: SNESConvergedReason snesreason;
350: e = eps->work[0];
351: y = eps->work[1];
352: if (power->nonlinear) Bx = eps->work[2];
353: else Bx = NULL;
355: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
356: if (power->nonlinear) {
357: PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps);
358: if (power->update) {
359: EPSPowerComputeInitialGuess_Update(eps);
360: }
361: } else {
362: DSGetLeadingDimension(eps->ds,&ld);
363: }
364: if (!power->update) {
365: EPSGetStartVector(eps,0,NULL);
366: }
367: if (power->nonlinear) {
368: BVGetColumn(eps->V,0,&v);
369: EPSPowerUpdateFunctionB(eps,v,Bx);
370: VecNorm(Bx,NORM_2,&norm);
371: Normalize(Bx,norm,NULL);
372: BVRestoreColumn(eps->V,0,&v);
373: }
375: STGetShift(eps->st,&sigma); /* original shift */
376: rho = sigma;
378: while (eps->reason == EPS_CONVERGED_ITERATING) {
379: eps->its++;
380: k = eps->nconv;
382: /* y = OP v */
383: BVGetColumn(eps->V,k,&v);
384: if (power->nonlinear) {
385: VecCopy(v,eps->work[3]);
386: EPSPowerApply_SNES(eps,v,y);
387: VecCopy(eps->work[3],v);
388: } else {
389: STApply(eps->st,v,y);
390: }
391: BVRestoreColumn(eps->V,k,&v);
393: /* purge previously converged eigenvectors */
394: if (power->nonlinear) {
395: EPSPowerUpdateFunctionB(eps,y,Bx);
396: VecNorm(Bx,NORM_2,&norm);
397: Normalize(Bx,norm,&sign);
398: } else {
399: DSGetArray(eps->ds,DS_MAT_A,&T);
400: BVSetActiveColumns(eps->V,0,k);
401: BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL);
402: }
404: /* theta = (v,y)_B */
405: BVSetActiveColumns(eps->V,k,k+1);
406: BVDotVec(eps->V,y,&theta);
407: if (!power->nonlinear) {
408: T[k+k*ld] = theta;
409: DSRestoreArray(eps->ds,DS_MAT_A,&T);
410: }
412: if (power->nonlinear) theta = norm*sign;
414: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
416: /* approximate eigenvalue is the Rayleigh quotient */
417: eps->eigr[eps->nconv] = theta;
419: /* compute relative error as ||y-theta v||_2/|theta| */
420: VecCopy(y,e);
421: BVGetColumn(eps->V,k,&v);
422: VecAXPY(e,power->nonlinear?-1.0:-theta,v);
423: BVRestoreColumn(eps->V,k,&v);
424: VecNorm(e,NORM_2,&relerr);
425: relerr /= PetscAbsScalar(theta);
427: } else { /* RQI */
429: /* delta = ||y||_B */
430: delta = norm;
432: /* compute relative error */
433: if (rho == 0.0) relerr = PETSC_MAX_REAL;
434: else relerr = 1.0 / (norm*PetscAbsScalar(rho));
436: /* approximate eigenvalue is the shift */
437: eps->eigr[eps->nconv] = rho;
439: /* compute new shift */
440: if (relerr<eps->tol) {
441: rho = sigma; /* if converged, restore original shift */
442: STSetShift(eps->st,rho);
443: } else {
444: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
445: if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
446: /* beta1 is the norm of the residual associated with R(v) */
447: BVGetColumn(eps->V,k,&v);
448: VecAXPY(v,-PetscConj(theta)/(delta*delta),y);
449: BVRestoreColumn(eps->V,k,&v);
450: BVScaleColumn(eps->V,k,1.0/delta);
451: BVNormColumn(eps->V,k,NORM_2,&norm1);
452: beta1 = norm1;
454: /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
455: STGetMatrix(eps->st,0,&A);
456: BVGetColumn(eps->V,k,&v);
457: MatMult(A,v,e);
458: VecDot(v,e,&alpha2);
459: BVRestoreColumn(eps->V,k,&v);
460: alpha2 = alpha2 / (beta1 * beta1);
462: /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
463: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
464: PetscStackCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
465: PetscFPTrapPop();
466: if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
467: else rho = rt2;
468: }
469: /* update operator according to new shift */
470: KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
471: STSetShift(eps->st,rho);
472: KSPGetConvergedReason(ksp,&reason);
473: if (reason) {
474: PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
475: rho *= 1+10*PETSC_MACHINE_EPSILON;
476: STSetShift(eps->st,rho);
477: KSPGetConvergedReason(ksp,&reason);
478: if (reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Second factorization failed");
479: }
480: KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
481: }
482: }
483: eps->errest[eps->nconv] = relerr;
485: /* normalize */
486: if (!power->nonlinear) { Normalize(y,norm,NULL); }
487: BVInsertVec(eps->V,k,y);
489: if (power->update) {
490: SNESGetConvergedReason(power->snes,&snesreason);
491: }
492: /* if relerr<tol, accept eigenpair */
493: if (relerr<eps->tol || (power->update && snesreason > 0)) {
494: eps->nconv = eps->nconv + 1;
495: if (eps->nconv<eps->nev) {
496: EPSGetStartVector(eps,eps->nconv,&breakdown);
497: if (breakdown) {
498: eps->reason = EPS_DIVERGED_BREAKDOWN;
499: PetscInfo(eps,"Unable to generate more start vectors\n");
500: break;
501: }
502: }
503: }
504: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
505: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
506: }
508: if (power->nonlinear) {
509: PetscObjectCompose((PetscObject)power->snes,"eps",NULL);
510: } else {
511: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
512: DSSetState(eps->ds,DS_STATE_RAW);
513: }
514: return(0);
515: #endif
516: }
518: PetscErrorCode EPSSolve_TS_Power(EPS eps)
519: {
520: PetscErrorCode ierr;
521: EPS_POWER *power = (EPS_POWER*)eps->data;
522: PetscInt k,ld;
523: Vec v,w,y,e,z;
524: KSP ksp;
525: PetscReal relerr,relerrl,delta;
526: PetscScalar theta,rho,alpha,sigma;
527: PetscBool breakdown;
528: KSPConvergedReason reason;
531: e = eps->work[0];
532: v = eps->work[1];
533: w = eps->work[2];
535: if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
536: DSGetLeadingDimension(eps->ds,&ld);
537: EPSGetStartVector(eps,0,NULL);
538: BVSetRandomColumn(eps->W,0);
539: BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL);
540: BVCopyVec(eps->V,0,v);
541: BVCopyVec(eps->W,0,w);
542: STGetShift(eps->st,&sigma); /* original shift */
543: rho = sigma;
545: while (eps->reason == EPS_CONVERGED_ITERATING) {
546: eps->its++;
547: k = eps->nconv;
549: /* y = OP v, z = OP' w */
550: BVGetColumn(eps->V,k,&y);
551: STApply(eps->st,v,y);
552: BVRestoreColumn(eps->V,k,&y);
553: BVGetColumn(eps->W,k,&z);
554: VecConjugate(w);
555: STApplyTranspose(eps->st,w,z);
556: VecConjugate(w);
557: VecConjugate(z);
558: BVRestoreColumn(eps->W,k,&z);
560: /* purge previously converged eigenvectors */
561: BVBiorthogonalizeColumn(eps->V,eps->W,k);
563: /* theta = (w,y)_B */
564: BVSetActiveColumns(eps->V,k,k+1);
565: BVDotVec(eps->V,w,&theta);
566: theta = PetscConj(theta);
568: if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */
570: /* approximate eigenvalue is the Rayleigh quotient */
571: eps->eigr[eps->nconv] = theta;
573: /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
574: BVCopyVec(eps->V,k,e);
575: VecAXPY(e,-theta,v);
576: VecNorm(e,NORM_2,&relerr);
577: BVCopyVec(eps->W,k,e);
578: VecAXPY(e,-PetscConj(theta),w);
579: VecNorm(e,NORM_2,&relerrl);
580: relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
581: }
583: /* normalize */
584: BVSetActiveColumns(eps->V,k,k+1);
585: BVGetColumn(eps->W,k,&z);
586: BVDotVec(eps->V,z,&alpha);
587: BVRestoreColumn(eps->W,k,&z);
588: delta = PetscSqrtReal(PetscAbsScalar(alpha));
589: if (delta==0.0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Breakdown in two-sided Power/RQI");
590: BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta));
591: BVScaleColumn(eps->W,k,1.0/delta);
592: BVCopyVec(eps->V,k,v);
593: BVCopyVec(eps->W,k,w);
595: if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */
597: /* compute relative error */
598: if (rho == 0.0) relerr = PETSC_MAX_REAL;
599: else relerr = 1.0 / PetscAbsScalar(delta*rho);
601: /* approximate eigenvalue is the shift */
602: eps->eigr[eps->nconv] = rho;
604: /* compute new shift */
605: if (relerr<eps->tol) {
606: rho = sigma; /* if converged, restore original shift */
607: STSetShift(eps->st,rho);
608: } else {
609: rho = rho + PetscConj(theta)/(delta*delta); /* Rayleigh quotient R(v) */
610: /* update operator according to new shift */
611: KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
612: STSetShift(eps->st,rho);
613: KSPGetConvergedReason(ksp,&reason);
614: if (reason) {
615: PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
616: rho *= 1+10*PETSC_MACHINE_EPSILON;
617: STSetShift(eps->st,rho);
618: KSPGetConvergedReason(ksp,&reason);
619: if (reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Second factorization failed");
620: }
621: KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
622: }
623: }
624: eps->errest[eps->nconv] = relerr;
626: /* if relerr<tol, accept eigenpair */
627: if (relerr<eps->tol) {
628: eps->nconv = eps->nconv + 1;
629: if (eps->nconv<eps->nev) {
630: EPSGetStartVector(eps,eps->nconv,&breakdown);
631: if (breakdown) {
632: eps->reason = EPS_DIVERGED_BREAKDOWN;
633: PetscInfo(eps,"Unable to generate more start vectors\n");
634: break;
635: }
636: BVSetRandomColumn(eps->W,eps->nconv);
637: BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL);
638: }
639: }
640: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
641: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
642: }
644: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
645: DSSetState(eps->ds,DS_STATE_RAW);
646: return(0);
647: }
649: PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
650: {
652: EPS_POWER *power = (EPS_POWER*)eps->data;
653: SNESConvergedReason snesreason;
656: if (power->update) {
657: SNESGetConvergedReason(power->snes,&snesreason);
658: if (snesreason < 0) {
659: *reason = EPS_DIVERGED_BREAKDOWN;
660: return(0);
661: }
662: }
663: EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx);
664: return(0);
665: }
667: PetscErrorCode EPSBackTransform_Power(EPS eps)
668: {
670: EPS_POWER *power = (EPS_POWER*)eps->data;
673: if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
674: else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) {
675: EPSBackTransform_Default(eps);
676: }
677: return(0);
678: }
680: PetscErrorCode EPSSetFromOptions_Power(PetscOptionItems *PetscOptionsObject,EPS eps)
681: {
682: PetscErrorCode ierr;
683: EPS_POWER *power = (EPS_POWER*)eps->data;
684: PetscBool flg,val;
685: EPSPowerShiftType shift;
688: PetscOptionsHead(PetscOptionsObject,"EPS Power Options");
690: PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg);
691: if (flg) { EPSPowerSetShiftType(eps,shift); }
693: PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg);
694: if (flg) { EPSPowerSetNonlinear(eps,val); }
696: PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg);
697: if (flg) { EPSPowerSetUpdate(eps,val); }
699: PetscOptionsTail();
700: return(0);
701: }
703: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
704: {
705: EPS_POWER *power = (EPS_POWER*)eps->data;
708: switch (shift) {
709: case EPS_POWER_SHIFT_CONSTANT:
710: case EPS_POWER_SHIFT_RAYLEIGH:
711: case EPS_POWER_SHIFT_WILKINSON:
712: if (power->shift_type != shift) {
713: power->shift_type = shift;
714: eps->state = EPS_STATE_INITIAL;
715: }
716: break;
717: default:
718: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
719: }
720: return(0);
721: }
723: /*@
724: EPSPowerSetShiftType - Sets the type of shifts used during the power
725: iteration. This can be used to emulate the Rayleigh Quotient Iteration
726: (RQI) method.
728: Logically Collective on EPS
730: Input Parameters:
731: + eps - the eigenproblem solver context
732: - shift - the type of shift
734: Options Database Key:
735: . -eps_power_shift_type - Sets the shift type (either 'constant' or
736: 'rayleigh' or 'wilkinson')
738: Notes:
739: By default, shifts are constant (EPS_POWER_SHIFT_CONSTANT) and the iteration
740: is the simple power method (or inverse iteration if a shift-and-invert
741: transformation is being used).
743: A variable shift can be specified (EPS_POWER_SHIFT_RAYLEIGH or
744: EPS_POWER_SHIFT_WILKINSON). In this case, the iteration behaves rather like
745: a cubic converging method such as RQI.
747: Level: advanced
749: .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
750: @*/
751: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
752: {
758: PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
759: return(0);
760: }
762: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
763: {
764: EPS_POWER *power = (EPS_POWER*)eps->data;
767: *shift = power->shift_type;
768: return(0);
769: }
771: /*@
772: EPSPowerGetShiftType - Gets the type of shifts used during the power
773: iteration.
775: Not Collective
777: Input Parameter:
778: . eps - the eigenproblem solver context
780: Input Parameter:
781: . shift - the type of shift
783: Level: advanced
785: .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
786: @*/
787: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
788: {
794: PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
795: return(0);
796: }
798: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
799: {
800: EPS_POWER *power = (EPS_POWER*)eps->data;
803: if (power->nonlinear != nonlinear) {
804: power->nonlinear = nonlinear;
805: eps->useds = PetscNot(nonlinear);
806: eps->state = EPS_STATE_INITIAL;
807: }
808: return(0);
809: }
811: /*@
812: EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.
814: Logically Collective on EPS
816: Input Parameters:
817: + eps - the eigenproblem solver context
818: - nonlinear - whether the problem is nonlinear or not
820: Options Database Key:
821: . -eps_power_nonlinear - Sets the nonlinear flag
823: Notes:
824: If this flag is set, the solver assumes that the problem is nonlinear,
825: that is, the operators that define the eigenproblem are not constant
826: matrices, but depend on the eigenvector: A(x)*x=lambda*B(x)*x. This is
827: different from the case of nonlinearity with respect to the eigenvalue
828: (use the NEP solver class for this kind of problems).
830: The way in which nonlinear operators are specified is very similar to
831: the case of PETSc's SNES solver. The difference is that the callback
832: functions are provided via composed functions "formFunction" and
833: "formJacobian" in each of the matrix objects passed as arguments of
834: EPSSetOperators(). The application context required for these functions
835: can be attached via a composed PetscContainer.
837: Level: advanced
839: .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
840: @*/
841: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
842: {
848: PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
849: return(0);
850: }
852: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
853: {
854: EPS_POWER *power = (EPS_POWER*)eps->data;
857: *nonlinear = power->nonlinear;
858: return(0);
859: }
861: /*@
862: EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.
864: Not Collective
866: Input Parameter:
867: . eps - the eigenproblem solver context
869: Input Parameter:
870: . nonlinear - the nonlinear flag
872: Level: advanced
874: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
875: @*/
876: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
877: {
883: PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
884: return(0);
885: }
887: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
888: {
889: EPS_POWER *power = (EPS_POWER*)eps->data;
892: if (!power->nonlinear) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems \n");
893: power->update = update;
894: eps->state = EPS_STATE_INITIAL;
895: return(0);
896: }
898: /*@
899: EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
900: for nonlinear problems. This potentially has a better convergence.
902: Logically Collective on EPS
904: Input Parameters:
905: + eps - the eigenproblem solver context
906: - update - whether the residual is updated monolithically or not
908: Options Database Key:
909: . -eps_power_update - Sets the update flag
911: Level: advanced
913: .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
914: @*/
915: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
916: {
922: PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
923: return(0);
924: }
926: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
927: {
928: EPS_POWER *power = (EPS_POWER*)eps->data;
931: *update = power->update;
932: return(0);
933: }
935: /*@
936: EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
937: for nonlinear problems.
939: Not Collective
941: Input Parameter:
942: . eps - the eigenproblem solver context
944: Input Parameter:
945: . update - the update flag
947: Level: advanced
949: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
950: @*/
951: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
952: {
958: PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
959: return(0);
960: }
962: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
963: {
965: EPS_POWER *power = (EPS_POWER*)eps->data;
968: PetscObjectReference((PetscObject)snes);
969: SNESDestroy(&power->snes);
970: power->snes = snes;
971: PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
972: eps->state = EPS_STATE_INITIAL;
973: return(0);
974: }
976: /*@
977: EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
978: eigenvalue solver (to be used in nonlinear inverse iteration).
980: Collective on EPS
982: Input Parameters:
983: + eps - the eigenvalue solver
984: - snes - the nonlinear solver object
986: Level: advanced
988: .seealso: EPSPowerGetSNES()
989: @*/
990: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
991: {
998: PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
999: return(0);
1000: }
1002: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1003: {
1005: EPS_POWER *power = (EPS_POWER*)eps->data;
1008: if (!power->snes) {
1009: SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes);
1010: PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1);
1011: SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix);
1012: SNESAppendOptionsPrefix(power->snes,"eps_power_");
1013: PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
1014: PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options);
1015: }
1016: *snes = power->snes;
1017: return(0);
1018: }
1020: /*@
1021: EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1022: with the eigenvalue solver.
1024: Not Collective
1026: Input Parameter:
1027: . eps - the eigenvalue solver
1029: Output Parameter:
1030: . snes - the nonlinear solver object
1032: Level: advanced
1034: .seealso: EPSPowerSetSNES()
1035: @*/
1036: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1037: {
1043: PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1044: return(0);
1045: }
1047: PetscErrorCode EPSReset_Power(EPS eps)
1048: {
1050: EPS_POWER *power = (EPS_POWER*)eps->data;
1053: if (power->snes) { SNESReset(power->snes); }
1054: return(0);
1055: }
1057: PetscErrorCode EPSDestroy_Power(EPS eps)
1058: {
1060: EPS_POWER *power = (EPS_POWER*)eps->data;
1063: if (power->nonlinear) {
1064: SNESDestroy(&power->snes);
1065: }
1066: PetscFree(eps->data);
1067: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL);
1068: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL);
1069: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL);
1070: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL);
1071: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL);
1072: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL);
1073: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL);
1074: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL);
1075: return(0);
1076: }
1078: PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1079: {
1081: EPS_POWER *power = (EPS_POWER*)eps->data;
1082: PetscBool isascii;
1085: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1086: if (isascii) {
1087: if (power->nonlinear) {
1088: PetscViewerASCIIPrintf(viewer," using nonlinear inverse iteration\n");
1089: if (power->update) {
1090: PetscViewerASCIIPrintf(viewer," updating the residual monolithically\n");
1091: }
1092: if (!power->snes) { EPSPowerGetSNES(eps,&power->snes); }
1093: PetscViewerASCIIPushTab(viewer);
1094: SNESView(power->snes,viewer);
1095: PetscViewerASCIIPopTab(viewer);
1096: } else {
1097: PetscViewerASCIIPrintf(viewer," %s shifts\n",EPSPowerShiftTypes[power->shift_type]);
1098: }
1099: }
1100: return(0);
1101: }
1103: PetscErrorCode EPSComputeVectors_Power(EPS eps)
1104: {
1106: EPS_POWER *power = (EPS_POWER*)eps->data;
1107: PetscReal norm;
1108: PetscInt i;
1111: if (eps->twosided) {
1112: EPSComputeVectors_Twosided(eps);
1113: /* normalize (no need to take care of 2x2 blocks */
1114: for (i=0;i<eps->nconv;i++) {
1115: BVNormColumn(eps->V,i,NORM_2,&norm);
1116: BVScaleColumn(eps->V,i,1.0/norm);
1117: BVNormColumn(eps->W,i,NORM_2,&norm);
1118: BVScaleColumn(eps->W,i,1.0/norm);
1119: }
1120: } else if (!power->nonlinear) {
1121: EPSComputeVectors_Schur(eps);
1122: }
1123: return(0);
1124: }
1126: PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1127: {
1129: EPS_POWER *power = (EPS_POWER*)eps->data;
1130: KSP ksp;
1131: PC pc;
1134: if (power->nonlinear) {
1135: STGetKSP(eps->st,&ksp);
1136: KSPGetPC(ksp,&pc);
1137: PCSetType(pc,PCNONE);
1138: }
1139: return(0);
1140: }
1142: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1143: {
1144: EPS_POWER *ctx;
1148: PetscNewLog(eps,&ctx);
1149: eps->data = (void*)ctx;
1151: eps->useds = PETSC_TRUE;
1152: eps->hasts = PETSC_TRUE;
1153: eps->categ = EPS_CATEGORY_OTHER;
1155: eps->ops->setup = EPSSetUp_Power;
1156: eps->ops->setfromoptions = EPSSetFromOptions_Power;
1157: eps->ops->reset = EPSReset_Power;
1158: eps->ops->destroy = EPSDestroy_Power;
1159: eps->ops->view = EPSView_Power;
1160: eps->ops->backtransform = EPSBackTransform_Power;
1161: eps->ops->computevectors = EPSComputeVectors_Power;
1162: eps->ops->setdefaultst = EPSSetDefaultST_Power;
1163: eps->stopping = EPSStopping_Power;
1165: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power);
1166: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power);
1167: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power);
1168: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power);
1169: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power);
1170: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power);
1171: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power);
1172: PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power);
1173: return(0);
1174: }