Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, 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 : Explicit linearization for polynomial eigenproblems
12 : */
13 :
14 : #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
15 : #include "linear.h"
16 :
17 998 : static PetscErrorCode MatMult_Linear_Shift(Mat M,Vec x,Vec y)
18 : {
19 998 : PEP_LINEAR *ctx;
20 998 : PEP pep;
21 998 : const PetscScalar *px;
22 998 : PetscScalar *py,a,sigma=0.0;
23 998 : PetscInt nmat,deg,i,m;
24 998 : Vec x1,x2,x3,y1,aux;
25 998 : PetscReal *ca,*cb,*cg;
26 998 : PetscBool flg;
27 :
28 998 : PetscFunctionBegin;
29 998 : PetscCall(MatShellGetContext(M,&ctx));
30 998 : pep = ctx->pep;
31 998 : PetscCall(STGetTransform(pep->st,&flg));
32 998 : if (!flg) PetscCall(STGetShift(pep->st,&sigma));
33 998 : nmat = pep->nmat;
34 998 : deg = nmat-1;
35 998 : m = pep->nloc;
36 998 : ca = pep->pbc;
37 998 : cb = pep->pbc+nmat;
38 998 : cg = pep->pbc+2*nmat;
39 998 : x1=ctx->w[0];x2=ctx->w[1];x3=ctx->w[2];y1=ctx->w[3];aux=ctx->w[4];
40 :
41 998 : PetscCall(VecSet(y,0.0));
42 998 : PetscCall(VecGetArrayRead(x,&px));
43 998 : PetscCall(VecGetArray(y,&py));
44 998 : a = 1.0;
45 :
46 : /* first block */
47 998 : PetscCall(VecPlaceArray(x2,px));
48 998 : PetscCall(VecPlaceArray(x3,px+m));
49 998 : PetscCall(VecPlaceArray(y1,py));
50 998 : PetscCall(VecAXPY(y1,cb[0]-sigma,x2));
51 998 : PetscCall(VecAXPY(y1,ca[0],x3));
52 998 : PetscCall(VecResetArray(x2));
53 998 : PetscCall(VecResetArray(x3));
54 998 : PetscCall(VecResetArray(y1));
55 :
56 : /* inner blocks */
57 1096 : for (i=1;i<deg-1;i++) {
58 98 : PetscCall(VecPlaceArray(x1,px+(i-1)*m));
59 98 : PetscCall(VecPlaceArray(x2,px+i*m));
60 98 : PetscCall(VecPlaceArray(x3,px+(i+1)*m));
61 98 : PetscCall(VecPlaceArray(y1,py+i*m));
62 98 : PetscCall(VecAXPY(y1,cg[i],x1));
63 98 : PetscCall(VecAXPY(y1,cb[i]-sigma,x2));
64 98 : PetscCall(VecAXPY(y1,ca[i],x3));
65 98 : PetscCall(VecResetArray(x1));
66 98 : PetscCall(VecResetArray(x2));
67 98 : PetscCall(VecResetArray(x3));
68 98 : PetscCall(VecResetArray(y1));
69 : }
70 :
71 : /* last block */
72 998 : PetscCall(VecPlaceArray(y1,py+(deg-1)*m));
73 3092 : for (i=0;i<deg;i++) {
74 2094 : PetscCall(VecPlaceArray(x1,px+i*m));
75 2094 : PetscCall(STMatMult(pep->st,i,x1,aux));
76 2094 : PetscCall(VecAXPY(y1,a,aux));
77 2094 : PetscCall(VecResetArray(x1));
78 2094 : a *= pep->sfactor;
79 : }
80 998 : PetscCall(VecCopy(y1,aux));
81 998 : PetscCall(STMatSolve(pep->st,aux,y1));
82 998 : PetscCall(VecScale(y1,-ca[deg-1]/a));
83 998 : PetscCall(VecPlaceArray(x1,px+(deg-2)*m));
84 998 : PetscCall(VecPlaceArray(x2,px+(deg-1)*m));
85 998 : PetscCall(VecAXPY(y1,cg[deg-1],x1));
86 998 : PetscCall(VecAXPY(y1,cb[deg-1]-sigma,x2));
87 998 : PetscCall(VecResetArray(x1));
88 998 : PetscCall(VecResetArray(x2));
89 998 : PetscCall(VecResetArray(y1));
90 :
91 998 : PetscCall(VecRestoreArrayRead(x,&px));
92 998 : PetscCall(VecRestoreArray(y,&py));
93 998 : PetscFunctionReturn(PETSC_SUCCESS);
94 : }
95 :
96 1360 : static PetscErrorCode MatMult_Linear_Sinvert(Mat M,Vec x,Vec y)
97 : {
98 1360 : PEP_LINEAR *ctx;
99 1360 : PEP pep;
100 1360 : const PetscScalar *px;
101 1360 : PetscScalar *py,a,sigma,t=1.0,tp=0.0,tt;
102 1360 : PetscInt nmat,deg,i,m;
103 1360 : Vec x1,y1,y2,y3,aux,aux2;
104 1360 : PetscReal *ca,*cb,*cg;
105 :
106 1360 : PetscFunctionBegin;
107 1360 : PetscCall(MatShellGetContext(M,&ctx));
108 1360 : pep = ctx->pep;
109 1360 : nmat = pep->nmat;
110 1360 : deg = nmat-1;
111 1360 : m = pep->nloc;
112 1360 : ca = pep->pbc;
113 1360 : cb = pep->pbc+nmat;
114 1360 : cg = pep->pbc+2*nmat;
115 1360 : 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 1360 : PetscCall(EPSGetTarget(ctx->eps,&sigma));
117 1360 : PetscCall(VecSet(y,0.0));
118 1360 : PetscCall(VecGetArrayRead(x,&px));
119 1360 : PetscCall(VecGetArray(y,&py));
120 1360 : a = pep->sfactor;
121 :
122 : /* first block */
123 1360 : PetscCall(VecPlaceArray(x1,px));
124 1360 : PetscCall(VecPlaceArray(y1,py+m));
125 1360 : PetscCall(VecCopy(x1,y1));
126 1360 : PetscCall(VecScale(y1,1.0/ca[0]));
127 1360 : PetscCall(VecResetArray(x1));
128 1360 : PetscCall(VecResetArray(y1));
129 :
130 : /* second block */
131 1360 : if (deg>2) {
132 497 : PetscCall(VecPlaceArray(x1,px+m));
133 497 : PetscCall(VecPlaceArray(y1,py+m));
134 497 : PetscCall(VecPlaceArray(y2,py+2*m));
135 497 : PetscCall(VecCopy(x1,y2));
136 497 : PetscCall(VecAXPY(y2,sigma-cb[1],y1));
137 497 : PetscCall(VecScale(y2,1.0/ca[1]));
138 497 : PetscCall(VecResetArray(x1));
139 497 : PetscCall(VecResetArray(y1));
140 497 : PetscCall(VecResetArray(y2));
141 : }
142 :
143 : /* inner blocks */
144 1857 : for (i=2;i<deg-1;i++) {
145 497 : PetscCall(VecPlaceArray(x1,px+i*m));
146 497 : PetscCall(VecPlaceArray(y1,py+(i-1)*m));
147 497 : PetscCall(VecPlaceArray(y2,py+i*m));
148 497 : PetscCall(VecPlaceArray(y3,py+(i+1)*m));
149 497 : PetscCall(VecCopy(x1,y3));
150 497 : PetscCall(VecAXPY(y3,sigma-cb[i],y2));
151 497 : PetscCall(VecAXPY(y3,-cg[i],y1));
152 497 : PetscCall(VecScale(y3,1.0/ca[i]));
153 497 : PetscCall(VecResetArray(x1));
154 497 : PetscCall(VecResetArray(y1));
155 497 : PetscCall(VecResetArray(y2));
156 497 : PetscCall(VecResetArray(y3));
157 : }
158 :
159 : /* last block */
160 1360 : PetscCall(VecPlaceArray(y1,py));
161 2354 : for (i=0;i<deg-2;i++) {
162 994 : PetscCall(VecPlaceArray(y2,py+(i+1)*m));
163 994 : PetscCall(STMatMult(pep->st,i+1,y2,aux));
164 994 : PetscCall(VecAXPY(y1,a,aux));
165 994 : PetscCall(VecResetArray(y2));
166 994 : a *= pep->sfactor;
167 : }
168 1360 : i = deg-2;
169 1360 : PetscCall(VecPlaceArray(y2,py+(i+1)*m));
170 1360 : PetscCall(VecPlaceArray(y3,py+i*m));
171 1360 : PetscCall(VecCopy(y2,aux2));
172 1360 : PetscCall(VecAXPY(aux2,cg[i+1]/ca[i+1],y3));
173 1360 : PetscCall(STMatMult(pep->st,i+1,aux2,aux));
174 1360 : PetscCall(VecAXPY(y1,a,aux));
175 1360 : PetscCall(VecResetArray(y2));
176 1360 : PetscCall(VecResetArray(y3));
177 1360 : a *= pep->sfactor;
178 1360 : i = deg-1;
179 1360 : PetscCall(VecPlaceArray(x1,px+i*m));
180 1360 : PetscCall(VecPlaceArray(y3,py+i*m));
181 1360 : PetscCall(VecCopy(x1,aux2));
182 1360 : PetscCall(VecAXPY(aux2,sigma-cb[i],y3));
183 1360 : PetscCall(VecScale(aux2,1.0/ca[i]));
184 1360 : PetscCall(STMatMult(pep->st,i+1,aux2,aux));
185 1360 : PetscCall(VecAXPY(y1,a,aux));
186 1360 : PetscCall(VecResetArray(x1));
187 1360 : PetscCall(VecResetArray(y3));
188 :
189 1360 : PetscCall(VecCopy(y1,aux));
190 1360 : PetscCall(STMatSolve(pep->st,aux,y1));
191 1360 : PetscCall(VecScale(y1,-1.0));
192 :
193 : /* final update */
194 3714 : for (i=1;i<deg;i++) {
195 2354 : PetscCall(VecPlaceArray(y2,py+i*m));
196 2354 : tt = t;
197 2354 : t = ((sigma-cb[i-1])*t-cg[i-1]*tp)/ca[i-1]; /* i-th basis polynomial */
198 2354 : tp = tt;
199 2354 : PetscCall(VecAXPY(y2,t,y1));
200 2354 : PetscCall(VecResetArray(y2));
201 : }
202 1360 : PetscCall(VecResetArray(y1));
203 :
204 1360 : PetscCall(VecRestoreArrayRead(x,&px));
205 1360 : PetscCall(VecRestoreArray(y,&py));
206 1360 : PetscFunctionReturn(PETSC_SUCCESS);
207 : }
208 :
209 41949 : static PetscErrorCode BackTransform_Linear(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
210 : {
211 41949 : PEP_LINEAR *ctx;
212 41949 : ST stctx;
213 :
214 41949 : PetscFunctionBegin;
215 41949 : PetscCall(STShellGetContext(st,&ctx));
216 41949 : PetscCall(PEPGetST(ctx->pep,&stctx));
217 41949 : PetscCall(STBackTransform(stctx,n,eigr,eigi));
218 41949 : PetscFunctionReturn(PETSC_SUCCESS);
219 : }
220 :
221 : /*
222 : Dummy backtransform operation
223 : */
224 172 : static PetscErrorCode BackTransform_Skip(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
225 : {
226 172 : PetscFunctionBegin;
227 172 : PetscFunctionReturn(PETSC_SUCCESS);
228 : }
229 :
230 2358 : static PetscErrorCode Apply_Linear(ST st,Vec x,Vec y)
231 : {
232 2358 : PEP_LINEAR *ctx;
233 :
234 2358 : PetscFunctionBegin;
235 2358 : PetscCall(STShellGetContext(st,&ctx));
236 2358 : PetscCall(MatMult(ctx->A,x,y));
237 2358 : PetscFunctionReturn(PETSC_SUCCESS);
238 : }
239 :
240 35 : static PetscErrorCode PEPSetUp_Linear(PEP pep)
241 : {
242 35 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
243 35 : ST st;
244 35 : PetscInt i=0,deg=pep->nmat-1;
245 35 : EPSWhich which = EPS_LARGEST_MAGNITUDE;
246 35 : EPSProblemType ptype;
247 35 : PetscBool trackall,istrivial,transf,sinv,ks;
248 35 : PetscScalar sigma,*epsarray,*peparray;
249 35 : Vec veps,w=NULL;
250 : /* function tables */
251 35 : 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 : };
256 :
257 35 : PetscFunctionBegin;
258 35 : PEPCheckShiftSinvert(pep);
259 35 : PEPCheckUnsupported(pep,PEP_FEATURE_STOPPING);
260 35 : PEPCheckIgnored(pep,PEP_FEATURE_CONVERGENCE);
261 35 : PetscCall(STGetTransform(pep->st,&transf));
262 35 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
263 35 : if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
264 35 : PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
265 35 : PetscCall(STSetUp(pep->st));
266 35 : if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
267 35 : PetscCall(EPSGetST(ctx->eps,&st));
268 35 : if (!transf && !ctx->usereps) PetscCall(EPSSetTarget(ctx->eps,pep->target));
269 35 : if (sinv && !transf && !ctx->usereps) PetscCall(STSetDefaultShift(st,pep->target));
270 : /* compute scale factor if not set by user */
271 35 : PetscCall(PEPComputeScaleFactor(pep));
272 :
273 35 : if (ctx->explicitmatrix) {
274 11 : PEPCheckQuadraticCondition(pep,PETSC_TRUE," (with explicit matrix)");
275 11 : PEPCheckUnsupportedCondition(pep,PEP_FEATURE_NONMONOMIAL,PETSC_TRUE," (with explicit matrix)");
276 11 : PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
277 11 : 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 11 : if (sinv && !transf) PetscCall(STSetType(st,STSINVERT));
279 11 : PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
280 11 : PetscCall(STGetMatrixTransformed(pep->st,0,&ctx->K));
281 11 : PetscCall(STGetMatrixTransformed(pep->st,1,&ctx->C));
282 11 : PetscCall(STGetMatrixTransformed(pep->st,2,&ctx->M));
283 11 : ctx->sfactor = pep->sfactor;
284 11 : ctx->dsfactor = pep->dsfactor;
285 :
286 11 : PetscCall(MatDestroy(&ctx->A));
287 11 : PetscCall(MatDestroy(&ctx->B));
288 11 : PetscCall(VecDestroy(&ctx->w[0]));
289 11 : PetscCall(VecDestroy(&ctx->w[1]));
290 11 : PetscCall(VecDestroy(&ctx->w[2]));
291 11 : PetscCall(VecDestroy(&ctx->w[3]));
292 :
293 11 : switch (pep->problem_type) {
294 : case PEP_GENERAL: i = 0; break;
295 3 : case PEP_HERMITIAN:
296 3 : case PEP_HYPERBOLIC: i = 1; break;
297 2 : case PEP_GYROSCOPIC: i = 2; break;
298 : }
299 :
300 11 : PetscCall((*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A));
301 11 : PetscCall((*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B));
302 :
303 : } else { /* implicit matrix */
304 24 : 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 24 : if (!((PetscObject)ctx->eps)->type_name) PetscCall(EPSSetType(ctx->eps,EPSKRYLOVSCHUR));
306 : else {
307 23 : PetscCall(PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks));
308 23 : PetscCheck(ks,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
309 : }
310 24 : 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 24 : PetscCall(STSetType(st,STSHELL));
312 24 : PetscCall(STShellSetContext(st,ctx));
313 24 : if (!transf) PetscCall(STShellSetBackTransform(st,BackTransform_Linear));
314 1 : else PetscCall(STShellSetBackTransform(st,BackTransform_Skip));
315 24 : PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]));
316 24 : PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]));
317 24 : PetscCall(MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]));
318 24 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A));
319 24 : if (sinv && !transf) PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert));
320 13 : else PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift));
321 24 : PetscCall(STShellSetApply(st,Apply_Linear));
322 24 : ctx->pep = pep;
323 :
324 24 : PetscCall(PEPBasisCoefficients(pep,pep->pbc));
325 24 : if (!transf) {
326 23 : PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
327 23 : if (sinv) PetscCall(PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL));
328 : else {
329 38 : for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
330 12 : pep->solvematcoeffs[deg] = 1.0;
331 : }
332 23 : PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
333 23 : PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
334 : }
335 24 : if (pep->sfactor!=1.0) {
336 0 : for (i=0;i<pep->nmat;i++) {
337 0 : pep->pbc[pep->nmat+i] /= pep->sfactor;
338 0 : pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
339 : }
340 : }
341 : }
342 :
343 35 : PetscCall(EPSSetOperators(ctx->eps,ctx->A,ctx->B));
344 35 : PetscCall(EPSGetProblemType(ctx->eps,&ptype));
345 35 : if (!ptype) {
346 33 : if (ctx->explicitmatrix) PetscCall(EPSSetProblemType(ctx->eps,EPS_GNHEP));
347 23 : else PetscCall(EPSSetProblemType(ctx->eps,EPS_NHEP));
348 : }
349 35 : if (!ctx->usereps) {
350 34 : if (transf) which = EPS_LARGEST_MAGNITUDE;
351 : else {
352 33 : switch (pep->which) {
353 : case PEP_LARGEST_MAGNITUDE: which = EPS_LARGEST_MAGNITUDE; break;
354 0 : case PEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
355 0 : case PEP_LARGEST_REAL: which = EPS_LARGEST_REAL; break;
356 0 : case PEP_SMALLEST_REAL: which = EPS_SMALLEST_REAL; break;
357 0 : case PEP_LARGEST_IMAGINARY: which = EPS_LARGEST_IMAGINARY; break;
358 0 : case PEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
359 15 : case PEP_TARGET_MAGNITUDE: which = EPS_TARGET_MAGNITUDE; break;
360 0 : case PEP_TARGET_REAL: which = EPS_TARGET_REAL; break;
361 0 : case PEP_TARGET_IMAGINARY: which = EPS_TARGET_IMAGINARY; break;
362 0 : case PEP_ALL: which = EPS_ALL; break;
363 1 : case PEP_WHICH_USER: which = EPS_WHICH_USER;
364 1 : PetscCall(EPSSetEigenvalueComparison(ctx->eps,pep->sc->comparison,pep->sc->comparisonctx));
365 : break;
366 : }
367 34 : }
368 34 : PetscCall(EPSSetWhichEigenpairs(ctx->eps,which));
369 :
370 34 : PetscCall(EPSSetDimensions(ctx->eps,pep->nev,pep->ncv,pep->mpd));
371 49 : PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(pep->tol),pep->max_it));
372 : }
373 35 : PetscCall(RGIsTrivial(pep->rg,&istrivial));
374 35 : if (!istrivial) {
375 1 : PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
376 1 : PetscCall(EPSSetRG(ctx->eps,pep->rg));
377 : }
378 : /* Transfer the trackall option from pep to eps */
379 35 : PetscCall(PEPGetTrackAll(pep,&trackall));
380 35 : PetscCall(EPSSetTrackAll(ctx->eps,trackall));
381 :
382 : /* temporary change of target */
383 35 : if (pep->sfactor!=1.0) {
384 1 : PetscCall(EPSGetTarget(ctx->eps,&sigma));
385 1 : PetscCall(EPSSetTarget(ctx->eps,sigma/pep->sfactor));
386 : }
387 :
388 : /* process initial vector */
389 35 : if (pep->nini<0) {
390 2 : PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps));
391 2 : PetscCall(VecGetArray(veps,&epsarray));
392 6 : for (i=0;i<deg;i++) {
393 4 : if (i<-pep->nini) {
394 4 : PetscCall(VecGetArray(pep->IS[i],&peparray));
395 4 : PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
396 4 : PetscCall(VecRestoreArray(pep->IS[i],&peparray));
397 : } else {
398 0 : if (!w) PetscCall(VecDuplicate(pep->IS[0],&w));
399 0 : PetscCall(VecSetRandom(w,NULL));
400 0 : PetscCall(VecGetArray(w,&peparray));
401 0 : PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
402 4 : PetscCall(VecRestoreArray(w,&peparray));
403 : }
404 : }
405 2 : PetscCall(VecRestoreArray(veps,&epsarray));
406 2 : PetscCall(EPSSetInitialSpace(ctx->eps,1,&veps));
407 2 : PetscCall(VecDestroy(&veps));
408 2 : PetscCall(VecDestroy(&w));
409 2 : PetscCall(SlepcBasisDestroy_Private(&pep->nini,&pep->IS));
410 : }
411 :
412 35 : PetscCall(EPSSetUp(ctx->eps));
413 35 : PetscCall(EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd));
414 35 : PetscCall(EPSGetTolerances(ctx->eps,NULL,&pep->max_it));
415 35 : PetscCall(PEPAllocateSolution(pep,0));
416 35 : PetscFunctionReturn(PETSC_SUCCESS);
417 : }
418 :
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 2 : static PetscErrorCode PEPLinearExtract_Residual(PEP pep,EPS eps)
430 : {
431 2 : PetscInt i,k;
432 2 : const PetscScalar *px;
433 2 : PetscScalar *er=pep->eigr,*ei=pep->eigi;
434 2 : PetscReal rn1,rn2;
435 2 : Vec xr,xi=NULL,wr;
436 2 : Mat A;
437 : #if !defined(PETSC_USE_COMPLEX)
438 : Vec wi;
439 : const PetscScalar *py;
440 : #endif
441 :
442 2 : PetscFunctionBegin;
443 : #if defined(PETSC_USE_COMPLEX)
444 2 : PetscCall(PEPSetWorkVecs(pep,2));
445 : #else
446 : PetscCall(PEPSetWorkVecs(pep,4));
447 : #endif
448 2 : PetscCall(EPSGetOperators(eps,&A,NULL));
449 2 : PetscCall(MatCreateVecs(A,&xr,NULL));
450 2 : 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 19 : for (i=0;i<pep->nconv;i++) {
456 17 : 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 17 : PetscCall(VecGetArrayRead(xr,&px));
489 17 : PetscCall(VecPlaceArray(wr,px));
490 17 : PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
491 17 : PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn1));
492 17 : PetscCall(BVInsertVec(pep->V,i,wr));
493 34 : for (k=1;k<pep->nmat-1;k++) {
494 17 : PetscCall(VecResetArray(wr));
495 17 : PetscCall(VecPlaceArray(wr,px+k*pep->nloc));
496 17 : PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
497 17 : PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn2));
498 17 : if (rn1>rn2) {
499 13 : PetscCall(BVInsertVec(pep->V,i,wr));
500 13 : rn1 = rn2;
501 : }
502 : }
503 17 : PetscCall(VecResetArray(wr));
504 17 : PetscCall(VecRestoreArrayRead(xr,&px));
505 : }
506 : }
507 2 : PetscCall(VecDestroy(&wr));
508 2 : PetscCall(VecDestroy(&xr));
509 : #if !defined(PETSC_USE_COMPLEX)
510 : PetscCall(VecDestroy(&wi));
511 : PetscCall(VecDestroy(&xi));
512 : #endif
513 2 : PetscFunctionReturn(PETSC_SUCCESS);
514 : }
515 :
516 : /*
517 : PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
518 : the first block.
519 : */
520 2 : static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
521 : {
522 2 : PetscInt i;
523 2 : const PetscScalar *px;
524 2 : Mat A;
525 2 : Vec xr,xi=NULL,w;
526 : #if !defined(PETSC_USE_COMPLEX)
527 : PetscScalar *ei=pep->eigi;
528 : #endif
529 :
530 2 : PetscFunctionBegin;
531 2 : PetscCall(EPSGetOperators(eps,&A,NULL));
532 2 : PetscCall(MatCreateVecs(A,&xr,NULL));
533 : #if !defined(PETSC_USE_COMPLEX)
534 : PetscCall(VecDuplicate(xr,&xi));
535 : #endif
536 2 : PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
537 19 : for (i=0;i<pep->nconv;i++) {
538 17 : 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 17 : PetscCall(VecGetArrayRead(xr,&px));
556 17 : PetscCall(VecPlaceArray(w,px));
557 17 : PetscCall(BVInsertVec(pep->V,i,w));
558 17 : PetscCall(VecResetArray(w));
559 17 : PetscCall(VecRestoreArrayRead(xr,&px));
560 : }
561 : }
562 2 : PetscCall(VecDestroy(&w));
563 2 : PetscCall(VecDestroy(&xr));
564 : #if !defined(PETSC_USE_COMPLEX)
565 : PetscCall(VecDestroy(&xi));
566 : #endif
567 2 : PetscFunctionReturn(PETSC_SUCCESS);
568 : }
569 :
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 31 : static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
580 : {
581 31 : PetscInt i,offset;
582 31 : const PetscScalar *px;
583 31 : PetscScalar *er=pep->eigr;
584 31 : Mat A;
585 31 : Vec xr,xi=NULL,w;
586 : #if !defined(PETSC_USE_COMPLEX)
587 : PetscScalar *ei=pep->eigi;
588 : #endif
589 :
590 31 : PetscFunctionBegin;
591 31 : PetscCall(EPSGetOperators(eps,&A,NULL));
592 31 : PetscCall(MatCreateVecs(A,&xr,NULL));
593 : #if !defined(PETSC_USE_COMPLEX)
594 : PetscCall(VecDuplicate(xr,&xi));
595 : #endif
596 31 : PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
597 187 : for (i=0;i<pep->nconv;i++) {
598 156 : PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,xr,xi));
599 156 : 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 156 : PetscCall(VecGetArrayRead(xr,&px));
618 156 : PetscCall(VecPlaceArray(w,px+offset));
619 156 : PetscCall(BVInsertVec(pep->V,i,w));
620 156 : PetscCall(VecResetArray(w));
621 156 : PetscCall(VecRestoreArrayRead(xr,&px));
622 : }
623 : }
624 31 : PetscCall(VecDestroy(&w));
625 31 : PetscCall(VecDestroy(&xr));
626 : #if !defined(PETSC_USE_COMPLEX)
627 : PetscCall(VecDestroy(&xi));
628 : #endif
629 31 : PetscFunctionReturn(PETSC_SUCCESS);
630 : }
631 :
632 35 : static PetscErrorCode PEPExtractVectors_Linear(PEP pep)
633 : {
634 35 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
635 :
636 35 : PetscFunctionBegin;
637 35 : switch (pep->extract) {
638 2 : case PEP_EXTRACT_NONE:
639 2 : PetscCall(PEPLinearExtract_None(pep,ctx->eps));
640 : break;
641 31 : case PEP_EXTRACT_NORM:
642 31 : PetscCall(PEPLinearExtract_Norm(pep,ctx->eps));
643 : break;
644 2 : case PEP_EXTRACT_RESIDUAL:
645 2 : PetscCall(PEPLinearExtract_Residual(pep,ctx->eps));
646 : break;
647 0 : case PEP_EXTRACT_STRUCTURED:
648 0 : SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Extraction not implemented in this solver");
649 : }
650 35 : PetscFunctionReturn(PETSC_SUCCESS);
651 : }
652 :
653 35 : static PetscErrorCode PEPSolve_Linear(PEP pep)
654 : {
655 35 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
656 35 : PetscScalar sigma;
657 35 : PetscBool flg;
658 35 : PetscInt i;
659 :
660 35 : PetscFunctionBegin;
661 35 : PetscCall(EPSSolve(ctx->eps));
662 35 : PetscCall(EPSGetConverged(ctx->eps,&pep->nconv));
663 35 : PetscCall(EPSGetIterationNumber(ctx->eps,&pep->its));
664 35 : PetscCall(EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason));
665 :
666 : /* recover eigenvalues */
667 225 : for (i=0;i<pep->nconv;i++) {
668 190 : PetscCall(EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL));
669 190 : pep->eigr[i] *= pep->sfactor;
670 190 : pep->eigi[i] *= pep->sfactor;
671 : }
672 :
673 : /* restore target */
674 35 : PetscCall(EPSGetTarget(ctx->eps,&sigma));
675 35 : PetscCall(EPSSetTarget(ctx->eps,sigma*pep->sfactor));
676 :
677 35 : PetscCall(STGetTransform(pep->st,&flg));
678 35 : if (flg) PetscTryTypeMethod(pep,backtransform);
679 35 : if (pep->sfactor!=1.0) {
680 : /* Restore original values */
681 4 : for (i=0;i<pep->nmat;i++) {
682 3 : pep->pbc[pep->nmat+i] *= pep->sfactor;
683 3 : pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
684 : }
685 1 : if (!flg && !ctx->explicitmatrix) PetscCall(STScaleShift(pep->st,pep->sfactor));
686 : }
687 35 : if (ctx->explicitmatrix || !flg) PetscCall(RGPopScale(pep->rg));
688 35 : PetscFunctionReturn(PETSC_SUCCESS);
689 : }
690 :
691 466 : static PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
692 : {
693 466 : PEP pep = (PEP)ctx;
694 :
695 466 : PetscFunctionBegin;
696 466 : PetscCall(PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest));
697 466 : PetscFunctionReturn(PETSC_SUCCESS);
698 : }
699 :
700 32 : static PetscErrorCode PEPSetFromOptions_Linear(PEP pep,PetscOptionItems *PetscOptionsObject)
701 : {
702 32 : PetscBool set,val;
703 32 : PetscInt k;
704 32 : PetscReal array[2]={0,0};
705 32 : PetscBool flg;
706 32 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
707 :
708 32 : PetscFunctionBegin;
709 32 : PetscOptionsHeadBegin(PetscOptionsObject,"PEP Linear Options");
710 :
711 32 : k = 2;
712 32 : PetscCall(PetscOptionsRealArray("-pep_linear_linearization","Parameters of the linearization","PEPLinearSetLinearization",array,&k,&flg));
713 32 : if (flg) PetscCall(PEPLinearSetLinearization(pep,array[0],array[1]));
714 :
715 32 : PetscCall(PetscOptionsBool("-pep_linear_explicitmatrix","Use explicit matrix in linearization","PEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set));
716 32 : if (set) PetscCall(PEPLinearSetExplicitMatrix(pep,val));
717 :
718 32 : PetscOptionsHeadEnd();
719 :
720 32 : if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
721 32 : PetscCall(EPSSetFromOptions(ctx->eps));
722 32 : PetscFunctionReturn(PETSC_SUCCESS);
723 : }
724 :
725 4 : static PetscErrorCode PEPLinearSetLinearization_Linear(PEP pep,PetscReal alpha,PetscReal beta)
726 : {
727 4 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
728 :
729 4 : PetscFunctionBegin;
730 4 : PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
731 4 : ctx->alpha = alpha;
732 4 : ctx->beta = beta;
733 4 : PetscFunctionReturn(PETSC_SUCCESS);
734 : }
735 :
736 : /*@
737 : PEPLinearSetLinearization - Set the coefficients that define
738 : the linearization of a quadratic eigenproblem.
739 :
740 : Logically Collective
741 :
742 : Input Parameters:
743 : + pep - polynomial eigenvalue solver
744 : . alpha - first parameter of the linearization
745 : - beta - second parameter of the linearization
746 :
747 : Options Database Key:
748 : . -pep_linear_linearization <alpha,beta> - Sets the coefficients
749 :
750 : Notes:
751 : Cannot pass zero for both alpha and beta. The default values are
752 : alpha=1 and beta=0.
753 :
754 : Level: advanced
755 :
756 : .seealso: PEPLinearGetLinearization()
757 : @*/
758 4 : PetscErrorCode PEPLinearSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
759 : {
760 4 : PetscFunctionBegin;
761 4 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
762 12 : PetscValidLogicalCollectiveReal(pep,alpha,2);
763 12 : PetscValidLogicalCollectiveReal(pep,beta,3);
764 4 : PetscTryMethod(pep,"PEPLinearSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
765 4 : PetscFunctionReturn(PETSC_SUCCESS);
766 : }
767 :
768 1 : static PetscErrorCode PEPLinearGetLinearization_Linear(PEP pep,PetscReal *alpha,PetscReal *beta)
769 : {
770 1 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
771 :
772 1 : PetscFunctionBegin;
773 1 : if (alpha) *alpha = ctx->alpha;
774 1 : if (beta) *beta = ctx->beta;
775 1 : PetscFunctionReturn(PETSC_SUCCESS);
776 : }
777 :
778 : /*@
779 : PEPLinearGetLinearization - Returns the coefficients that define
780 : the linearization of a quadratic eigenproblem.
781 :
782 : Not Collective
783 :
784 : Input Parameter:
785 : . pep - polynomial eigenvalue solver
786 :
787 : Output Parameters:
788 : + alpha - the first parameter of the linearization
789 : - beta - the second parameter of the linearization
790 :
791 : Level: advanced
792 :
793 : .seealso: PEPLinearSetLinearization()
794 : @*/
795 1 : PetscErrorCode PEPLinearGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
796 : {
797 1 : PetscFunctionBegin;
798 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
799 1 : PetscUseMethod(pep,"PEPLinearGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
800 1 : PetscFunctionReturn(PETSC_SUCCESS);
801 : }
802 :
803 11 : static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
804 : {
805 11 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
806 :
807 11 : PetscFunctionBegin;
808 11 : if (ctx->explicitmatrix != explicitmatrix) {
809 11 : ctx->explicitmatrix = explicitmatrix;
810 11 : pep->state = PEP_STATE_INITIAL;
811 : }
812 11 : PetscFunctionReturn(PETSC_SUCCESS);
813 : }
814 :
815 : /*@
816 : PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
817 : linearization of the problem must be built explicitly.
818 :
819 : Logically Collective
820 :
821 : Input Parameters:
822 : + pep - polynomial eigenvalue solver
823 : - explicitmat - boolean flag indicating if the matrices are built explicitly
824 :
825 : Options Database Key:
826 : . -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag
827 :
828 : Level: advanced
829 :
830 : .seealso: PEPLinearGetExplicitMatrix()
831 : @*/
832 11 : PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmat)
833 : {
834 11 : PetscFunctionBegin;
835 11 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
836 33 : PetscValidLogicalCollectiveBool(pep,explicitmat,2);
837 11 : PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmat));
838 11 : PetscFunctionReturn(PETSC_SUCCESS);
839 : }
840 :
841 1 : static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmat)
842 : {
843 1 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
844 :
845 1 : PetscFunctionBegin;
846 1 : *explicitmat = ctx->explicitmatrix;
847 1 : PetscFunctionReturn(PETSC_SUCCESS);
848 : }
849 :
850 : /*@
851 : PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
852 : A and B for the linearization are built explicitly.
853 :
854 : Not Collective
855 :
856 : Input Parameter:
857 : . pep - polynomial eigenvalue solver
858 :
859 : Output Parameter:
860 : . explicitmat - the mode flag
861 :
862 : Level: advanced
863 :
864 : .seealso: PEPLinearSetExplicitMatrix()
865 : @*/
866 1 : PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmat)
867 : {
868 1 : PetscFunctionBegin;
869 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
870 1 : PetscAssertPointer(explicitmat,2);
871 1 : PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmat));
872 1 : PetscFunctionReturn(PETSC_SUCCESS);
873 : }
874 :
875 1 : static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
876 : {
877 1 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
878 :
879 1 : PetscFunctionBegin;
880 1 : PetscCall(PetscObjectReference((PetscObject)eps));
881 1 : PetscCall(EPSDestroy(&ctx->eps));
882 1 : ctx->eps = eps;
883 1 : ctx->usereps = PETSC_TRUE;
884 1 : pep->state = PEP_STATE_INITIAL;
885 1 : PetscFunctionReturn(PETSC_SUCCESS);
886 : }
887 :
888 : /*@
889 : PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
890 : polynomial eigenvalue solver.
891 :
892 : Collective
893 :
894 : Input Parameters:
895 : + pep - polynomial eigenvalue solver
896 : - eps - the eigensolver object
897 :
898 : Level: advanced
899 :
900 : .seealso: PEPLinearGetEPS()
901 : @*/
902 1 : PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
903 : {
904 1 : PetscFunctionBegin;
905 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
906 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
907 1 : PetscCheckSameComm(pep,1,eps,2);
908 1 : PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
909 1 : PetscFunctionReturn(PETSC_SUCCESS);
910 : }
911 :
912 33 : static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
913 : {
914 33 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
915 :
916 33 : PetscFunctionBegin;
917 33 : if (!ctx->eps) {
918 33 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps));
919 33 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1));
920 33 : PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix));
921 33 : PetscCall(EPSAppendOptionsPrefix(ctx->eps,"pep_linear_"));
922 33 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)pep)->options));
923 33 : PetscCall(EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL));
924 : }
925 33 : *eps = ctx->eps;
926 33 : PetscFunctionReturn(PETSC_SUCCESS);
927 : }
928 :
929 : /*@
930 : PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
931 : to the polynomial eigenvalue solver.
932 :
933 : Collective
934 :
935 : Input Parameter:
936 : . pep - polynomial eigenvalue solver
937 :
938 : Output Parameter:
939 : . eps - the eigensolver object
940 :
941 : Level: advanced
942 :
943 : .seealso: PEPLinearSetEPS()
944 : @*/
945 33 : PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
946 : {
947 33 : PetscFunctionBegin;
948 33 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
949 33 : PetscAssertPointer(eps,2);
950 33 : PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
951 33 : PetscFunctionReturn(PETSC_SUCCESS);
952 : }
953 :
954 0 : static PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
955 : {
956 0 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
957 0 : PetscBool isascii;
958 :
959 0 : PetscFunctionBegin;
960 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
961 0 : if (isascii) {
962 0 : if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
963 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %s matrices\n",ctx->explicitmatrix? "explicit": "implicit"));
964 0 : PetscCall(PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
965 0 : PetscCall(PetscViewerASCIIPushTab(viewer));
966 0 : PetscCall(EPSView(ctx->eps,viewer));
967 0 : PetscCall(PetscViewerASCIIPopTab(viewer));
968 : }
969 0 : PetscFunctionReturn(PETSC_SUCCESS);
970 : }
971 :
972 35 : static PetscErrorCode PEPReset_Linear(PEP pep)
973 : {
974 35 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
975 :
976 35 : PetscFunctionBegin;
977 35 : if (!ctx->eps) PetscCall(EPSReset(ctx->eps));
978 35 : PetscCall(MatDestroy(&ctx->A));
979 35 : PetscCall(MatDestroy(&ctx->B));
980 35 : PetscCall(VecDestroy(&ctx->w[0]));
981 35 : PetscCall(VecDestroy(&ctx->w[1]));
982 35 : PetscCall(VecDestroy(&ctx->w[2]));
983 35 : PetscCall(VecDestroy(&ctx->w[3]));
984 35 : PetscCall(VecDestroy(&ctx->w[4]));
985 35 : PetscCall(VecDestroy(&ctx->w[5]));
986 35 : PetscFunctionReturn(PETSC_SUCCESS);
987 : }
988 :
989 34 : static PetscErrorCode PEPDestroy_Linear(PEP pep)
990 : {
991 34 : PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
992 :
993 34 : PetscFunctionBegin;
994 34 : PetscCall(EPSDestroy(&ctx->eps));
995 34 : PetscCall(PetscFree(pep->data));
996 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",NULL));
997 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",NULL));
998 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL));
999 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL));
1000 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL));
1001 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL));
1002 34 : PetscFunctionReturn(PETSC_SUCCESS);
1003 : }
1004 :
1005 34 : SLEPC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
1006 : {
1007 34 : PEP_LINEAR *ctx;
1008 :
1009 34 : PetscFunctionBegin;
1010 34 : PetscCall(PetscNew(&ctx));
1011 34 : pep->data = (void*)ctx;
1012 :
1013 34 : pep->lineariz = PETSC_TRUE;
1014 34 : ctx->explicitmatrix = PETSC_FALSE;
1015 34 : ctx->alpha = 1.0;
1016 34 : ctx->beta = 0.0;
1017 :
1018 34 : pep->ops->solve = PEPSolve_Linear;
1019 34 : pep->ops->setup = PEPSetUp_Linear;
1020 34 : pep->ops->setfromoptions = PEPSetFromOptions_Linear;
1021 34 : pep->ops->destroy = PEPDestroy_Linear;
1022 34 : pep->ops->reset = PEPReset_Linear;
1023 34 : pep->ops->view = PEPView_Linear;
1024 34 : pep->ops->backtransform = PEPBackTransform_Default;
1025 34 : pep->ops->computevectors = PEPComputeVectors_Default;
1026 34 : pep->ops->extractvectors = PEPExtractVectors_Linear;
1027 :
1028 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",PEPLinearSetLinearization_Linear));
1029 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",PEPLinearGetLinearization_Linear));
1030 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear));
1031 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear));
1032 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear));
1033 34 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear));
1034 34 : PetscFunctionReturn(PETSC_SUCCESS);
1035 : }
|
