Actual source code: power.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }