Actual source code: linear.c

slepc-3.21.0 2024-03-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:    Explicit linearization for polynomial eigenproblems
 12: */

 14: #include <slepc/private/pepimpl.h>
 15: #include "linear.h"

 17: static PetscErrorCode MatMult_Linear_Shift(Mat M,Vec x,Vec y)
 18: {
 19:   PEP_LINEAR        *ctx;
 20:   PEP               pep;
 21:   const PetscScalar *px;
 22:   PetscScalar       *py,a,sigma=0.0;
 23:   PetscInt          nmat,deg,i,m;
 24:   Vec               x1,x2,x3,y1,aux;
 25:   PetscReal         *ca,*cb,*cg;
 26:   PetscBool         flg;

 28:   PetscFunctionBegin;
 29:   PetscCall(MatShellGetContext(M,&ctx));
 30:   pep = ctx->pep;
 31:   PetscCall(STGetTransform(pep->st,&flg));
 32:   if (!flg) PetscCall(STGetShift(pep->st,&sigma));
 33:   nmat = pep->nmat;
 34:   deg = nmat-1;
 35:   m = pep->nloc;
 36:   ca = pep->pbc;
 37:   cb = pep->pbc+nmat;
 38:   cg = pep->pbc+2*nmat;
 39:   x1=ctx->w[0];x2=ctx->w[1];x3=ctx->w[2];y1=ctx->w[3];aux=ctx->w[4];

 41:   PetscCall(VecSet(y,0.0));
 42:   PetscCall(VecGetArrayRead(x,&px));
 43:   PetscCall(VecGetArray(y,&py));
 44:   a = 1.0;

 46:   /* first block */
 47:   PetscCall(VecPlaceArray(x2,px));
 48:   PetscCall(VecPlaceArray(x3,px+m));
 49:   PetscCall(VecPlaceArray(y1,py));
 50:   PetscCall(VecAXPY(y1,cb[0]-sigma,x2));
 51:   PetscCall(VecAXPY(y1,ca[0],x3));
 52:   PetscCall(VecResetArray(x2));
 53:   PetscCall(VecResetArray(x3));
 54:   PetscCall(VecResetArray(y1));

 56:   /* inner blocks */
 57:   for (i=1;i<deg-1;i++) {
 58:     PetscCall(VecPlaceArray(x1,px+(i-1)*m));
 59:     PetscCall(VecPlaceArray(x2,px+i*m));
 60:     PetscCall(VecPlaceArray(x3,px+(i+1)*m));
 61:     PetscCall(VecPlaceArray(y1,py+i*m));
 62:     PetscCall(VecAXPY(y1,cg[i],x1));
 63:     PetscCall(VecAXPY(y1,cb[i]-sigma,x2));
 64:     PetscCall(VecAXPY(y1,ca[i],x3));
 65:     PetscCall(VecResetArray(x1));
 66:     PetscCall(VecResetArray(x2));
 67:     PetscCall(VecResetArray(x3));
 68:     PetscCall(VecResetArray(y1));
 69:   }

 71:   /* last block */
 72:   PetscCall(VecPlaceArray(y1,py+(deg-1)*m));
 73:   for (i=0;i<deg;i++) {
 74:     PetscCall(VecPlaceArray(x1,px+i*m));
 75:     PetscCall(STMatMult(pep->st,i,x1,aux));
 76:     PetscCall(VecAXPY(y1,a,aux));
 77:     PetscCall(VecResetArray(x1));
 78:     a *= pep->sfactor;
 79:   }
 80:   PetscCall(VecCopy(y1,aux));
 81:   PetscCall(STMatSolve(pep->st,aux,y1));
 82:   PetscCall(VecScale(y1,-ca[deg-1]/a));
 83:   PetscCall(VecPlaceArray(x1,px+(deg-2)*m));
 84:   PetscCall(VecPlaceArray(x2,px+(deg-1)*m));
 85:   PetscCall(VecAXPY(y1,cg[deg-1],x1));
 86:   PetscCall(VecAXPY(y1,cb[deg-1]-sigma,x2));
 87:   PetscCall(VecResetArray(x1));
 88:   PetscCall(VecResetArray(x2));
 89:   PetscCall(VecResetArray(y1));

 91:   PetscCall(VecRestoreArrayRead(x,&px));
 92:   PetscCall(VecRestoreArray(y,&py));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

 96: static PetscErrorCode MatMult_Linear_Sinvert(Mat M,Vec x,Vec y)
 97: {
 98:   PEP_LINEAR        *ctx;
 99:   PEP               pep;
100:   const PetscScalar *px;
101:   PetscScalar       *py,a,sigma,t=1.0,tp=0.0,tt;
102:   PetscInt          nmat,deg,i,m;
103:   Vec               x1,y1,y2,y3,aux,aux2;
104:   PetscReal         *ca,*cb,*cg;

106:   PetscFunctionBegin;
107:   PetscCall(MatShellGetContext(M,&ctx));
108:   pep = ctx->pep;
109:   nmat = pep->nmat;
110:   deg = nmat-1;
111:   m = pep->nloc;
112:   ca = pep->pbc;
113:   cb = pep->pbc+nmat;
114:   cg = pep->pbc+2*nmat;
115:   x1=ctx->w[0];y1=ctx->w[1];y2=ctx->w[2];y3=ctx->w[3];aux=ctx->w[4];aux2=ctx->w[5];
116:   PetscCall(EPSGetTarget(ctx->eps,&sigma));
117:   PetscCall(VecSet(y,0.0));
118:   PetscCall(VecGetArrayRead(x,&px));
119:   PetscCall(VecGetArray(y,&py));
120:   a = pep->sfactor;

122:   /* first block */
123:   PetscCall(VecPlaceArray(x1,px));
124:   PetscCall(VecPlaceArray(y1,py+m));
125:   PetscCall(VecCopy(x1,y1));
126:   PetscCall(VecScale(y1,1.0/ca[0]));
127:   PetscCall(VecResetArray(x1));
128:   PetscCall(VecResetArray(y1));

130:   /* second block */
131:   if (deg>2) {
132:     PetscCall(VecPlaceArray(x1,px+m));
133:     PetscCall(VecPlaceArray(y1,py+m));
134:     PetscCall(VecPlaceArray(y2,py+2*m));
135:     PetscCall(VecCopy(x1,y2));
136:     PetscCall(VecAXPY(y2,sigma-cb[1],y1));
137:     PetscCall(VecScale(y2,1.0/ca[1]));
138:     PetscCall(VecResetArray(x1));
139:     PetscCall(VecResetArray(y1));
140:     PetscCall(VecResetArray(y2));
141:   }

143:   /* inner blocks */
144:   for (i=2;i<deg-1;i++) {
145:     PetscCall(VecPlaceArray(x1,px+i*m));
146:     PetscCall(VecPlaceArray(y1,py+(i-1)*m));
147:     PetscCall(VecPlaceArray(y2,py+i*m));
148:     PetscCall(VecPlaceArray(y3,py+(i+1)*m));
149:     PetscCall(VecCopy(x1,y3));
150:     PetscCall(VecAXPY(y3,sigma-cb[i],y2));
151:     PetscCall(VecAXPY(y3,-cg[i],y1));
152:     PetscCall(VecScale(y3,1.0/ca[i]));
153:     PetscCall(VecResetArray(x1));
154:     PetscCall(VecResetArray(y1));
155:     PetscCall(VecResetArray(y2));
156:     PetscCall(VecResetArray(y3));
157:   }

159:   /* last block */
160:   PetscCall(VecPlaceArray(y1,py));
161:   for (i=0;i<deg-2;i++) {
162:     PetscCall(VecPlaceArray(y2,py+(i+1)*m));
163:     PetscCall(STMatMult(pep->st,i+1,y2,aux));
164:     PetscCall(VecAXPY(y1,a,aux));
165:     PetscCall(VecResetArray(y2));
166:     a *= pep->sfactor;
167:   }
168:   i = deg-2;
169:   PetscCall(VecPlaceArray(y2,py+(i+1)*m));
170:   PetscCall(VecPlaceArray(y3,py+i*m));
171:   PetscCall(VecCopy(y2,aux2));
172:   PetscCall(VecAXPY(aux2,cg[i+1]/ca[i+1],y3));
173:   PetscCall(STMatMult(pep->st,i+1,aux2,aux));
174:   PetscCall(VecAXPY(y1,a,aux));
175:   PetscCall(VecResetArray(y2));
176:   PetscCall(VecResetArray(y3));
177:   a *= pep->sfactor;
178:   i = deg-1;
179:   PetscCall(VecPlaceArray(x1,px+i*m));
180:   PetscCall(VecPlaceArray(y3,py+i*m));
181:   PetscCall(VecCopy(x1,aux2));
182:   PetscCall(VecAXPY(aux2,sigma-cb[i],y3));
183:   PetscCall(VecScale(aux2,1.0/ca[i]));
184:   PetscCall(STMatMult(pep->st,i+1,aux2,aux));
185:   PetscCall(VecAXPY(y1,a,aux));
186:   PetscCall(VecResetArray(x1));
187:   PetscCall(VecResetArray(y3));

189:   PetscCall(VecCopy(y1,aux));
190:   PetscCall(STMatSolve(pep->st,aux,y1));
191:   PetscCall(VecScale(y1,-1.0));

193:   /* final update */
194:   for (i=1;i<deg;i++) {
195:     PetscCall(VecPlaceArray(y2,py+i*m));
196:     tt = t;
197:     t = ((sigma-cb[i-1])*t-cg[i-1]*tp)/ca[i-1]; /* i-th basis polynomial */
198:     tp = tt;
199:     PetscCall(VecAXPY(y2,t,y1));
200:     PetscCall(VecResetArray(y2));
201:   }
202:   PetscCall(VecResetArray(y1));

204:   PetscCall(VecRestoreArrayRead(x,&px));
205:   PetscCall(VecRestoreArray(y,&py));
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: static PetscErrorCode BackTransform_Linear(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
210: {
211:   PEP_LINEAR     *ctx;
212:   ST             stctx;

214:   PetscFunctionBegin;
215:   PetscCall(STShellGetContext(st,&ctx));
216:   PetscCall(PEPGetST(ctx->pep,&stctx));
217:   PetscCall(STBackTransform(stctx,n,eigr,eigi));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /*
222:    Dummy backtransform operation
223:  */
224: static PetscErrorCode BackTransform_Skip(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
225: {
226:   PetscFunctionBegin;
227:   PetscFunctionReturn(PETSC_SUCCESS);
228: }

230: static PetscErrorCode Apply_Linear(ST st,Vec x,Vec y)
231: {
232:   PEP_LINEAR     *ctx;

234:   PetscFunctionBegin;
235:   PetscCall(STShellGetContext(st,&ctx));
236:   PetscCall(MatMult(ctx->A,x,y));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: static PetscErrorCode PEPSetUp_Linear(PEP pep)
241: {
242:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
243:   ST             st;
244:   PetscInt       i=0,deg=pep->nmat-1;
245:   EPSWhich       which = EPS_LARGEST_MAGNITUDE;
246:   EPSProblemType ptype;
247:   PetscBool      trackall,istrivial,transf,sinv,ks;
248:   PetscScalar    sigma,*epsarray,*peparray;
249:   Vec            veps,w=NULL;
250:   /* function tables */
251:   PetscErrorCode (*fcreate[][2])(MPI_Comm,PEP_LINEAR*,Mat*) = {
252:     { MatCreateExplicit_Linear_NA, MatCreateExplicit_Linear_NB },
253:     { MatCreateExplicit_Linear_SA, MatCreateExplicit_Linear_SB },
254:     { MatCreateExplicit_Linear_HA, MatCreateExplicit_Linear_HB },
255:   };

257:   PetscFunctionBegin;
258:   PEPCheckShiftSinvert(pep);
259:   PEPCheckUnsupported(pep,PEP_FEATURE_STOPPING);
260:   PEPCheckIgnored(pep,PEP_FEATURE_CONVERGENCE);
261:   PetscCall(STGetTransform(pep->st,&transf));
262:   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
263:   if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
264:   PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
265:   PetscCall(STSetUp(pep->st));
266:   if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
267:   PetscCall(EPSGetST(ctx->eps,&st));
268:   if (!transf && !ctx->usereps) PetscCall(EPSSetTarget(ctx->eps,pep->target));
269:   if (sinv && !transf && !ctx->usereps) PetscCall(STSetDefaultShift(st,pep->target));
270:   /* compute scale factor if not set by user */
271:   PetscCall(PEPComputeScaleFactor(pep));

273:   if (ctx->explicitmatrix) {
274:     PEPCheckQuadraticCondition(pep,PETSC_TRUE," (with explicit matrix)");
275:     PEPCheckUnsupportedCondition(pep,PEP_FEATURE_NONMONOMIAL,PETSC_TRUE," (with explicit matrix)");
276:     PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
277:     PetscCheck(pep->scale!=PEP_SCALE_DIAGONAL && pep->scale!=PEP_SCALE_BOTH,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Diagonal scaling not allowed in PEPLINEAR with explicit matrices");
278:     if (sinv && !transf) PetscCall(STSetType(st,STSINVERT));
279:     PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
280:     PetscCall(STGetMatrixTransformed(pep->st,0,&ctx->K));
281:     PetscCall(STGetMatrixTransformed(pep->st,1,&ctx->C));
282:     PetscCall(STGetMatrixTransformed(pep->st,2,&ctx->M));
283:     ctx->sfactor = pep->sfactor;
284:     ctx->dsfactor = pep->dsfactor;

286:     PetscCall(MatDestroy(&ctx->A));
287:     PetscCall(MatDestroy(&ctx->B));
288:     PetscCall(VecDestroy(&ctx->w[0]));
289:     PetscCall(VecDestroy(&ctx->w[1]));
290:     PetscCall(VecDestroy(&ctx->w[2]));
291:     PetscCall(VecDestroy(&ctx->w[3]));

293:     switch (pep->problem_type) {
294:       case PEP_GENERAL:    i = 0; break;
295:       case PEP_HERMITIAN:
296:       case PEP_HYPERBOLIC: i = 1; break;
297:       case PEP_GYROSCOPIC: i = 2; break;
298:     }

300:     PetscCall((*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A));
301:     PetscCall((*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B));

303:   } else {   /* implicit matrix */
304:     PetscCheck(pep->problem_type==PEP_GENERAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Must use the explicit matrix option if problem type is not general");
305:     if (!((PetscObject)ctx->eps)->type_name) PetscCall(EPSSetType(ctx->eps,EPSKRYLOVSCHUR));
306:     else {
307:       PetscCall(PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks));
308:       PetscCheck(ks,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
309:     }
310:     PetscCheck(ctx->alpha==1.0 && ctx->beta==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option does not support setting alpha,beta parameters of the linearization");
311:     PetscCall(STSetType(st,STSHELL));
312:     PetscCall(STShellSetContext(st,ctx));
313:     if (!transf) PetscCall(STShellSetBackTransform(st,BackTransform_Linear));
314:     else PetscCall(STShellSetBackTransform(st,BackTransform_Skip));
315:     PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]));
316:     PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]));
317:     PetscCall(MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]));
318:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A));
319:     if (sinv && !transf) PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert));
320:     else PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift));
321:     PetscCall(STShellSetApply(st,Apply_Linear));
322:     ctx->pep = pep;

324:     PetscCall(PEPBasisCoefficients(pep,pep->pbc));
325:     if (!transf) {
326:       PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
327:       if (sinv) PetscCall(PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL));
328:       else {
329:         for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
330:         pep->solvematcoeffs[deg] = 1.0;
331:       }
332:       PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
333:       PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
334:     }
335:     if (pep->sfactor!=1.0) {
336:       for (i=0;i<pep->nmat;i++) {
337:         pep->pbc[pep->nmat+i] /= pep->sfactor;
338:         pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
339:       }
340:     }
341:   }

343:   PetscCall(EPSSetOperators(ctx->eps,ctx->A,ctx->B));
344:   PetscCall(EPSGetProblemType(ctx->eps,&ptype));
345:   if (!ptype) {
346:     if (ctx->explicitmatrix) PetscCall(EPSSetProblemType(ctx->eps,EPS_GNHEP));
347:     else PetscCall(EPSSetProblemType(ctx->eps,EPS_NHEP));
348:   }
349:   if (!ctx->usereps) {
350:     if (transf) which = EPS_LARGEST_MAGNITUDE;
351:     else {
352:       switch (pep->which) {
353:         case PEP_LARGEST_MAGNITUDE:  which = EPS_LARGEST_MAGNITUDE; break;
354:         case PEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
355:         case PEP_LARGEST_REAL:       which = EPS_LARGEST_REAL; break;
356:         case PEP_SMALLEST_REAL:      which = EPS_SMALLEST_REAL; break;
357:         case PEP_LARGEST_IMAGINARY:  which = EPS_LARGEST_IMAGINARY; break;
358:         case PEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
359:         case PEP_TARGET_MAGNITUDE:   which = EPS_TARGET_MAGNITUDE; break;
360:         case PEP_TARGET_REAL:        which = EPS_TARGET_REAL; break;
361:         case PEP_TARGET_IMAGINARY:   which = EPS_TARGET_IMAGINARY; break;
362:         case PEP_ALL:                which = EPS_ALL; break;
363:         case PEP_WHICH_USER:         which = EPS_WHICH_USER;
364:           PetscCall(EPSSetEigenvalueComparison(ctx->eps,pep->sc->comparison,pep->sc->comparisonctx));
365:           break;
366:       }
367:     }
368:     PetscCall(EPSSetWhichEigenpairs(ctx->eps,which));

370:     PetscCall(EPSSetDimensions(ctx->eps,pep->nev,pep->ncv,pep->mpd));
371:     PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(pep->tol),pep->max_it));
372:   }
373:   PetscCall(RGIsTrivial(pep->rg,&istrivial));
374:   if (!istrivial) {
375:     PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
376:     PetscCall(EPSSetRG(ctx->eps,pep->rg));
377:   }
378:   /* Transfer the trackall option from pep to eps */
379:   PetscCall(PEPGetTrackAll(pep,&trackall));
380:   PetscCall(EPSSetTrackAll(ctx->eps,trackall));

382:   /* temporary change of target */
383:   if (pep->sfactor!=1.0) {
384:     PetscCall(EPSGetTarget(ctx->eps,&sigma));
385:     PetscCall(EPSSetTarget(ctx->eps,sigma/pep->sfactor));
386:   }

388:   /* process initial vector */
389:   if (pep->nini<0) {
390:     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps));
391:     PetscCall(VecGetArray(veps,&epsarray));
392:     for (i=0;i<deg;i++) {
393:       if (i<-pep->nini) {
394:         PetscCall(VecGetArray(pep->IS[i],&peparray));
395:         PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
396:         PetscCall(VecRestoreArray(pep->IS[i],&peparray));
397:       } else {
398:         if (!w) PetscCall(VecDuplicate(pep->IS[0],&w));
399:         PetscCall(VecSetRandom(w,NULL));
400:         PetscCall(VecGetArray(w,&peparray));
401:         PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
402:         PetscCall(VecRestoreArray(w,&peparray));
403:       }
404:     }
405:     PetscCall(VecRestoreArray(veps,&epsarray));
406:     PetscCall(EPSSetInitialSpace(ctx->eps,1,&veps));
407:     PetscCall(VecDestroy(&veps));
408:     PetscCall(VecDestroy(&w));
409:     PetscCall(SlepcBasisDestroy_Private(&pep->nini,&pep->IS));
410:   }

412:   PetscCall(EPSSetUp(ctx->eps));
413:   PetscCall(EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd));
414:   PetscCall(EPSGetTolerances(ctx->eps,NULL,&pep->max_it));
415:   PetscCall(PEPAllocateSolution(pep,0));
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: /*
420:    PEPLinearExtract_Residual - Auxiliary routine that copies the solution of the
421:    linear eigenproblem to the PEP object. The eigenvector of the generalized
422:    problem is supposed to be
423:                                z = [  x  ]
424:                                    [ l*x ]
425:    The eigenvector is taken from z(1:n) or z(n+1:2*n) depending on the explicitly
426:    computed residual norm.
427:    Finally, x is normalized so that ||x||_2 = 1.
428: */
429: static PetscErrorCode PEPLinearExtract_Residual(PEP pep,EPS eps)
430: {
431:   PetscInt          i,k;
432:   const PetscScalar *px;
433:   PetscScalar       *er=pep->eigr,*ei=pep->eigi;
434:   PetscReal         rn1,rn2;
435:   Vec               xr,xi=NULL,wr;
436:   Mat               A;
437: #if !defined(PETSC_USE_COMPLEX)
438:   Vec               wi;
439:   const PetscScalar *py;
440: #endif

442:   PetscFunctionBegin;
443: #if defined(PETSC_USE_COMPLEX)
444:   PetscCall(PEPSetWorkVecs(pep,2));
445: #else
446:   PetscCall(PEPSetWorkVecs(pep,4));
447: #endif
448:   PetscCall(EPSGetOperators(eps,&A,NULL));
449:   PetscCall(MatCreateVecs(A,&xr,NULL));
450:   PetscCall(MatCreateVecsEmpty(pep->A[0],&wr,NULL));
451: #if !defined(PETSC_USE_COMPLEX)
452:   PetscCall(VecDuplicate(xr,&xi));
453:   PetscCall(VecDuplicateEmpty(wr,&wi));
454: #endif
455:   for (i=0;i<pep->nconv;i++) {
456:     PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,xr,xi));
457: #if !defined(PETSC_USE_COMPLEX)
458:     if (ei[i]!=0.0) {   /* complex conjugate pair */
459:       PetscCall(VecGetArrayRead(xr,&px));
460:       PetscCall(VecGetArrayRead(xi,&py));
461:       PetscCall(VecPlaceArray(wr,px));
462:       PetscCall(VecPlaceArray(wi,py));
463:       PetscCall(VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL));
464:       PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn1));
465:       PetscCall(BVInsertVec(pep->V,i,wr));
466:       PetscCall(BVInsertVec(pep->V,i+1,wi));
467:       for (k=1;k<pep->nmat-1;k++) {
468:         PetscCall(VecResetArray(wr));
469:         PetscCall(VecResetArray(wi));
470:         PetscCall(VecPlaceArray(wr,px+k*pep->nloc));
471:         PetscCall(VecPlaceArray(wi,py+k*pep->nloc));
472:         PetscCall(VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL));
473:         PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn2));
474:         if (rn1>rn2) {
475:           PetscCall(BVInsertVec(pep->V,i,wr));
476:           PetscCall(BVInsertVec(pep->V,i+1,wi));
477:           rn1 = rn2;
478:         }
479:       }
480:       PetscCall(VecResetArray(wr));
481:       PetscCall(VecResetArray(wi));
482:       PetscCall(VecRestoreArrayRead(xr,&px));
483:       PetscCall(VecRestoreArrayRead(xi,&py));
484:       i++;
485:     } else   /* real eigenvalue */
486: #endif
487:     {
488:       PetscCall(VecGetArrayRead(xr,&px));
489:       PetscCall(VecPlaceArray(wr,px));
490:       PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
491:       PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn1));
492:       PetscCall(BVInsertVec(pep->V,i,wr));
493:       for (k=1;k<pep->nmat-1;k++) {
494:         PetscCall(VecResetArray(wr));
495:         PetscCall(VecPlaceArray(wr,px+k*pep->nloc));
496:         PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
497:         PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn2));
498:         if (rn1>rn2) {
499:           PetscCall(BVInsertVec(pep->V,i,wr));
500:           rn1 = rn2;
501:         }
502:       }
503:       PetscCall(VecResetArray(wr));
504:       PetscCall(VecRestoreArrayRead(xr,&px));
505:     }
506:   }
507:   PetscCall(VecDestroy(&wr));
508:   PetscCall(VecDestroy(&xr));
509: #if !defined(PETSC_USE_COMPLEX)
510:   PetscCall(VecDestroy(&wi));
511:   PetscCall(VecDestroy(&xi));
512: #endif
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: /*
517:    PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
518:    the first block.
519: */
520: static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
521: {
522:   PetscInt          i;
523:   const PetscScalar *px;
524:   Mat               A;
525:   Vec               xr,xi=NULL,w;
526: #if !defined(PETSC_USE_COMPLEX)
527:   PetscScalar       *ei=pep->eigi;
528: #endif

530:   PetscFunctionBegin;
531:   PetscCall(EPSGetOperators(eps,&A,NULL));
532:   PetscCall(MatCreateVecs(A,&xr,NULL));
533: #if !defined(PETSC_USE_COMPLEX)
534:   PetscCall(VecDuplicate(xr,&xi));
535: #endif
536:   PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
537:   for (i=0;i<pep->nconv;i++) {
538:     PetscCall(EPSGetEigenvector(eps,i,xr,xi));
539: #if !defined(PETSC_USE_COMPLEX)
540:     if (ei[i]!=0.0) {   /* complex conjugate pair */
541:       PetscCall(VecGetArrayRead(xr,&px));
542:       PetscCall(VecPlaceArray(w,px));
543:       PetscCall(BVInsertVec(pep->V,i,w));
544:       PetscCall(VecResetArray(w));
545:       PetscCall(VecRestoreArrayRead(xr,&px));
546:       PetscCall(VecGetArrayRead(xi,&px));
547:       PetscCall(VecPlaceArray(w,px));
548:       PetscCall(BVInsertVec(pep->V,i+1,w));
549:       PetscCall(VecResetArray(w));
550:       PetscCall(VecRestoreArrayRead(xi,&px));
551:       i++;
552:     } else   /* real eigenvalue */
553: #endif
554:     {
555:       PetscCall(VecGetArrayRead(xr,&px));
556:       PetscCall(VecPlaceArray(w,px));
557:       PetscCall(BVInsertVec(pep->V,i,w));
558:       PetscCall(VecResetArray(w));
559:       PetscCall(VecRestoreArrayRead(xr,&px));
560:     }
561:   }
562:   PetscCall(VecDestroy(&w));
563:   PetscCall(VecDestroy(&xr));
564: #if !defined(PETSC_USE_COMPLEX)
565:   PetscCall(VecDestroy(&xi));
566: #endif
567:   PetscFunctionReturn(PETSC_SUCCESS);
568: }

570: /*
571:    PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
572:    linear eigenproblem to the PEP object. The eigenvector of the generalized
573:    problem is supposed to be
574:                                z = [  x  ]
575:                                    [ l*x ]
576:    If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
577:    Finally, x is normalized so that ||x||_2 = 1.
578: */
579: static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
580: {
581:   PetscInt          i,offset;
582:   const PetscScalar *px;
583:   PetscScalar       *er=pep->eigr;
584:   Mat               A;
585:   Vec               xr,xi=NULL,w;
586: #if !defined(PETSC_USE_COMPLEX)
587:   PetscScalar       *ei=pep->eigi;
588: #endif

590:   PetscFunctionBegin;
591:   PetscCall(EPSGetOperators(eps,&A,NULL));
592:   PetscCall(MatCreateVecs(A,&xr,NULL));
593: #if !defined(PETSC_USE_COMPLEX)
594:   PetscCall(VecDuplicate(xr,&xi));
595: #endif
596:   PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
597:   for (i=0;i<pep->nconv;i++) {
598:     PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,xr,xi));
599:     if (SlepcAbsEigenvalue(er[i],ei[i])>1.0) offset = (pep->nmat-2)*pep->nloc;
600:     else offset = 0;
601: #if !defined(PETSC_USE_COMPLEX)
602:     if (ei[i]!=0.0) {   /* complex conjugate pair */
603:       PetscCall(VecGetArrayRead(xr,&px));
604:       PetscCall(VecPlaceArray(w,px+offset));
605:       PetscCall(BVInsertVec(pep->V,i,w));
606:       PetscCall(VecResetArray(w));
607:       PetscCall(VecRestoreArrayRead(xr,&px));
608:       PetscCall(VecGetArrayRead(xi,&px));
609:       PetscCall(VecPlaceArray(w,px+offset));
610:       PetscCall(BVInsertVec(pep->V,i+1,w));
611:       PetscCall(VecResetArray(w));
612:       PetscCall(VecRestoreArrayRead(xi,&px));
613:       i++;
614:     } else /* real eigenvalue */
615: #endif
616:     {
617:       PetscCall(VecGetArrayRead(xr,&px));
618:       PetscCall(VecPlaceArray(w,px+offset));
619:       PetscCall(BVInsertVec(pep->V,i,w));
620:       PetscCall(VecResetArray(w));
621:       PetscCall(VecRestoreArrayRead(xr,&px));
622:     }
623:   }
624:   PetscCall(VecDestroy(&w));
625:   PetscCall(VecDestroy(&xr));
626: #if !defined(PETSC_USE_COMPLEX)
627:   PetscCall(VecDestroy(&xi));
628: #endif
629:   PetscFunctionReturn(PETSC_SUCCESS);
630: }

632: static PetscErrorCode PEPExtractVectors_Linear(PEP pep)
633: {
634:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

636:   PetscFunctionBegin;
637:   switch (pep->extract) {
638:   case PEP_EXTRACT_NONE:
639:     PetscCall(PEPLinearExtract_None(pep,ctx->eps));
640:     break;
641:   case PEP_EXTRACT_NORM:
642:     PetscCall(PEPLinearExtract_Norm(pep,ctx->eps));
643:     break;
644:   case PEP_EXTRACT_RESIDUAL:
645:     PetscCall(PEPLinearExtract_Residual(pep,ctx->eps));
646:     break;
647:   case PEP_EXTRACT_STRUCTURED:
648:     SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Extraction not implemented in this solver");
649:   }
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: static PetscErrorCode PEPSolve_Linear(PEP pep)
654: {
655:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
656:   PetscScalar    sigma;
657:   PetscBool      flg;
658:   PetscInt       i;

660:   PetscFunctionBegin;
661:   PetscCall(EPSSolve(ctx->eps));
662:   PetscCall(EPSGetConverged(ctx->eps,&pep->nconv));
663:   PetscCall(EPSGetIterationNumber(ctx->eps,&pep->its));
664:   PetscCall(EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason));

666:   /* recover eigenvalues */
667:   for (i=0;i<pep->nconv;i++) {
668:     PetscCall(EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL));
669:     pep->eigr[i] *= pep->sfactor;
670:     pep->eigi[i] *= pep->sfactor;
671:   }

673:   /* restore target */
674:   PetscCall(EPSGetTarget(ctx->eps,&sigma));
675:   PetscCall(EPSSetTarget(ctx->eps,sigma*pep->sfactor));

677:   PetscCall(STGetTransform(pep->st,&flg));
678:   if (flg) PetscTryTypeMethod(pep,backtransform);
679:   if (pep->sfactor!=1.0) {
680:     /* Restore original values */
681:     for (i=0;i<pep->nmat;i++) {
682:       pep->pbc[pep->nmat+i] *= pep->sfactor;
683:       pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
684:     }
685:     if (!flg && !ctx->explicitmatrix) PetscCall(STScaleShift(pep->st,pep->sfactor));
686:   }
687:   if (ctx->explicitmatrix || !flg) PetscCall(RGPopScale(pep->rg));
688:   PetscFunctionReturn(PETSC_SUCCESS);
689: }

691: static PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
692: {
693:   PEP            pep = (PEP)ctx;

695:   PetscFunctionBegin;
696:   PetscCall(PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest));
697:   PetscFunctionReturn(PETSC_SUCCESS);
698: }

700: static PetscErrorCode PEPSetFromOptions_Linear(PEP pep,PetscOptionItems *PetscOptionsObject)
701: {
702:   PetscBool      set,val;
703:   PetscInt       k;
704:   PetscReal      array[2]={0,0};
705:   PetscBool      flg;
706:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

708:   PetscFunctionBegin;
709:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP Linear Options");

711:     k = 2;
712:     PetscCall(PetscOptionsRealArray("-pep_linear_linearization","Parameters of the linearization","PEPLinearSetLinearization",array,&k,&flg));
713:     if (flg) PetscCall(PEPLinearSetLinearization(pep,array[0],array[1]));

715:     PetscCall(PetscOptionsBool("-pep_linear_explicitmatrix","Use explicit matrix in linearization","PEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set));
716:     if (set) PetscCall(PEPLinearSetExplicitMatrix(pep,val));

718:   PetscOptionsHeadEnd();

720:   if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
721:   PetscCall(EPSSetFromOptions(ctx->eps));
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: static PetscErrorCode PEPLinearSetLinearization_Linear(PEP pep,PetscReal alpha,PetscReal beta)
726: {
727:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

729:   PetscFunctionBegin;
730:   PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
731:   ctx->alpha = alpha;
732:   ctx->beta  = beta;
733:   PetscFunctionReturn(PETSC_SUCCESS);
734: }

736: /*@
737:    PEPLinearSetLinearization - Set the coefficients that define
738:    the linearization of a quadratic eigenproblem.

740:    Logically Collective

742:    Input Parameters:
743: +  pep   - polynomial eigenvalue solver
744: .  alpha - first parameter of the linearization
745: -  beta  - second parameter of the linearization

747:    Options Database Key:
748: .  -pep_linear_linearization <alpha,beta> - Sets the coefficients

750:    Notes:
751:    Cannot pass zero for both alpha and beta. The default values are
752:    alpha=1 and beta=0.

754:    Level: advanced

756: .seealso: PEPLinearGetLinearization()
757: @*/
758: PetscErrorCode PEPLinearSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
759: {
760:   PetscFunctionBegin;
764:   PetscTryMethod(pep,"PEPLinearSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
765:   PetscFunctionReturn(PETSC_SUCCESS);
766: }

768: static PetscErrorCode PEPLinearGetLinearization_Linear(PEP pep,PetscReal *alpha,PetscReal *beta)
769: {
770:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

772:   PetscFunctionBegin;
773:   if (alpha) *alpha = ctx->alpha;
774:   if (beta)  *beta  = ctx->beta;
775:   PetscFunctionReturn(PETSC_SUCCESS);
776: }

778: /*@
779:    PEPLinearGetLinearization - Returns the coefficients that define
780:    the linearization of a quadratic eigenproblem.

782:    Not Collective

784:    Input Parameter:
785: .  pep  - polynomial eigenvalue solver

787:    Output Parameters:
788: +  alpha - the first parameter of the linearization
789: -  beta  - the second parameter of the linearization

791:    Level: advanced

793: .seealso: PEPLinearSetLinearization()
794: @*/
795: PetscErrorCode PEPLinearGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
796: {
797:   PetscFunctionBegin;
799:   PetscUseMethod(pep,"PEPLinearGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
800:   PetscFunctionReturn(PETSC_SUCCESS);
801: }

803: static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
804: {
805:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

807:   PetscFunctionBegin;
808:   if (ctx->explicitmatrix != explicitmatrix) {
809:     ctx->explicitmatrix = explicitmatrix;
810:     pep->state = PEP_STATE_INITIAL;
811:   }
812:   PetscFunctionReturn(PETSC_SUCCESS);
813: }

815: /*@
816:    PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
817:    linearization of the problem must be built explicitly.

819:    Logically Collective

821:    Input Parameters:
822: +  pep         - polynomial eigenvalue solver
823: -  explicitmat - boolean flag indicating if the matrices are built explicitly

825:    Options Database Key:
826: .  -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag

828:    Level: advanced

830: .seealso: PEPLinearGetExplicitMatrix()
831: @*/
832: PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmat)
833: {
834:   PetscFunctionBegin;
837:   PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmat));
838:   PetscFunctionReturn(PETSC_SUCCESS);
839: }

841: static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmat)
842: {
843:   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;

845:   PetscFunctionBegin;
846:   *explicitmat = ctx->explicitmatrix;
847:   PetscFunctionReturn(PETSC_SUCCESS);
848: }

850: /*@
851:    PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
852:    A and B for the linearization are built explicitly.

854:    Not Collective

856:    Input Parameter:
857: .  pep  - polynomial eigenvalue solver

859:    Output Parameter:
860: .  explicitmat - the mode flag

862:    Level: advanced

864: .seealso: PEPLinearSetExplicitMatrix()
865: @*/
866: PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmat)
867: {
868:   PetscFunctionBegin;
870:   PetscAssertPointer(explicitmat,2);
871:   PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmat));
872:   PetscFunctionReturn(PETSC_SUCCESS);
873: }

875: static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
876: {
877:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

879:   PetscFunctionBegin;
880:   PetscCall(PetscObjectReference((PetscObject)eps));
881:   PetscCall(EPSDestroy(&ctx->eps));
882:   ctx->eps     = eps;
883:   ctx->usereps = PETSC_TRUE;
884:   pep->state   = PEP_STATE_INITIAL;
885:   PetscFunctionReturn(PETSC_SUCCESS);
886: }

888: /*@
889:    PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
890:    polynomial eigenvalue solver.

892:    Collective

894:    Input Parameters:
895: +  pep - polynomial eigenvalue solver
896: -  eps - the eigensolver object

898:    Level: advanced

900: .seealso: PEPLinearGetEPS()
901: @*/
902: PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
903: {
904:   PetscFunctionBegin;
907:   PetscCheckSameComm(pep,1,eps,2);
908:   PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
909:   PetscFunctionReturn(PETSC_SUCCESS);
910: }

912: static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
913: {
914:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

916:   PetscFunctionBegin;
917:   if (!ctx->eps) {
918:     PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps));
919:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1));
920:     PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix));
921:     PetscCall(EPSAppendOptionsPrefix(ctx->eps,"pep_linear_"));
922:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)pep)->options));
923:     PetscCall(EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL));
924:   }
925:   *eps = ctx->eps;
926:   PetscFunctionReturn(PETSC_SUCCESS);
927: }

929: /*@
930:    PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
931:    to the polynomial eigenvalue solver.

933:    Collective

935:    Input Parameter:
936: .  pep - polynomial eigenvalue solver

938:    Output Parameter:
939: .  eps - the eigensolver object

941:    Level: advanced

943: .seealso: PEPLinearSetEPS()
944: @*/
945: PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
946: {
947:   PetscFunctionBegin;
949:   PetscAssertPointer(eps,2);
950:   PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
951:   PetscFunctionReturn(PETSC_SUCCESS);
952: }

954: static PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
955: {
956:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
957:   PetscBool      isascii;

959:   PetscFunctionBegin;
960:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
961:   if (isascii) {
962:     if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
963:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrices\n",ctx->explicitmatrix? "explicit": "implicit"));
964:     PetscCall(PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
965:     PetscCall(PetscViewerASCIIPushTab(viewer));
966:     PetscCall(EPSView(ctx->eps,viewer));
967:     PetscCall(PetscViewerASCIIPopTab(viewer));
968:   }
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

972: static PetscErrorCode PEPReset_Linear(PEP pep)
973: {
974:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

976:   PetscFunctionBegin;
977:   if (!ctx->eps) PetscCall(EPSReset(ctx->eps));
978:   PetscCall(MatDestroy(&ctx->A));
979:   PetscCall(MatDestroy(&ctx->B));
980:   PetscCall(VecDestroy(&ctx->w[0]));
981:   PetscCall(VecDestroy(&ctx->w[1]));
982:   PetscCall(VecDestroy(&ctx->w[2]));
983:   PetscCall(VecDestroy(&ctx->w[3]));
984:   PetscCall(VecDestroy(&ctx->w[4]));
985:   PetscCall(VecDestroy(&ctx->w[5]));
986:   PetscFunctionReturn(PETSC_SUCCESS);
987: }

989: static PetscErrorCode PEPDestroy_Linear(PEP pep)
990: {
991:   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;

993:   PetscFunctionBegin;
994:   PetscCall(EPSDestroy(&ctx->eps));
995:   PetscCall(PetscFree(pep->data));
996:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",NULL));
997:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",NULL));
998:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL));
999:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL));
1000:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL));
1001:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL));
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

1005: SLEPC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
1006: {
1007:   PEP_LINEAR     *ctx;

1009:   PetscFunctionBegin;
1010:   PetscCall(PetscNew(&ctx));
1011:   pep->data = (void*)ctx;

1013:   pep->lineariz       = PETSC_TRUE;
1014:   ctx->explicitmatrix = PETSC_FALSE;
1015:   ctx->alpha          = 1.0;
1016:   ctx->beta           = 0.0;

1018:   pep->ops->solve          = PEPSolve_Linear;
1019:   pep->ops->setup          = PEPSetUp_Linear;
1020:   pep->ops->setfromoptions = PEPSetFromOptions_Linear;
1021:   pep->ops->destroy        = PEPDestroy_Linear;
1022:   pep->ops->reset          = PEPReset_Linear;
1023:   pep->ops->view           = PEPView_Linear;
1024:   pep->ops->backtransform  = PEPBackTransform_Default;
1025:   pep->ops->computevectors = PEPComputeVectors_Default;
1026:   pep->ops->extractvectors = PEPExtractVectors_Linear;

1028:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",PEPLinearSetLinearization_Linear));
1029:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",PEPLinearGetLinearization_Linear));
1030:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear));
1031:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear));
1032:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear));
1033:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear));
1034:   PetscFunctionReturn(PETSC_SUCCESS);
1035: }