Actual source code: power.c

slepc-3.14.2 2021-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc eigensolver: "power"

 13:    Method: Power Iteration

 15:    Algorithm:

 17:        This solver implements the power iteration for finding dominant
 18:        eigenpairs. It also includes the following well-known methods:
 19:        - Inverse Iteration: when used in combination with shift-and-invert
 20:          spectral transformation.
 21:        - Rayleigh Quotient Iteration (RQI): also with shift-and-invert plus
 22:          a variable shift.

 24:        It can also be used for nonlinear inverse iteration on the problem
 25:        A(x)*x=lambda*B(x)*x, where A and B are not constant but depend on x.

 27:    References:

 29:        [1] "Single Vector Iteration Methods in SLEPc", SLEPc Technical Report
 30:            STR-2, available at https://slepc.upv.es.
 31: */

 33: #include <slepc/private/epsimpl.h>
 34: #include <slepcblaslapack.h>
 35: /* petsc headers */
 36: #include <petscdm.h>
 37: #include <petscsnes.h>

 39: static PetscErrorCode EPSPowerFormFunction_Update(SNES,Vec,Vec,void*);
 40: PetscErrorCode EPSSolve_Power(EPS);
 41: PetscErrorCode EPSSolve_TS_Power(EPS);

 43: typedef struct {
 44:   EPSPowerShiftType shift_type;
 45:   PetscBool         nonlinear;
 46:   PetscBool         update;
 47:   SNES              snes;
 48:   PetscErrorCode    (*formFunctionB)(SNES,Vec,Vec,void*);
 49:   void              *formFunctionBctx;
 50:   PetscErrorCode    (*formFunctionA)(SNES,Vec,Vec,void*);
 51:   void              *formFunctionActx;
 52:   PetscErrorCode    (*formFunctionAB)(SNES,Vec,Vec,Vec,void*);
 53:   PetscInt          idx;  /* index of the first nonzero entry in the iteration vector */
 54:   PetscMPIInt       p;    /* process id of the owner of idx */
 55:   PetscReal         norm0; /* norm of initial vector */
 56: } EPS_POWER;

 58: static PetscErrorCode SNESMonitor_PowerUpdate(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
 59: {
 61:   EPS            eps;

 64:   PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
 65:   if (!eps) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
 66:   /* Call EPS monitor on each SNES iteration */
 67:   EPSMonitor(eps,its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
 68:   return(0);
 69: }

 71: PetscErrorCode EPSSetUp_Power(EPS eps)
 72: {
 74:   EPS_POWER      *power = (EPS_POWER*)eps->data;
 75:   STMatMode      mode;
 76:   Mat            A,B,P;
 77:   Vec            res;
 78:   PetscContainer container;
 79:   PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
 80:   PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
 81:   void           *ctx;

 84:   if (eps->ncv!=PETSC_DEFAULT) {
 85:     if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
 86:   } else eps->ncv = eps->nev;
 87:   if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
 88:   if (eps->max_it==PETSC_DEFAULT) {
 89:     /* SNES will directly return the solution for us, and we need to do only one iteration */
 90:     if (power->nonlinear && power->update) eps->max_it = 1;
 91:     else eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
 92:   }
 93:   if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
 94:   if (eps->which!=EPS_LARGEST_MAGNITUDE && eps->which!=EPS_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
 95:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
 96:     if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
 97:     EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," (with variable shifts)");
 98:     STGetMatMode(eps->st,&mode);
 99:     if (mode == ST_MATMODE_INPLACE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
100:   }
101:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE);
102:   EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
103:   EPSAllocateSolution(eps,0);
104:   EPS_SetInnerProduct(eps);

106:   if (power->nonlinear) {
107:     if (eps->nev>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
108:     EPSSetWorkVecs(eps,3);
109:     if (power->update && eps->max_it>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"More than one iteration is not allowed for Newton eigensolver (SNES)");

111:     /* Set up SNES solver */
112:     /* If SNES was setup, we need to reset it so that it will be consistent with the current EPS */
113:     if (power->snes) { SNESReset(power->snes); }
114:     else { EPSPowerGetSNES(eps,&power->snes); }

116:     /* For nonlinear solver (Newton), we should scale the initial vector back.
117:        The initial vector will be scaled in EPSSetUp. */
118:     if (eps->IS) {
119:       VecNorm((eps->IS)[0],NORM_2,&power->norm0);
120:     }

122:     EPSGetOperators(eps,&A,&B);

124:     PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA);
125:     if (!formFunctionA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
126:     PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container);
127:     if (container) {
128:       PetscContainerGetPointer(container,&ctx);
129:     } else ctx = NULL;
130:     MatCreateVecs(A,&res,NULL);
131:     power->formFunctionA = formFunctionA;
132:     power->formFunctionActx = ctx;
133:     if (power->update) {
134:       SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx);
135:       PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB);
136:       SNESMonitorSet(power->snes,SNESMonitor_PowerUpdate,NULL,NULL);
137:     }
138:     else { SNESSetFunction(power->snes,res,formFunctionA,ctx); }
139:     VecDestroy(&res);

141:     PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA);
142:     if (!formJacobianA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
143:     PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container);
144:     if (container) {
145:       PetscContainerGetPointer(container,&ctx);
146:     } else ctx = NULL;
147:     /* This allows users to compute a different preconditioning matrix than A.
148:      * It is useful when A and B are shell matrices.
149:      */
150:     STPrecondGetMatForPC(eps->st,&P);
151:     /* If users do not provide a matrix, we simply use A */
152:     SNESSetJacobian(power->snes,A,P? P:A,formJacobianA,ctx);
153:     SNESSetFromOptions(power->snes);
154:     SNESSetUp(power->snes);
155:     if (B) {
156:       PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB);
157:       PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container);
158:       if (power->formFunctionB && container) {
159:         PetscContainerGetPointer(container,&power->formFunctionBctx);
160:       } else power->formFunctionBctx = NULL;
161:     }
162:   } else {
163:     if (eps->twosided) {
164:       EPSSetWorkVecs(eps,3);
165:     } else {
166:       EPSSetWorkVecs(eps,2);
167:     }
168:     DSSetType(eps->ds,DSNHEP);
169:     DSAllocate(eps->ds,eps->nev);
170:   }
171:   /* dispatch solve method */
172:   if (eps->twosided) {
173:     if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
174:     if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
175:     eps->ops->solve = EPSSolve_TS_Power;
176:   } else eps->ops->solve = EPSSolve_Power;
177:   return(0);
178: }

180: /*
181:    Find the index of the first nonzero entry of x, and the process that owns it.
182: */
183: static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscMPIInt *p)
184: {
185:   PetscErrorCode    ierr;
186:   PetscInt          i,first,last,N;
187:   PetscLayout       map;
188:   const PetscScalar *xx;

191:   VecGetSize(x,&N);
192:   VecGetOwnershipRange(x,&first,&last);
193:   VecGetArrayRead(x,&xx);
194:   for (i=first;i<last;i++) {
195:     if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
196:   }
197:   VecRestoreArrayRead(x,&xx);
198:   if (i==last) i=N;
199:   MPIU_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x));
200:   if (*idx==N) SETERRQ(PetscObjectComm((PetscObject)x),1,"Zero vector found");
201:   VecGetLayout(x,&map);
202:   PetscLayoutFindOwner(map,*idx,p);
203:   return(0);
204: }

206: /*
207:    Normalize a vector x with respect to a given norm as well as the
208:    sign of the first nonzero entry (assumed to be idx in process p).
209: */
210: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscMPIInt p,PetscScalar *sign)
211: {
212:   PetscErrorCode    ierr;
213:   PetscScalar       alpha=1.0;
214:   PetscInt          first,last;
215:   PetscReal         absv;
216:   const PetscScalar *xx;

219:   VecGetOwnershipRange(x,&first,&last);
220:   if (idx>=first && idx<last) {
221:     VecGetArrayRead(x,&xx);
222:     absv = PetscAbsScalar(xx[idx-first]);
223:     if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
224:     VecRestoreArrayRead(x,&xx);
225:   }
226:   MPI_Bcast(&alpha,1,MPIU_SCALAR,p,PetscObjectComm((PetscObject)x));
227:   if (sign) *sign = alpha;
228:   alpha *= norm;
229:   VecScale(x,1.0/alpha);
230:   return(0);
231: }

233: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
234: {
236:   EPS_POWER      *power = (EPS_POWER*)eps->data;
237:   Mat            B;

240:   STResetMatrixState(eps->st);
241:   EPSGetOperators(eps,NULL,&B);
242:   if (B) {
243:     if (power->formFunctionB) {
244:       (*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx);
245:     } else {
246:       MatMult(B,x,Bx);
247:     }
248:   } else {
249:     VecCopy(x,Bx);
250:   }
251:   return(0);
252: }

254: static PetscErrorCode EPSPowerUpdateFunctionA(EPS eps,Vec x,Vec Ax)
255: {
257:   EPS_POWER      *power = (EPS_POWER*)eps->data;
258:   Mat            A;

261:   STResetMatrixState(eps->st);
262:   EPSGetOperators(eps,&A,NULL);
263:   if (A) {
264:     if (power->formFunctionA) {
265:       (*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx);
266:     } else {
267:       MatMult(A,x,Ax);
268:     }
269:   } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
270:   return(0);
271: }

273: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
274: {
276:   EPS            eps;
277:   PetscReal      bx;
278:   Vec            Bx;
279:   PetscScalar    sign;
280:   EPS_POWER      *power;

283:   PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
284:   if (!eps) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
285:   STResetMatrixState(eps->st);
286:   power = (EPS_POWER*)eps->data;
287:   Bx = eps->work[2];
288:   if (power->formFunctionAB) {
289:     (*power->formFunctionAB)(snes,x,y,Bx,ctx);
290:   } else {
291:     EPSPowerUpdateFunctionA(eps,x,y);
292:     EPSPowerUpdateFunctionB(eps,x,Bx);
293:   }
294:   VecNorm(Bx,NORM_2,&bx);
295:   Normalize(Bx,bx,power->idx,power->p,&sign);
296:   VecAXPY(y,-1.0,Bx);
297:   /* Keep tracking eigenvalue update. It would be useful when we want to monitor solver progress via snes monitor. */
298:   eps->eigr[(eps->nconv < eps->nev)? eps->nconv:(eps->nconv-1)] = 1.0/(bx*sign);
299:   return(0);
300: }

302: /*
303:    Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
304: */
305: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
306: {
308:   EPS_POWER      *power = (EPS_POWER*)eps->data;
309:   Vec            Bx;

312:   VecCopy(x,y);
313:   if (power->update) {
314:     SNESSolve(power->snes,NULL,y);
315:   } else {
316:     Bx = eps->work[2];
317:     SNESSolve(power->snes,Bx,y);
318:   }
319:   return(0);
320: }

322: /*
323:    Use nonlinear inverse power to compute an initial guess
324: */
325: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
326: {
327:   EPS            powereps;
328:   Mat            A,B,P;
329:   Vec            v1,v2;
330:   SNES           snes;
331:   DM             dm,newdm;

335:   EPSCreate(PetscObjectComm((PetscObject)eps),&powereps);
336:   EPSGetOperators(eps,&A,&B);
337:   EPSSetType(powereps,EPSPOWER);
338:   EPSSetOperators(powereps,A,B);
339:   EPSSetTolerances(powereps,1e-6,4);
340:   EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix);
341:   EPSAppendOptionsPrefix(powereps,"init_");
342:   EPSSetProblemType(powereps,EPS_GNHEP);
343:   EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE);
344:   EPSPowerSetNonlinear(powereps,PETSC_TRUE);
345:   STPrecondGetMatForPC(eps->st,&P);
346:   /* attach dm to initial solve */
347:   EPSPowerGetSNES(eps,&snes);
348:   SNESGetDM(snes,&dm);
349:   /* use  dmshell to temporarily store snes context */
350:   DMCreate(PetscObjectComm((PetscObject)eps),&newdm);
351:   DMSetType(newdm,DMSHELL);
352:   DMSetUp(newdm);
353:   DMCopyDMSNES(dm,newdm);
354:   EPSPowerGetSNES(powereps,&snes);
355:   SNESSetDM(snes,dm);
356:   EPSSetFromOptions(powereps);
357:   if (P) {
358:     STPrecondSetMatForPC(powereps->st,P);
359:   }
360:   EPSSolve(powereps);
361:   BVGetColumn(eps->V,0,&v2);
362:   BVGetColumn(powereps->V,0,&v1);
363:   VecCopy(v1,v2);
364:   BVRestoreColumn(powereps->V,0,&v1);
365:   BVRestoreColumn(eps->V,0,&v2);
366:   EPSDestroy(&powereps);
367:   /* restore context back to the old nonlinear solver */
368:   DMCopyDMSNES(newdm,dm);
369:   DMDestroy(&newdm);
370:   return(0);
371: }

373: PetscErrorCode EPSSolve_Power(EPS eps)
374: {
375:   PetscErrorCode      ierr;
376:   EPS_POWER           *power = (EPS_POWER*)eps->data;
377:   PetscInt            k,ld;
378:   Vec                 v,y,e,Bx;
379:   Mat                 A;
380:   KSP                 ksp;
381:   PetscReal           relerr,norm,norm1,rt1,rt2,cs1;
382:   PetscScalar         theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
383:   PetscBool           breakdown;
384:   KSPConvergedReason  reason;
385:   SNESConvergedReason snesreason;

388:   e = eps->work[0];
389:   y = eps->work[1];
390:   if (power->nonlinear) Bx = eps->work[2];
391:   else Bx = NULL;

393:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
394:   if (power->nonlinear) {
395:     PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps);
396:     /* Compute an intial guess only when users do not provide one */
397:     if (power->update && !eps->nini) {
398:       EPSPowerComputeInitialGuess_Update(eps);
399:     }
400:   } else {
401:     DSGetLeadingDimension(eps->ds,&ld);
402:   }
403:   if (!power->update) {
404:     EPSGetStartVector(eps,0,NULL);
405:   }
406:   if (power->nonlinear) {
407:     BVGetColumn(eps->V,0,&v);
408:     if (eps->nini) {
409:       /* We scale the initial vector back if the initial vector was provided by users */
410:       VecScale(v,power->norm0);
411:     }
412:     EPSPowerUpdateFunctionB(eps,v,Bx);
413:     VecNorm(Bx,NORM_2,&norm);
414:     FirstNonzeroIdx(Bx,&power->idx,&power->p);
415:     Normalize(Bx,norm,power->idx,power->p,NULL);
416:     BVRestoreColumn(eps->V,0,&v);
417:   }

419:   STGetShift(eps->st,&sigma);    /* original shift */
420:   rho = sigma;

422:   while (eps->reason == EPS_CONVERGED_ITERATING) {
423:     eps->its++;
424:     k = eps->nconv;

426:     /* y = OP v */
427:     BVGetColumn(eps->V,k,&v);
428:     if (power->nonlinear) {
429:       EPSPowerApply_SNES(eps,v,y);
430:     } else {
431:       STApply(eps->st,v,y);
432:     }
433:     BVRestoreColumn(eps->V,k,&v);

435:     /* purge previously converged eigenvectors */
436:     if (power->nonlinear) {
437:       /* We do not need to call this for Newton eigensolver since eigenvalue is
438:        * updated in function evaluations.
439:        */
440:       if (!power->update) {
441:         EPSPowerUpdateFunctionB(eps,y,Bx);
442:         VecNorm(Bx,NORM_2,&norm);
443:         Normalize(Bx,norm,power->idx,power->p,&sign);
444:       }
445:     } else {
446:       DSGetArray(eps->ds,DS_MAT_A,&T);
447:       BVSetActiveColumns(eps->V,0,k);
448:       BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL);
449:     }

451:     /* theta = (v,y)_B */
452:     BVSetActiveColumns(eps->V,k,k+1);
453:     BVDotVec(eps->V,y,&theta);
454:     if (!power->nonlinear) {
455:       T[k+k*ld] = theta;
456:       DSRestoreArray(eps->ds,DS_MAT_A,&T);
457:     }

459:     /* Eigenvalue is already stored in function evaluations.
460:      * Assign eigenvalue to theta to make the rest of the code consistent
461:      */
462:     if (power->update) theta = eps->eigr[eps->nconv];
463:     else if (power->nonlinear) theta = 1.0/norm*sign; /* Eigenvalue: 1/|Bx| */

465:     if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */

467:       /* approximate eigenvalue is the Rayleigh quotient */
468:       eps->eigr[eps->nconv] = theta;

470:       /**
471:        * If the Newton method (update, SNES) is used, we do not compute "relerr"
472:        * since SNES determines the convergence.
473:        */
474:       if (power->update) relerr = 0.;
475:       else {
476:         /* compute relative error as ||y-theta v||_2/|theta| */
477:         VecCopy(y,e);
478:         BVGetColumn(eps->V,k,&v);
479:         VecAXPY(e,power->nonlinear?-1.0:-theta,v);
480:         BVRestoreColumn(eps->V,k,&v);
481:         VecNorm(e,NORM_2,&relerr);
482:         if (power->nonlinear) relerr *= PetscAbsScalar(theta);
483:         else relerr /= PetscAbsScalar(theta);
484:       }

486:     } else {  /* RQI */

488:       /* delta = ||y||_B */
489:       delta = norm;

491:       /* compute relative error */
492:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
493:       else relerr = 1.0 / (norm*PetscAbsScalar(rho));

495:       /* approximate eigenvalue is the shift */
496:       eps->eigr[eps->nconv] = rho;

498:       /* compute new shift */
499:       if (relerr<eps->tol) {
500:         rho = sigma;  /* if converged, restore original shift */
501:         STSetShift(eps->st,rho);
502:       } else {
503:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
504:         if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
505:           /* beta1 is the norm of the residual associated with R(v) */
506:           BVGetColumn(eps->V,k,&v);
507:           VecAXPY(v,-PetscConj(theta)/(delta*delta),y);
508:           BVRestoreColumn(eps->V,k,&v);
509:           BVScaleColumn(eps->V,k,1.0/delta);
510:           BVNormColumn(eps->V,k,NORM_2,&norm1);
511:           beta1 = norm1;

513:           /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
514:           STGetMatrix(eps->st,0,&A);
515:           BVGetColumn(eps->V,k,&v);
516:           MatMult(A,v,e);
517:           VecDot(v,e,&alpha2);
518:           BVRestoreColumn(eps->V,k,&v);
519:           alpha2 = alpha2 / (beta1 * beta1);

521:           /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
522:           PetscFPTrapPush(PETSC_FP_TRAP_OFF);
523:           PetscStackCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
524:           PetscFPTrapPop();
525:           if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
526:           else rho = rt2;
527:         }
528:         /* update operator according to new shift */
529:         KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
530:         STSetShift(eps->st,rho);
531:         KSPGetConvergedReason(ksp,&reason);
532:         if (reason) {
533:           PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
534:           rho *= 1+10*PETSC_MACHINE_EPSILON;
535:           STSetShift(eps->st,rho);
536:           KSPGetConvergedReason(ksp,&reason);
537:           if (reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
538:         }
539:         KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
540:       }
541:     }
542:     eps->errest[eps->nconv] = relerr;

544:     /* normalize */
545:     if (!power->nonlinear) { Normalize(y,norm,power->idx,power->p,NULL); }
546:     BVInsertVec(eps->V,k,y);

548:     if (power->update) {
549:       SNESGetConvergedReason(power->snes,&snesreason);
550:       /* For Newton eigensolver, we are ready to return once SNES converged. */
551:       if (snesreason>0) eps->nconv = 1;
552:     } else if (relerr<eps->tol) {   /* accept eigenpair */
553:       eps->nconv = eps->nconv + 1;
554:       if (eps->nconv<eps->nev) {
555:         EPSGetStartVector(eps,eps->nconv,&breakdown);
556:         if (breakdown) {
557:           eps->reason = EPS_DIVERGED_BREAKDOWN;
558:           PetscInfo(eps,"Unable to generate more start vectors\n");
559:           break;
560:         }
561:       }
562:     }
563:     /* For Newton eigensolver, monitor will be called from SNES monitor */
564:     if (!power->update) {
565:       EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
566:     }
567:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);

569:     /**
570:      * When a customized stopping test is used, and reason can be set to be converged (EPS_CONVERGED_USER).
571:      * In this case, we need to increase eps->nconv to "1" so users can retrieve the solution.
572:      */
573:     if (power->nonlinear && eps->reason>0) eps->nconv = 1;
574:   }

576:   if (power->nonlinear) {
577:     PetscObjectCompose((PetscObject)power->snes,"eps",NULL);
578:   } else {
579:     DSSetDimensions(eps->ds,eps->nconv,0,0,0);
580:     DSSetState(eps->ds,DS_STATE_RAW);
581:   }
582:   return(0);
583: }

585: PetscErrorCode EPSSolve_TS_Power(EPS eps)
586: {
587:   PetscErrorCode     ierr;
588:   EPS_POWER          *power = (EPS_POWER*)eps->data;
589:   PetscInt           k,ld;
590:   Vec                v,w,y,e,z;
591:   KSP                ksp;
592:   PetscReal          relerr=1.0,relerrl,delta;
593:   PetscScalar        theta,rho,alpha,sigma;
594:   PetscBool          breakdown,breakdownl;
595:   KSPConvergedReason reason;

598:   e = eps->work[0];
599:   v = eps->work[1];
600:   w = eps->work[2];

602:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
603:   DSGetLeadingDimension(eps->ds,&ld);
604:   EPSGetStartVector(eps,0,NULL);
605:   EPSGetLeftStartVector(eps,0,NULL);
606:   BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL);
607:   BVCopyVec(eps->V,0,v);
608:   BVCopyVec(eps->W,0,w);
609:   STGetShift(eps->st,&sigma);    /* original shift */
610:   rho = sigma;

612:   while (eps->reason == EPS_CONVERGED_ITERATING) {
613:     eps->its++;
614:     k = eps->nconv;

616:     /* y = OP v, z = OP' w */
617:     BVGetColumn(eps->V,k,&y);
618:     STApply(eps->st,v,y);
619:     BVRestoreColumn(eps->V,k,&y);
620:     BVGetColumn(eps->W,k,&z);
621:     STApplyHermitianTranspose(eps->st,w,z);
622:     BVRestoreColumn(eps->W,k,&z);

624:     /* purge previously converged eigenvectors */
625:     BVBiorthogonalizeColumn(eps->V,eps->W,k);

627:     /* theta = (w,y)_B */
628:     BVSetActiveColumns(eps->V,k,k+1);
629:     BVDotVec(eps->V,w,&theta);
630:     theta = PetscConj(theta);

632:     if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) { /* direct & inverse iteration */

634:       /* approximate eigenvalue is the Rayleigh quotient */
635:       eps->eigr[eps->nconv] = theta;

637:       /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
638:       BVCopyVec(eps->V,k,e);
639:       VecAXPY(e,-theta,v);
640:       VecNorm(e,NORM_2,&relerr);
641:       BVCopyVec(eps->W,k,e);
642:       VecAXPY(e,-PetscConj(theta),w);
643:       VecNorm(e,NORM_2,&relerrl);
644:       relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
645:     }

647:     /* normalize */
648:     BVSetActiveColumns(eps->V,k,k+1);
649:     BVGetColumn(eps->W,k,&z);
650:     BVDotVec(eps->V,z,&alpha);
651:     BVRestoreColumn(eps->W,k,&z);
652:     delta = PetscSqrtReal(PetscAbsScalar(alpha));
653:     if (delta==0.0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Breakdown in two-sided Power/RQI");
654:     BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta));
655:     BVScaleColumn(eps->W,k,1.0/delta);
656:     BVCopyVec(eps->V,k,v);
657:     BVCopyVec(eps->W,k,w);

659:     if (power->shift_type == EPS_POWER_SHIFT_RAYLEIGH) { /* RQI */

661:       /* compute relative error */
662:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
663:       else relerr = 1.0 / PetscAbsScalar(delta*rho);

665:       /* approximate eigenvalue is the shift */
666:       eps->eigr[eps->nconv] = rho;

668:       /* compute new shift */
669:       if (relerr<eps->tol) {
670:         rho = sigma;  /* if converged, restore original shift */
671:         STSetShift(eps->st,rho);
672:       } else {
673:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
674:         /* update operator according to new shift */
675:         KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
676:         STSetShift(eps->st,rho);
677:         KSPGetConvergedReason(ksp,&reason);
678:         if (reason) {
679:           PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
680:           rho *= 1+10*PETSC_MACHINE_EPSILON;
681:           STSetShift(eps->st,rho);
682:           KSPGetConvergedReason(ksp,&reason);
683:           if (reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
684:         }
685:         KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
686:       }
687:     }
688:     eps->errest[eps->nconv] = relerr;

690:     /* if relerr<tol, accept eigenpair */
691:     if (relerr<eps->tol) {
692:       eps->nconv = eps->nconv + 1;
693:       if (eps->nconv<eps->nev) {
694:         EPSGetStartVector(eps,eps->nconv,&breakdown);
695:         EPSGetLeftStartVector(eps,eps->nconv,&breakdownl);
696:         if (breakdown || breakdownl) {
697:           eps->reason = EPS_DIVERGED_BREAKDOWN;
698:           PetscInfo(eps,"Unable to generate more start vectors\n");
699:           break;
700:         }
701:         BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL);
702:       }
703:     }
704:     EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
705:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
706:   }

708:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
709:   DSSetState(eps->ds,DS_STATE_RAW);
710:   return(0);
711: }

713: PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
714: {
716:   EPS_POWER      *power = (EPS_POWER*)eps->data;
717:   SNESConvergedReason snesreason;

720:   if (power->update) {
721:     SNESGetConvergedReason(power->snes,&snesreason);
722:     if (snesreason < 0) {
723:       *reason = EPS_DIVERGED_BREAKDOWN;
724:       return(0);
725:     }
726:   }
727:   EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx);
728:   return(0);
729: }

731: PetscErrorCode EPSBackTransform_Power(EPS eps)
732: {
734:   EPS_POWER      *power = (EPS_POWER*)eps->data;

737:   if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
738:   else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) {
739:     EPSBackTransform_Default(eps);
740:   }
741:   return(0);
742: }

744: PetscErrorCode EPSSetFromOptions_Power(PetscOptionItems *PetscOptionsObject,EPS eps)
745: {
746:   PetscErrorCode    ierr;
747:   EPS_POWER         *power = (EPS_POWER*)eps->data;
748:   PetscBool         flg,val;
749:   EPSPowerShiftType shift;

752:   PetscOptionsHead(PetscOptionsObject,"EPS Power Options");

754:     PetscOptionsEnum("-eps_power_shift_type","Shift type","EPSPowerSetShiftType",EPSPowerShiftTypes,(PetscEnum)power->shift_type,(PetscEnum*)&shift,&flg);
755:     if (flg) { EPSPowerSetShiftType(eps,shift); }

757:     PetscOptionsBool("-eps_power_nonlinear","Use nonlinear inverse iteration","EPSPowerSetNonlinear",power->nonlinear,&val,&flg);
758:     if (flg) { EPSPowerSetNonlinear(eps,val); }

760:     PetscOptionsBool("-eps_power_update","Update residual monolithically","EPSPowerSetUpdate",power->update,&val,&flg);
761:     if (flg) { EPSPowerSetUpdate(eps,val); }

763:   PetscOptionsTail();
764:   return(0);
765: }

767: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
768: {
769:   EPS_POWER *power = (EPS_POWER*)eps->data;

772:   switch (shift) {
773:     case EPS_POWER_SHIFT_CONSTANT:
774:     case EPS_POWER_SHIFT_RAYLEIGH:
775:     case EPS_POWER_SHIFT_WILKINSON:
776:       if (power->shift_type != shift) {
777:         power->shift_type = shift;
778:         eps->state = EPS_STATE_INITIAL;
779:       }
780:       break;
781:     default:
782:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
783:   }
784:   return(0);
785: }

787: /*@
788:    EPSPowerSetShiftType - Sets the type of shifts used during the power
789:    iteration. This can be used to emulate the Rayleigh Quotient Iteration
790:    (RQI) method.

792:    Logically Collective on eps

794:    Input Parameters:
795: +  eps - the eigenproblem solver context
796: -  shift - the type of shift

798:    Options Database Key:
799: .  -eps_power_shift_type - Sets the shift type (either 'constant' or
800:                            'rayleigh' or 'wilkinson')

802:    Notes:
803:    By default, shifts are constant (EPS_POWER_SHIFT_CONSTANT) and the iteration
804:    is the simple power method (or inverse iteration if a shift-and-invert
805:    transformation is being used).

807:    A variable shift can be specified (EPS_POWER_SHIFT_RAYLEIGH or
808:    EPS_POWER_SHIFT_WILKINSON). In this case, the iteration behaves rather like
809:    a cubic converging method such as RQI.

811:    Level: advanced

813: .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
814: @*/
815: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
816: {

822:   PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
823:   return(0);
824: }

826: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
827: {
828:   EPS_POWER *power = (EPS_POWER*)eps->data;

831:   *shift = power->shift_type;
832:   return(0);
833: }

835: /*@
836:    EPSPowerGetShiftType - Gets the type of shifts used during the power
837:    iteration.

839:    Not Collective

841:    Input Parameter:
842: .  eps - the eigenproblem solver context

844:    Input Parameter:
845: .  shift - the type of shift

847:    Level: advanced

849: .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
850: @*/
851: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
852: {

858:   PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
859:   return(0);
860: }

862: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
863: {
864:   EPS_POWER *power = (EPS_POWER*)eps->data;

867:   if (power->nonlinear != nonlinear) {
868:     power->nonlinear = nonlinear;
869:     eps->useds = PetscNot(nonlinear);
870:     eps->ops->setupsort = nonlinear? NULL: EPSSetUpSort_Default;
871:     eps->state = EPS_STATE_INITIAL;
872:   }
873:   return(0);
874: }

876: /*@
877:    EPSPowerSetNonlinear - Sets a flag to indicate that the problem is nonlinear.

879:    Logically Collective on eps

881:    Input Parameters:
882: +  eps - the eigenproblem solver context
883: -  nonlinear - whether the problem is nonlinear or not

885:    Options Database Key:
886: .  -eps_power_nonlinear - Sets the nonlinear flag

888:    Notes:
889:    If this flag is set, the solver assumes that the problem is nonlinear,
890:    that is, the operators that define the eigenproblem are not constant
891:    matrices, but depend on the eigenvector: A(x)*x=lambda*B(x)*x. This is
892:    different from the case of nonlinearity with respect to the eigenvalue
893:    (use the NEP solver class for this kind of problems).

895:    The way in which nonlinear operators are specified is very similar to
896:    the case of PETSc's SNES solver. The difference is that the callback
897:    functions are provided via composed functions "formFunction" and
898:    "formJacobian" in each of the matrix objects passed as arguments of
899:    EPSSetOperators(). The application context required for these functions
900:    can be attached via a composed PetscContainer.

902:    Level: advanced

904: .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
905: @*/
906: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
907: {

913:   PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
914:   return(0);
915: }

917: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
918: {
919:   EPS_POWER *power = (EPS_POWER*)eps->data;

922:   *nonlinear = power->nonlinear;
923:   return(0);
924: }

926: /*@
927:    EPSPowerGetNonlinear - Returns a flag indicating if the problem is nonlinear.

929:    Not Collective

931:    Input Parameter:
932: .  eps - the eigenproblem solver context

934:    Input Parameter:
935: .  nonlinear - the nonlinear flag

937:    Level: advanced

939: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
940: @*/
941: PetscErrorCode EPSPowerGetNonlinear(EPS eps,PetscBool *nonlinear)
942: {

948:   PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
949:   return(0);
950: }

952: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
953: {
954:   EPS_POWER *power = (EPS_POWER*)eps->data;

957:   if (!power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
958:   power->update = update;
959:   eps->state = EPS_STATE_INITIAL;
960:   return(0);
961: }

963: /*@
964:    EPSPowerSetUpdate - Sets a flag to indicate that the residual is updated monolithically
965:    for nonlinear problems. This potentially has a better convergence.

967:    Logically Collective on eps

969:    Input Parameters:
970: +  eps - the eigenproblem solver context
971: -  update - whether the residual is updated monolithically or not

973:    Options Database Key:
974: .  -eps_power_update - Sets the update flag

976:    Level: advanced

978: .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
979: @*/
980: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
981: {

987:   PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
988:   return(0);
989: }

991: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
992: {
993:   EPS_POWER *power = (EPS_POWER*)eps->data;

996:   *update = power->update;
997:   return(0);
998: }

1000: /*@
1001:    EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
1002:    for nonlinear problems.

1004:    Not Collective

1006:    Input Parameter:
1007: .  eps - the eigenproblem solver context

1009:    Input Parameter:
1010: .  update - the update flag

1012:    Level: advanced

1014: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
1015: @*/
1016: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
1017: {

1023:   PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
1024:   return(0);
1025: }

1027: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
1028: {
1030:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1033:   PetscObjectReference((PetscObject)snes);
1034:   SNESDestroy(&power->snes);
1035:   power->snes = snes;
1036:   PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
1037:   eps->state = EPS_STATE_INITIAL;
1038:   return(0);
1039: }

1041: /*@
1042:    EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
1043:    eigenvalue solver (to be used in nonlinear inverse iteration).

1045:    Collective on eps

1047:    Input Parameters:
1048: +  eps  - the eigenvalue solver
1049: -  snes - the nonlinear solver object

1051:    Level: advanced

1053: .seealso: EPSPowerGetSNES()
1054: @*/
1055: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1056: {

1063:   PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1064:   return(0);
1065: }

1067: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1068: {
1070:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1073:   if (!power->snes) {
1074:     SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes);
1075:     PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1);
1076:     SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix);
1077:     SNESAppendOptionsPrefix(power->snes,"eps_power_");
1078:     PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
1079:     PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options);
1080:   }
1081:   *snes = power->snes;
1082:   return(0);
1083: }

1085: /*@
1086:    EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1087:    with the eigenvalue solver.

1089:    Not Collective

1091:    Input Parameter:
1092: .  eps - the eigenvalue solver

1094:    Output Parameter:
1095: .  snes - the nonlinear solver object

1097:    Level: advanced

1099: .seealso: EPSPowerSetSNES()
1100: @*/
1101: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1102: {

1108:   PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1109:   return(0);
1110: }

1112: PetscErrorCode EPSReset_Power(EPS eps)
1113: {
1115:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1118:   if (power->snes) { SNESReset(power->snes); }
1119:   return(0);
1120: }

1122: PetscErrorCode EPSDestroy_Power(EPS eps)
1123: {
1125:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1128:   if (power->nonlinear) {
1129:     SNESDestroy(&power->snes);
1130:   }
1131:   PetscFree(eps->data);
1132:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL);
1133:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL);
1134:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL);
1135:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL);
1136:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL);
1137:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL);
1138:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL);
1139:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL);
1140:   return(0);
1141: }

1143: PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1144: {
1146:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1147:   PetscBool      isascii;

1150:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1151:   if (isascii) {
1152:     if (power->nonlinear) {
1153:       PetscViewerASCIIPrintf(viewer,"  using nonlinear inverse iteration\n");
1154:       if (power->update) {
1155:         PetscViewerASCIIPrintf(viewer,"  updating the residual monolithically\n");
1156:       }
1157:       if (!power->snes) { EPSPowerGetSNES(eps,&power->snes); }
1158:       PetscViewerASCIIPushTab(viewer);
1159:       SNESView(power->snes,viewer);
1160:       PetscViewerASCIIPopTab(viewer);
1161:     } else {
1162:       PetscViewerASCIIPrintf(viewer,"  %s shifts\n",EPSPowerShiftTypes[power->shift_type]);
1163:     }
1164:   }
1165:   return(0);
1166: }

1168: PetscErrorCode EPSComputeVectors_Power(EPS eps)
1169: {
1171:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1174:   if (eps->twosided) {
1175:     EPSComputeVectors_Twosided(eps);
1176:     BVNormalize(eps->V,NULL);
1177:     BVNormalize(eps->W,NULL);
1178:   } else if (!power->nonlinear) {
1179:     EPSComputeVectors_Schur(eps);
1180:   }
1181:   return(0);
1182: }

1184: PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1185: {
1187:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1188:   KSP            ksp;
1189:   PC             pc;

1192:   if (power->nonlinear) {
1193:     eps->categ=EPS_CATEGORY_PRECOND;
1194:     STGetKSP(eps->st,&ksp);
1195:     /* Set ST as STPRECOND so it can carry one preconditioning matrix
1196:      * It is useful when A and B are shell matrices
1197:      */
1198:     STSetType(eps->st,STPRECOND);
1199:     KSPGetPC(ksp,&pc);
1200:     PCSetType(pc,PCNONE);
1201:   }
1202:   return(0);
1203: }

1205: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1206: {
1207:   EPS_POWER      *ctx;

1211:   PetscNewLog(eps,&ctx);
1212:   eps->data = (void*)ctx;

1214:   eps->useds = PETSC_TRUE;
1215:   eps->categ = EPS_CATEGORY_OTHER;

1217:   eps->ops->setup          = EPSSetUp_Power;
1218:   eps->ops->setupsort      = EPSSetUpSort_Default;
1219:   eps->ops->setfromoptions = EPSSetFromOptions_Power;
1220:   eps->ops->reset          = EPSReset_Power;
1221:   eps->ops->destroy        = EPSDestroy_Power;
1222:   eps->ops->view           = EPSView_Power;
1223:   eps->ops->backtransform  = EPSBackTransform_Power;
1224:   eps->ops->computevectors = EPSComputeVectors_Power;
1225:   eps->ops->setdefaultst   = EPSSetDefaultST_Power;
1226:   eps->stopping            = EPSStopping_Power;

1228:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power);
1229:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power);
1230:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power);
1231:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power);
1232:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power);
1233:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power);
1234:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power);
1235:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power);
1236:   return(0);
1237: }