Actual source code: power.c

slepc-3.12.2 2020-01-13
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:   PetscInt          idx;  /* index of the first nonzero entry in the iteration vector */
 55:   PetscInt          p;    /* process id of the owner of idx */
 56: } EPS_POWER;

 58: PetscErrorCode EPSSetUp_Power(EPS eps)
 59: {
 61:   EPS_POWER      *power = (EPS_POWER*)eps->data;
 62:   PetscBool      flg,istrivial;
 63:   STMatMode      mode;
 64:   Mat            A,B;
 65:   Vec            res;
 66:   PetscContainer container;
 67:   PetscErrorCode (*formFunctionA)(SNES,Vec,Vec,void*);
 68:   PetscErrorCode (*formJacobianA)(SNES,Vec,Mat,Mat,void*);
 69:   void           *ctx;
 70:   SNESLineSearch linesearch;

 73:   if (eps->ncv) {
 74:     if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
 75:   } else eps->ncv = eps->nev;
 76:   if (eps->mpd) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
 77:   if (!eps->max_it) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
 78:   if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
 79:   if (eps->which!=EPS_LARGEST_MAGNITUDE && eps->which !=EPS_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
 80:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) {
 81:     if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts not allowed in nonlinear problems");
 82:     PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STCAYLEY,"");
 83:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Variable shifts only allowed in shift-and-invert or Cayley ST");
 84:     STGetMatMode(eps->st,&mode);
 85:     if (mode == ST_MATMODE_INPLACE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"ST matrix mode inplace does not work with variable shifts");
 86:   }
 87:   if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
 88:   if (eps->balance!=EPS_BALANCE_NONE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Balancing not supported in this solver");
 89:   if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
 90:   RGIsTrivial(eps->rg,&istrivial);
 91:   if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");
 92:   EPSAllocateSolution(eps,0);
 93:   EPS_SetInnerProduct(eps);

 95:   if (power->nonlinear) {
 96:     if (eps->nev>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration cannot compute more than one eigenvalue");
 97:     EPSSetWorkVecs(eps,4);

 99:     /* set up SNES solver */
100:     if (!power->snes) { EPSPowerGetSNES(eps,&power->snes); }
101:     EPSGetOperators(eps,&A,&B);

103:     PetscObjectQueryFunction((PetscObject)A,"formFunction",&formFunctionA);
104:     if (!formFunctionA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formFunction' in matrix A");
105:     PetscObjectQuery((PetscObject)A,"formFunctionCtx",(PetscObject*)&container);
106:     if (container) {
107:       PetscContainerGetPointer(container,&ctx);
108:     } else ctx = NULL;
109:     MatCreateVecs(A,&res,NULL);
110:     power->formFunctionA = formFunctionA;
111:     power->formFunctionActx = ctx;
112:     if (power->update) {
113:       SNESSetFunction(power->snes,res,EPSPowerFormFunction_Update,ctx);
114:       PetscObjectQueryFunction((PetscObject)A,"formFunctionAB",&power->formFunctionAB);
115:     }
116:     else { SNESSetFunction(power->snes,res,formFunctionA,ctx); }
117:     VecDestroy(&res);

119:     PetscObjectQueryFunction((PetscObject)A,"formJacobian",&formJacobianA);
120:     if (!formJacobianA) SETERRQ(PetscObjectComm((PetscObject)eps),1,"For nonlinear inverse iteration you must compose a callback function 'formJacobian' in matrix A");
121:     PetscObjectQuery((PetscObject)A,"formJacobianCtx",(PetscObject*)&container);
122:     if (container) {
123:       PetscContainerGetPointer(container,&ctx);
124:     } else ctx = NULL;
125:     SNESSetJacobian(power->snes,A,A,formJacobianA,ctx);
126:     SNESSetFromOptions(power->snes);
127:     SNESGetLineSearch(power->snes,&linesearch);
128:     if (power->update) {
129:       SNESLineSearchSetPostCheck(linesearch,SNESLineSearchPostheckFunction,ctx);
130:     }
131:     SNESSetUp(power->snes);
132:     if (B) {
133:       PetscObjectQueryFunction((PetscObject)B,"formFunction",&power->formFunctionB);
134:       PetscObjectQuery((PetscObject)B,"formFunctionCtx",(PetscObject*)&container);
135:       if (power->formFunctionB && container) {
136:         PetscContainerGetPointer(container,&power->formFunctionBctx);
137:       } else power->formFunctionBctx = NULL;
138:     }
139:   } else {
140:     if (eps->twosided) {
141:       EPSSetWorkVecs(eps,3);
142:     } else {
143:       EPSSetWorkVecs(eps,2);
144:     }
145:     DSSetType(eps->ds,DSNHEP);
146:     DSAllocate(eps->ds,eps->nev);
147:   }
148:   /* dispatch solve method */
149:   if (eps->twosided) {
150:     if (power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Nonlinear inverse iteration does not have two-sided variant");
151:     if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided variant does not support Wilkinson shifts");
152:     eps->ops->solve = EPSSolve_TS_Power;
153:   } else eps->ops->solve = EPSSolve_Power;
154:   return(0);
155: }

157: /*
158:    Find the index of the first nonzero entry of x, and the process that owns it.
159: */
160: static PetscErrorCode FirstNonzeroIdx(Vec x,PetscInt *idx,PetscInt *p)
161: {
162:   PetscErrorCode    ierr;
163:   PetscInt          i,first,last,N;
164:   PetscLayout       map;
165:   const PetscScalar *xx;

168:   VecGetSize(x,&N);
169:   VecGetOwnershipRange(x,&first,&last);
170:   VecGetArrayRead(x,&xx);
171:   for (i=first;i<last;i++) {
172:     if (PetscAbsScalar(xx[i-first])>10*PETSC_MACHINE_EPSILON) break;
173:   }
174:   VecRestoreArrayRead(x,&xx);
175:   if (i==last) i=N;
176:   MPI_Allreduce(&i,idx,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)x));
177:   if (*idx==N) SETERRQ(PetscObjectComm((PetscObject)x),1,"Zero vector found");
178:   VecGetLayout(x,&map);
179:   PetscLayoutFindOwner(map,*idx,p);
180:   return(0);
181: }

183: /*
184:    Normalize a vector x with respect to a given norm as well as the
185:    sign of the first nonzero entry (assumed to be idx in process p).
186: */
187: static PetscErrorCode Normalize(Vec x,PetscReal norm,PetscInt idx,PetscInt p,PetscScalar *sign)
188: {
189:   PetscErrorCode    ierr;
190:   PetscScalar       alpha=1.0;
191:   PetscInt          first,last;
192:   PetscMPIInt       pp;
193:   PetscReal         absv;
194:   const PetscScalar *xx;

197:   VecGetOwnershipRange(x,&first,&last);
198:   if (idx>=first && idx<last) {
199:     VecGetArrayRead(x,&xx);
200:     absv = PetscAbsScalar(xx[idx-first]);
201:     if (absv>10*PETSC_MACHINE_EPSILON) alpha = xx[idx-first]/absv;
202:     VecRestoreArrayRead(x,&xx);
203:   }
204:   PetscMPIIntCast(p,&pp);
205:   MPI_Bcast(&alpha,1,MPIU_SCALAR,pp,PetscObjectComm((PetscObject)x));
206:   if (sign) *sign = alpha;
207:   alpha *= norm;
208:   VecScale(x,1.0/alpha);
209:   return(0);
210: }

212: static PetscErrorCode EPSPowerUpdateFunctionB(EPS eps,Vec x,Vec Bx)
213: {
215:   EPS_POWER      *power = (EPS_POWER*)eps->data;
216:   Mat            B;

219:   STResetMatrixState(eps->st);
220:   EPSGetOperators(eps,NULL,&B);
221:   if (B) {
222:     if (power->formFunctionB) {
223:       (*power->formFunctionB)(power->snes,x,Bx,power->formFunctionBctx);
224:     } else {
225:       MatMult(B,x,Bx);
226:     }
227:   } else {
228:     VecCopy(x,Bx);
229:   }
230:   return(0);
231: }

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

240:   STResetMatrixState(eps->st);
241:   EPSGetOperators(eps,&A,NULL);
242:   if (A) {
243:     if (power->formFunctionA) {
244:       (*power->formFunctionA)(power->snes,x,Ax,power->formFunctionActx);
245:     } else {
246:       MatMult(A,x,Ax);
247:     }
248:   } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_NULL,"Matrix A is required for an eigenvalue problem");
249:   return(0);
250: }

252: static PetscErrorCode SNESLineSearchPostheckFunction(SNESLineSearch linesearch,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
253: {
255:   SNES           snes;
256:   EPS            eps;
257:   Vec            oldx;

260:   SNESLineSearchGetSNES(linesearch,&snes);
261:   PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
262:   if (!eps) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
263:   oldx = eps->work[3];
264:   VecCopy(x,oldx);
265:   return(0);
266: }

268: static PetscErrorCode EPSPowerFormFunction_Update(SNES snes,Vec x,Vec y,void *ctx)
269: {
271:   EPS            eps;
272:   PetscReal      bx;
273:   Vec            Bx;
274:   EPS_POWER      *power;

277:   PetscObjectQuery((PetscObject)snes,"eps",(PetscObject *)&eps);
278:   if (!eps) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_NULL,"No composed EPS");
279:   power = (EPS_POWER*)eps->data;
280:   Bx = eps->work[2];
281:   if (power->formFunctionAB) {
282:     (*power->formFunctionAB)(snes,x,y,Bx,ctx);
283:   } else {
284:     EPSPowerUpdateFunctionA(eps,x,y);
285:     EPSPowerUpdateFunctionB(eps,x,Bx);
286:   }
287:   VecNorm(Bx,NORM_2,&bx);
288:   Normalize(Bx,bx,power->idx,power->p,NULL);
289:   VecAXPY(y,-1.0,Bx);
290:   return(0);
291: }

293: /*
294:    Use SNES to compute y = A^{-1}*B*x for the nonlinear problem
295: */
296: static PetscErrorCode EPSPowerApply_SNES(EPS eps,Vec x,Vec y)
297: {
299:   EPS_POWER      *power = (EPS_POWER*)eps->data;
300:   Vec            Bx;

303:   VecCopy(x,y);
304:   if (power->update) {
305:     SNESSolve(power->snes,NULL,y);
306:   } else {
307:     Bx = eps->work[2];
308:     SNESSolve(power->snes,Bx,y);
309:   }
310:   return(0);
311: }

313: /*
314:    Use nonlinear inverse power to compute an initial guess
315: */
316: static PetscErrorCode EPSPowerComputeInitialGuess_Update(EPS eps)
317: {
318:   EPS            powereps;
319:   Mat            A,B;
320:   Vec            v1,v2;
321:   SNES           snes;
322:   DM             dm,newdm;

326:   EPSCreate(PetscObjectComm((PetscObject)eps),&powereps);
327:   EPSGetOperators(eps,&A,&B);
328:   EPSSetType(powereps,EPSPOWER);
329:   EPSSetOperators(powereps,A,B);
330:   EPSSetTolerances(powereps,1e-6,4);
331:   EPSSetOptionsPrefix(powereps,((PetscObject)eps)->prefix);
332:   EPSAppendOptionsPrefix(powereps,"init_");
333:   EPSSetProblemType(powereps,EPS_GNHEP);
334:   EPSSetWhichEigenpairs(powereps,EPS_TARGET_MAGNITUDE);
335:   EPSPowerSetNonlinear(powereps,PETSC_TRUE);
336:   STSetType(powereps->st,STSINVERT);
337:   /* attach dm to initial solve */
338:   EPSPowerGetSNES(eps,&snes);
339:   SNESGetDM(snes,&dm);
340:   /* use  dmshell to temporarily store snes context */
341:   DMCreate(PetscObjectComm((PetscObject)eps),&newdm);
342:   DMSetType(newdm,DMSHELL);
343:   DMSetUp(newdm);
344:   DMCopyDMSNES(dm,newdm);
345:   EPSPowerGetSNES(powereps,&snes);
346:   SNESSetDM(snes,dm);
347:   EPSSetFromOptions(powereps);
348:   EPSSolve(powereps);
349:   BVGetColumn(eps->V,0,&v2);
350:   BVGetColumn(powereps->V,0,&v1);
351:   VecCopy(v1,v2);
352:   BVRestoreColumn(powereps->V,0,&v1);
353:   BVRestoreColumn(eps->V,0,&v2);
354:   EPSDestroy(&powereps);
355:   /* restore context back to the old nonlinear solver */
356:   DMCopyDMSNES(newdm,dm);
357:   DMDestroy(&newdm);
358:   return(0);
359: }

361: PetscErrorCode EPSSolve_Power(EPS eps)
362: {
363: #if defined(SLEPC_MISSING_LAPACK_LAEV2)
365:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LAEV2 - Lapack routine is unavailable");
366: #else
367:   PetscErrorCode      ierr;
368:   EPS_POWER           *power = (EPS_POWER*)eps->data;
369:   PetscInt            k,ld;
370:   Vec                 v,y,e,Bx;
371:   Mat                 A;
372:   KSP                 ksp;
373:   PetscReal           relerr,norm,norm1,rt1,rt2,cs1;
374:   PetscScalar         theta,rho,delta,sigma,alpha2,beta1,sn1,*T,sign;
375:   PetscBool           breakdown;
376:   KSPConvergedReason  reason;
377:   SNESConvergedReason snesreason;

380:   e = eps->work[0];
381:   y = eps->work[1];
382:   if (power->nonlinear) Bx = eps->work[2];
383:   else Bx = NULL;

385:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
386:   if (power->nonlinear) {
387:     PetscObjectCompose((PetscObject)power->snes,"eps",(PetscObject)eps);
388:     if (power->update) {
389:       EPSPowerComputeInitialGuess_Update(eps);
390:     }
391:   } else {
392:     DSGetLeadingDimension(eps->ds,&ld);
393:   }
394:   if (!power->update) {
395:     EPSGetStartVector(eps,0,NULL);
396:   }
397:   if (power->nonlinear) {
398:     BVGetColumn(eps->V,0,&v);
399:     EPSPowerUpdateFunctionB(eps,v,Bx);
400:     VecNorm(Bx,NORM_2,&norm);
401:     FirstNonzeroIdx(Bx,&power->idx,&power->p);
402:     Normalize(Bx,norm,power->idx,power->p,NULL);
403:     BVRestoreColumn(eps->V,0,&v);
404:   }

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

409:   while (eps->reason == EPS_CONVERGED_ITERATING) {
410:     eps->its++;
411:     k = eps->nconv;

413:     /* y = OP v */
414:     BVGetColumn(eps->V,k,&v);
415:     if (power->nonlinear) {
416:       VecCopy(v,eps->work[3]);
417:       EPSPowerApply_SNES(eps,v,y);
418:       VecCopy(eps->work[3],v);
419:     } else {
420:       STApply(eps->st,v,y);
421:     }
422:     BVRestoreColumn(eps->V,k,&v);

424:     /* purge previously converged eigenvectors */
425:     if (power->nonlinear) {
426:       EPSPowerUpdateFunctionB(eps,y,Bx);
427:       VecNorm(Bx,NORM_2,&norm);
428:       Normalize(Bx,norm,power->idx,power->p,&sign);
429:     } else {
430:       DSGetArray(eps->ds,DS_MAT_A,&T);
431:       BVSetActiveColumns(eps->V,0,k);
432:       BVOrthogonalizeVec(eps->V,y,T+k*ld,&norm,NULL);
433:     }

435:     /* theta = (v,y)_B */
436:     BVSetActiveColumns(eps->V,k,k+1);
437:     BVDotVec(eps->V,y,&theta);
438:     if (!power->nonlinear) {
439:       T[k+k*ld] = theta;
440:       DSRestoreArray(eps->ds,DS_MAT_A,&T);
441:     }

443:     if (power->nonlinear) theta = norm*sign;

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

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

450:       /* compute relative error as ||y-theta v||_2/|theta| */
451:       VecCopy(y,e);
452:       BVGetColumn(eps->V,k,&v);
453:       VecAXPY(e,power->nonlinear?-1.0:-theta,v);
454:       BVRestoreColumn(eps->V,k,&v);
455:       VecNorm(e,NORM_2,&relerr);
456:       relerr /= PetscAbsScalar(theta);

458:     } else {  /* RQI */

460:       /* delta = ||y||_B */
461:       delta = norm;

463:       /* compute relative error */
464:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
465:       else relerr = 1.0 / (norm*PetscAbsScalar(rho));

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

470:       /* compute new shift */
471:       if (relerr<eps->tol) {
472:         rho = sigma;  /* if converged, restore original shift */
473:         STSetShift(eps->st,rho);
474:       } else {
475:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
476:         if (power->shift_type == EPS_POWER_SHIFT_WILKINSON) {
477:           /* beta1 is the norm of the residual associated with R(v) */
478:           BVGetColumn(eps->V,k,&v);
479:           VecAXPY(v,-PetscConj(theta)/(delta*delta),y);
480:           BVRestoreColumn(eps->V,k,&v);
481:           BVScaleColumn(eps->V,k,1.0/delta);
482:           BVNormColumn(eps->V,k,NORM_2,&norm1);
483:           beta1 = norm1;

485:           /* alpha2 = (e'*A*e)/(beta1*beta1), where e is the residual */
486:           STGetMatrix(eps->st,0,&A);
487:           BVGetColumn(eps->V,k,&v);
488:           MatMult(A,v,e);
489:           VecDot(v,e,&alpha2);
490:           BVRestoreColumn(eps->V,k,&v);
491:           alpha2 = alpha2 / (beta1 * beta1);

493:           /* choose the eigenvalue of [rho beta1; beta1 alpha2] closest to rho */
494:           PetscFPTrapPush(PETSC_FP_TRAP_OFF);
495:           PetscStackCallBLAS("LAPACKlaev2",LAPACKlaev2_(&rho,&beta1,&alpha2,&rt1,&rt2,&cs1,&sn1));
496:           PetscFPTrapPop();
497:           if (PetscAbsScalar(rt1-rho) < PetscAbsScalar(rt2-rho)) rho = rt1;
498:           else rho = rt2;
499:         }
500:         /* update operator according to new shift */
501:         KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
502:         STSetShift(eps->st,rho);
503:         KSPGetConvergedReason(ksp,&reason);
504:         if (reason) {
505:           PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
506:           rho *= 1+10*PETSC_MACHINE_EPSILON;
507:           STSetShift(eps->st,rho);
508:           KSPGetConvergedReason(ksp,&reason);
509:           if (reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
510:         }
511:         KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
512:       }
513:     }
514:     eps->errest[eps->nconv] = relerr;

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

520:     if (power->update) {
521:       SNESGetConvergedReason(power->snes,&snesreason);
522:     }
523:     /* if relerr<tol, accept eigenpair */
524:     if (relerr<eps->tol || (power->update && snesreason > 0)) {
525:       eps->nconv = eps->nconv + 1;
526:       if (eps->nconv<eps->nev) {
527:         EPSGetStartVector(eps,eps->nconv,&breakdown);
528:         if (breakdown) {
529:           eps->reason = EPS_DIVERGED_BREAKDOWN;
530:           PetscInfo(eps,"Unable to generate more start vectors\n");
531:           break;
532:         }
533:       }
534:     }
535:     EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
536:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
537:   }

539:   if (power->nonlinear) {
540:     PetscObjectCompose((PetscObject)power->snes,"eps",NULL);
541:   } else {
542:     DSSetDimensions(eps->ds,eps->nconv,0,0,0);
543:     DSSetState(eps->ds,DS_STATE_RAW);
544:   }
545:   return(0);
546: #endif
547: }

549: PetscErrorCode EPSSolve_TS_Power(EPS eps)
550: {
551:   PetscErrorCode     ierr;
552:   EPS_POWER          *power = (EPS_POWER*)eps->data;
553:   PetscInt           k,ld;
554:   Vec                v,w,y,e,z;
555:   KSP                ksp;
556:   PetscReal          relerr=1.0,relerrl,delta;
557:   PetscScalar        theta,rho,alpha,sigma;
558:   PetscBool          breakdown;
559:   KSPConvergedReason reason;

562:   e = eps->work[0];
563:   v = eps->work[1];
564:   w = eps->work[2];

566:   if (power->shift_type != EPS_POWER_SHIFT_CONSTANT) { STGetKSP(eps->st,&ksp); }
567:   DSGetLeadingDimension(eps->ds,&ld);
568:   EPSGetStartVector(eps,0,NULL);
569:   BVSetRandomColumn(eps->W,0);
570:   BVBiorthonormalizeColumn(eps->V,eps->W,0,NULL);
571:   BVCopyVec(eps->V,0,v);
572:   BVCopyVec(eps->W,0,w);
573:   STGetShift(eps->st,&sigma);    /* original shift */
574:   rho = sigma;

576:   while (eps->reason == EPS_CONVERGED_ITERATING) {
577:     eps->its++;
578:     k = eps->nconv;

580:     /* y = OP v, z = OP' w */
581:     BVGetColumn(eps->V,k,&y);
582:     STApply(eps->st,v,y);
583:     BVRestoreColumn(eps->V,k,&y);
584:     BVGetColumn(eps->W,k,&z);
585:     VecConjugate(w);
586:     STApplyTranspose(eps->st,w,z);
587:     VecConjugate(w);
588:     VecConjugate(z);
589:     BVRestoreColumn(eps->W,k,&z);

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

594:     /* theta = (w,y)_B */
595:     BVSetActiveColumns(eps->V,k,k+1);
596:     BVDotVec(eps->V,w,&theta);
597:     theta = PetscConj(theta);

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

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

604:       /* compute relative errors as ||y-theta v||_2/|theta| and ||z-conj(theta) w||_2/|theta|*/
605:       BVCopyVec(eps->V,k,e);
606:       VecAXPY(e,-theta,v);
607:       VecNorm(e,NORM_2,&relerr);
608:       BVCopyVec(eps->W,k,e);
609:       VecAXPY(e,-PetscConj(theta),w);
610:       VecNorm(e,NORM_2,&relerrl);
611:       relerr = PetscMax(relerr,relerrl)/PetscAbsScalar(theta);
612:     }

614:     /* normalize */
615:     BVSetActiveColumns(eps->V,k,k+1);
616:     BVGetColumn(eps->W,k,&z);
617:     BVDotVec(eps->V,z,&alpha);
618:     BVRestoreColumn(eps->W,k,&z);
619:     delta = PetscSqrtReal(PetscAbsScalar(alpha));
620:     if (delta==0.0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Breakdown in two-sided Power/RQI");
621:     BVScaleColumn(eps->V,k,1.0/PetscConj(alpha/delta));
622:     BVScaleColumn(eps->W,k,1.0/delta);
623:     BVCopyVec(eps->V,k,v);
624:     BVCopyVec(eps->W,k,w);

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

628:       /* compute relative error */
629:       if (rho == 0.0) relerr = PETSC_MAX_REAL;
630:       else relerr = 1.0 / PetscAbsScalar(delta*rho);

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

635:       /* compute new shift */
636:       if (relerr<eps->tol) {
637:         rho = sigma;  /* if converged, restore original shift */
638:         STSetShift(eps->st,rho);
639:       } else {
640:         rho = rho + PetscConj(theta)/(delta*delta);  /* Rayleigh quotient R(v) */
641:         /* update operator according to new shift */
642:         KSPSetErrorIfNotConverged(ksp,PETSC_FALSE);
643:         STSetShift(eps->st,rho);
644:         KSPGetConvergedReason(ksp,&reason);
645:         if (reason) {
646:           PetscInfo(eps,"Factorization failed, repeat with a perturbed shift\n");
647:           rho *= 1+10*PETSC_MACHINE_EPSILON;
648:           STSetShift(eps->st,rho);
649:           KSPGetConvergedReason(ksp,&reason);
650:           if (reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Second factorization failed");
651:         }
652:         KSPSetErrorIfNotConverged(ksp,PETSC_TRUE);
653:       }
654:     }
655:     eps->errest[eps->nconv] = relerr;

657:     /* if relerr<tol, accept eigenpair */
658:     if (relerr<eps->tol) {
659:       eps->nconv = eps->nconv + 1;
660:       if (eps->nconv<eps->nev) {
661:         EPSGetStartVector(eps,eps->nconv,&breakdown);
662:         if (breakdown) {
663:           eps->reason = EPS_DIVERGED_BREAKDOWN;
664:           PetscInfo(eps,"Unable to generate more start vectors\n");
665:           break;
666:         }
667:         BVSetRandomColumn(eps->W,eps->nconv);
668:         BVBiorthonormalizeColumn(eps->V,eps->W,eps->nconv,NULL);
669:       }
670:     }
671:     EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,PetscMin(eps->nconv+1,eps->nev));
672:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
673:   }

675:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
676:   DSSetState(eps->ds,DS_STATE_RAW);
677:   return(0);
678: }

680: PetscErrorCode EPSStopping_Power(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
681: {
683:   EPS_POWER      *power = (EPS_POWER*)eps->data;
684:   SNESConvergedReason snesreason;

687:   if (power->update) {
688:     SNESGetConvergedReason(power->snes,&snesreason);
689:     if (snesreason < 0) {
690:       *reason = EPS_DIVERGED_BREAKDOWN;
691:       return(0);
692:     }
693:   }
694:   EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,ctx);
695:   return(0);
696: }

698: PetscErrorCode EPSBackTransform_Power(EPS eps)
699: {
701:   EPS_POWER      *power = (EPS_POWER*)eps->data;

704:   if (power->nonlinear) eps->eigr[0] = 1.0/eps->eigr[0];
705:   else if (power->shift_type == EPS_POWER_SHIFT_CONSTANT) {
706:     EPSBackTransform_Default(eps);
707:   }
708:   return(0);
709: }

711: PetscErrorCode EPSSetFromOptions_Power(PetscOptionItems *PetscOptionsObject,EPS eps)
712: {
713:   PetscErrorCode    ierr;
714:   EPS_POWER         *power = (EPS_POWER*)eps->data;
715:   PetscBool         flg,val;
716:   EPSPowerShiftType shift;

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

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

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

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

730:   PetscOptionsTail();
731:   return(0);
732: }

734: static PetscErrorCode EPSPowerSetShiftType_Power(EPS eps,EPSPowerShiftType shift)
735: {
736:   EPS_POWER *power = (EPS_POWER*)eps->data;

739:   switch (shift) {
740:     case EPS_POWER_SHIFT_CONSTANT:
741:     case EPS_POWER_SHIFT_RAYLEIGH:
742:     case EPS_POWER_SHIFT_WILKINSON:
743:       if (power->shift_type != shift) {
744:         power->shift_type = shift;
745:         eps->state = EPS_STATE_INITIAL;
746:       }
747:       break;
748:     default:
749:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid shift type");
750:   }
751:   return(0);
752: }

754: /*@
755:    EPSPowerSetShiftType - Sets the type of shifts used during the power
756:    iteration. This can be used to emulate the Rayleigh Quotient Iteration
757:    (RQI) method.

759:    Logically Collective on eps

761:    Input Parameters:
762: +  eps - the eigenproblem solver context
763: -  shift - the type of shift

765:    Options Database Key:
766: .  -eps_power_shift_type - Sets the shift type (either 'constant' or
767:                            'rayleigh' or 'wilkinson')

769:    Notes:
770:    By default, shifts are constant (EPS_POWER_SHIFT_CONSTANT) and the iteration
771:    is the simple power method (or inverse iteration if a shift-and-invert
772:    transformation is being used).

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

778:    Level: advanced

780: .seealso: EPSPowerGetShiftType(), STSetShift(), EPSPowerShiftType
781: @*/
782: PetscErrorCode EPSPowerSetShiftType(EPS eps,EPSPowerShiftType shift)
783: {

789:   PetscTryMethod(eps,"EPSPowerSetShiftType_C",(EPS,EPSPowerShiftType),(eps,shift));
790:   return(0);
791: }

793: static PetscErrorCode EPSPowerGetShiftType_Power(EPS eps,EPSPowerShiftType *shift)
794: {
795:   EPS_POWER *power = (EPS_POWER*)eps->data;

798:   *shift = power->shift_type;
799:   return(0);
800: }

802: /*@
803:    EPSPowerGetShiftType - Gets the type of shifts used during the power
804:    iteration.

806:    Not Collective

808:    Input Parameter:
809: .  eps - the eigenproblem solver context

811:    Input Parameter:
812: .  shift - the type of shift

814:    Level: advanced

816: .seealso: EPSPowerSetShiftType(), EPSPowerShiftType
817: @*/
818: PetscErrorCode EPSPowerGetShiftType(EPS eps,EPSPowerShiftType *shift)
819: {

825:   PetscUseMethod(eps,"EPSPowerGetShiftType_C",(EPS,EPSPowerShiftType*),(eps,shift));
826:   return(0);
827: }

829: static PetscErrorCode EPSPowerSetNonlinear_Power(EPS eps,PetscBool nonlinear)
830: {
831:   EPS_POWER *power = (EPS_POWER*)eps->data;

834:   if (power->nonlinear != nonlinear) {
835:     power->nonlinear = nonlinear;
836:     eps->useds = PetscNot(nonlinear);
837:     eps->state = EPS_STATE_INITIAL;
838:   }
839:   return(0);
840: }

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

845:    Logically Collective on eps

847:    Input Parameters:
848: +  eps - the eigenproblem solver context
849: -  nonlinear - whether the problem is nonlinear or not

851:    Options Database Key:
852: .  -eps_power_nonlinear - Sets the nonlinear flag

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

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

868:    Level: advanced

870: .seealso: EPSPowerGetNonlinear(), EPSSetOperators()
871: @*/
872: PetscErrorCode EPSPowerSetNonlinear(EPS eps,PetscBool nonlinear)
873: {

879:   PetscTryMethod(eps,"EPSPowerSetNonlinear_C",(EPS,PetscBool),(eps,nonlinear));
880:   return(0);
881: }

883: static PetscErrorCode EPSPowerGetNonlinear_Power(EPS eps,PetscBool *nonlinear)
884: {
885:   EPS_POWER *power = (EPS_POWER*)eps->data;

888:   *nonlinear = power->nonlinear;
889:   return(0);
890: }

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

895:    Not Collective

897:    Input Parameter:
898: .  eps - the eigenproblem solver context

900:    Input Parameter:
901: .  nonlinear - the nonlinear flag

903:    Level: advanced

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

914:   PetscUseMethod(eps,"EPSPowerGetNonlinear_C",(EPS,PetscBool*),(eps,nonlinear));
915:   return(0);
916: }

918: static PetscErrorCode EPSPowerSetUpdate_Power(EPS eps,PetscBool update)
919: {
920:   EPS_POWER *power = (EPS_POWER*)eps->data;

923:   if (!power->nonlinear) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"This option does not make sense for linear problems");
924:   power->update = update;
925:   eps->state = EPS_STATE_INITIAL;
926:   return(0);
927: }

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

933:    Logically Collective on eps

935:    Input Parameters:
936: +  eps - the eigenproblem solver context
937: -  update - whether the residual is updated monolithically or not

939:    Options Database Key:
940: .  -eps_power_update - Sets the update flag

942:    Level: advanced

944: .seealso: EPSPowerGetUpdate(), EPSPowerGetNonlinear(), EPSSetOperators()
945: @*/
946: PetscErrorCode EPSPowerSetUpdate(EPS eps,PetscBool update)
947: {

953:   PetscTryMethod(eps,"EPSPowerSetUpdate_C",(EPS,PetscBool),(eps,update));
954:   return(0);
955: }

957: static PetscErrorCode EPSPowerGetUpdate_Power(EPS eps,PetscBool *update)
958: {
959:   EPS_POWER *power = (EPS_POWER*)eps->data;

962:   *update = power->update;
963:   return(0);
964: }

966: /*@
967:    EPSPowerGetUpdate - Returns a flag indicating if the residual is updated monolithically
968:    for nonlinear problems.

970:    Not Collective

972:    Input Parameter:
973: .  eps - the eigenproblem solver context

975:    Input Parameter:
976: .  update - the update flag

978:    Level: advanced

980: .seealso: EPSPowerSetUpdate(), EPSPowerSetNonlinear()
981: @*/
982: PetscErrorCode EPSPowerGetUpdate(EPS eps,PetscBool *update)
983: {

989:   PetscUseMethod(eps,"EPSPowerGetUpdate_C",(EPS,PetscBool*),(eps,update));
990:   return(0);
991: }

993: static PetscErrorCode EPSPowerSetSNES_Power(EPS eps,SNES snes)
994: {
996:   EPS_POWER      *power = (EPS_POWER*)eps->data;

999:   PetscObjectReference((PetscObject)snes);
1000:   SNESDestroy(&power->snes);
1001:   power->snes = snes;
1002:   PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
1003:   eps->state = EPS_STATE_INITIAL;
1004:   return(0);
1005: }

1007: /*@
1008:    EPSPowerSetSNES - Associate a nonlinear solver object (SNES) to the
1009:    eigenvalue solver (to be used in nonlinear inverse iteration).

1011:    Collective on eps

1013:    Input Parameters:
1014: +  eps  - the eigenvalue solver
1015: -  snes - the nonlinear solver object

1017:    Level: advanced

1019: .seealso: EPSPowerGetSNES()
1020: @*/
1021: PetscErrorCode EPSPowerSetSNES(EPS eps,SNES snes)
1022: {

1029:   PetscTryMethod(eps,"EPSPowerSetSNES_C",(EPS,SNES),(eps,snes));
1030:   return(0);
1031: }

1033: static PetscErrorCode EPSPowerGetSNES_Power(EPS eps,SNES *snes)
1034: {
1036:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1039:   if (!power->snes) {
1040:     SNESCreate(PetscObjectComm((PetscObject)eps),&power->snes);
1041:     PetscObjectIncrementTabLevel((PetscObject)power->snes,(PetscObject)eps,1);
1042:     SNESSetOptionsPrefix(power->snes,((PetscObject)eps)->prefix);
1043:     SNESAppendOptionsPrefix(power->snes,"eps_power_");
1044:     PetscLogObjectParent((PetscObject)eps,(PetscObject)power->snes);
1045:     PetscObjectSetOptions((PetscObject)power->snes,((PetscObject)eps)->options);
1046:   }
1047:   *snes = power->snes;
1048:   return(0);
1049: }

1051: /*@
1052:    EPSPowerGetSNES - Retrieve the nonlinear solver object (SNES) associated
1053:    with the eigenvalue solver.

1055:    Not Collective

1057:    Input Parameter:
1058: .  eps - the eigenvalue solver

1060:    Output Parameter:
1061: .  snes - the nonlinear solver object

1063:    Level: advanced

1065: .seealso: EPSPowerSetSNES()
1066: @*/
1067: PetscErrorCode EPSPowerGetSNES(EPS eps,SNES *snes)
1068: {

1074:   PetscUseMethod(eps,"EPSPowerGetSNES_C",(EPS,SNES*),(eps,snes));
1075:   return(0);
1076: }

1078: PetscErrorCode EPSReset_Power(EPS eps)
1079: {
1081:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1084:   if (power->snes) { SNESReset(power->snes); }
1085:   return(0);
1086: }

1088: PetscErrorCode EPSDestroy_Power(EPS eps)
1089: {
1091:   EPS_POWER      *power = (EPS_POWER*)eps->data;

1094:   if (power->nonlinear) {
1095:     SNESDestroy(&power->snes);
1096:   }
1097:   PetscFree(eps->data);
1098:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",NULL);
1099:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",NULL);
1100:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",NULL);
1101:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",NULL);
1102:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",NULL);
1103:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",NULL);
1104:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",NULL);
1105:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",NULL);
1106:   return(0);
1107: }

1109: PetscErrorCode EPSView_Power(EPS eps,PetscViewer viewer)
1110: {
1112:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1113:   PetscBool      isascii;

1116:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1117:   if (isascii) {
1118:     if (power->nonlinear) {
1119:       PetscViewerASCIIPrintf(viewer,"  using nonlinear inverse iteration\n");
1120:       if (power->update) {
1121:         PetscViewerASCIIPrintf(viewer,"  updating the residual monolithically\n");
1122:       }
1123:       if (!power->snes) { EPSPowerGetSNES(eps,&power->snes); }
1124:       PetscViewerASCIIPushTab(viewer);
1125:       SNESView(power->snes,viewer);
1126:       PetscViewerASCIIPopTab(viewer);
1127:     } else {
1128:       PetscViewerASCIIPrintf(viewer,"  %s shifts\n",EPSPowerShiftTypes[power->shift_type]);
1129:     }
1130:   }
1131:   return(0);
1132: }

1134: PetscErrorCode EPSComputeVectors_Power(EPS eps)
1135: {
1137:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1138:   PetscReal      norm;
1139:   PetscInt       i;

1142:   if (eps->twosided) {
1143:     EPSComputeVectors_Twosided(eps);
1144:     /* normalize (no need to take care of 2x2 blocks */
1145:     for (i=0;i<eps->nconv;i++) {
1146:       BVNormColumn(eps->V,i,NORM_2,&norm);
1147:       BVScaleColumn(eps->V,i,1.0/norm);
1148:       BVNormColumn(eps->W,i,NORM_2,&norm);
1149:       BVScaleColumn(eps->W,i,1.0/norm);
1150:     }
1151:   } else if (!power->nonlinear) {
1152:     EPSComputeVectors_Schur(eps);
1153:   }
1154:   return(0);
1155: }

1157: PetscErrorCode EPSSetDefaultST_Power(EPS eps)
1158: {
1160:   EPS_POWER      *power = (EPS_POWER*)eps->data;
1161:   KSP            ksp;
1162:   PC             pc;

1165:   if (power->nonlinear) {
1166:     STGetKSP(eps->st,&ksp);
1167:     KSPGetPC(ksp,&pc);
1168:     PCSetType(pc,PCNONE);
1169:   }
1170:   return(0);
1171: }

1173: SLEPC_EXTERN PetscErrorCode EPSCreate_Power(EPS eps)
1174: {
1175:   EPS_POWER      *ctx;

1179:   PetscNewLog(eps,&ctx);
1180:   eps->data = (void*)ctx;

1182:   eps->useds = PETSC_TRUE;
1183:   eps->hasts = PETSC_TRUE;
1184:   eps->categ = EPS_CATEGORY_OTHER;

1186:   eps->ops->setup          = EPSSetUp_Power;
1187:   eps->ops->setfromoptions = EPSSetFromOptions_Power;
1188:   eps->ops->reset          = EPSReset_Power;
1189:   eps->ops->destroy        = EPSDestroy_Power;
1190:   eps->ops->view           = EPSView_Power;
1191:   eps->ops->backtransform  = EPSBackTransform_Power;
1192:   eps->ops->computevectors = EPSComputeVectors_Power;
1193:   eps->ops->setdefaultst   = EPSSetDefaultST_Power;
1194:   eps->stopping            = EPSStopping_Power;

1196:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetShiftType_C",EPSPowerSetShiftType_Power);
1197:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetShiftType_C",EPSPowerGetShiftType_Power);
1198:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetNonlinear_C",EPSPowerSetNonlinear_Power);
1199:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetNonlinear_C",EPSPowerGetNonlinear_Power);
1200:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetUpdate_C",EPSPowerSetUpdate_Power);
1201:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetUpdate_C",EPSPowerGetUpdate_Power);
1202:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerSetSNES_C",EPSPowerSetSNES_Power);
1203:   PetscObjectComposeFunction((PetscObject)eps,"EPSPowerGetSNES_C",EPSPowerGetSNES_Power);
1204:   return(0);
1205: }