Actual source code: linear.c
slepc-3.21.1 2024-04-26
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: }