LCOV - code coverage report
Current view: top level - eps/impls/power - power.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 553 599 92.3 %
Date: 2020-07-05 02:54:12 Functions: 33 36 91.7 %

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

Generated by: LCOV version 1.13