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 : SLEPc polynomial eigensolver: "jd"
12 :
13 : Method: Jacobi-Davidson
14 :
15 : Algorithm:
16 :
17 : Jacobi-Davidson for polynomial eigenvalue problems.
18 :
19 : References:
20 :
21 : [1] C. Campos and J.E. Roman, "A polynomial Jacobi-Davidson solver
22 : with support for non-monomial bases and deflation", BIT Numer.
23 : Math. 60:295-318, 2020.
24 :
25 : [2] G.L.G. Sleijpen et al., "Jacobi-Davidson type methods for
26 : generalized eigenproblems and polynomial eigenproblems", BIT
27 : 36(3):595-633, 1996.
28 :
29 : [3] Feng-Nan Hwang, Zih-Hao Wei, Tsung-Ming Huang, Weichung Wang,
30 : "A Parallel Additive Schwarz Preconditioned Jacobi-Davidson
31 : Algorithm for Polynomial Eigenvalue Problems in Quantum Dot
32 : Simulation", J. Comput. Phys. 229(8):2932-2947, 2010.
33 : */
34 :
35 : #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
36 : #include <slepcblaslapack.h>
37 :
38 : static PetscBool cited = PETSC_FALSE;
39 : static const char citation[] =
40 : "@Article{slepc-slice-qep,\n"
41 : " author = \"C. Campos and J. E. Roman\",\n"
42 : " title = \"A polynomial {Jacobi-Davidson} solver with support for non-monomial bases and deflation\",\n"
43 : " journal = \"{BIT} Numer. Math.\",\n"
44 : " volume = \"60\",\n"
45 : " pages = \"295--318\",\n"
46 : " year = \"2020,\"\n"
47 : " doi = \"https://doi.org/10.1007/s10543-019-00778-z\"\n"
48 : "}\n";
49 :
50 : typedef struct {
51 : PetscReal keep; /* restart parameter */
52 : PetscReal fix; /* fix parameter */
53 : PetscBool reusepc; /* flag indicating whether pc is rebuilt or not */
54 : BV V; /* work basis vectors to store the search space */
55 : BV W; /* work basis vectors to store the test space */
56 : BV *TV; /* work basis vectors to store T*V (each TV[i] is the coefficient for \lambda^i of T*V for the extended T) */
57 : BV *AX; /* work basis vectors to store A_i*X for locked eigenvectors */
58 : BV N[2]; /* auxiliary work BVs */
59 : BV X; /* locked eigenvectors */
60 : PetscScalar *T; /* matrix of the invariant pair */
61 : PetscScalar *Tj; /* matrix containing the powers of the invariant pair matrix */
62 : PetscScalar *XpX; /* X^H*X */
63 : PetscInt ld; /* leading dimension for Tj and XpX */
64 : PC pcshell; /* preconditioner including basic precond+projector */
65 : Mat Pshell; /* auxiliary shell matrix */
66 : PetscInt nlock; /* number of locked vectors in the invariant pair */
67 : Vec vtempl; /* reference nested vector */
68 : PetscInt midx; /* minimality index */
69 : PetscInt mmidx; /* maximum allowed minimality index */
70 : PEPJDProjection proj; /* projection type (orthogonal, harmonic) */
71 : } PEP_JD;
72 :
73 : typedef struct {
74 : PEP pep;
75 : PC pc; /* basic preconditioner */
76 : Vec Bp[2]; /* preconditioned residual of derivative polynomial, B\p */
77 : Vec u[2]; /* Ritz vector */
78 : PetscScalar gamma[2]; /* precomputed scalar u'*B\p */
79 : PetscScalar theta;
80 : PetscScalar *M;
81 : PetscScalar *ps;
82 : PetscInt ld;
83 : Vec *work;
84 : Mat PPr;
85 : BV X;
86 : PetscInt n;
87 : } PEP_JD_PCSHELL;
88 :
89 : typedef struct {
90 : Mat Pr,Pi; /* matrix polynomial evaluated at theta */
91 : PEP pep;
92 : Vec *work;
93 : PetscScalar theta[2];
94 : } PEP_JD_MATSHELL;
95 :
96 : /*
97 : Duplicate and resize auxiliary basis
98 : */
99 55 : static PetscErrorCode PEPJDDuplicateBasis(PEP pep,BV *basis)
100 : {
101 55 : PEP_JD *pjd = (PEP_JD*)pep->data;
102 55 : PetscInt nloc,m;
103 55 : BVType type;
104 55 : BVOrthogType otype;
105 55 : BVOrthogRefineType oref;
106 55 : PetscReal oeta;
107 55 : BVOrthogBlockType oblock;
108 55 : VecType vtype;
109 :
110 55 : PetscFunctionBegin;
111 55 : if (pjd->ld>1) {
112 55 : PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),basis));
113 55 : PetscCall(BVGetSizes(pep->V,&nloc,NULL,&m));
114 55 : nloc += pjd->ld-1;
115 55 : PetscCall(BVSetSizes(*basis,nloc,PETSC_DECIDE,m));
116 55 : PetscCall(BVGetType(pep->V,&type));
117 55 : PetscCall(BVSetType(*basis,type));
118 55 : PetscCall(BVGetVecType(pep->V,&vtype));
119 55 : PetscCall(BVSetVecType(*basis,vtype));
120 55 : PetscCall(BVGetOrthogonalization(pep->V,&otype,&oref,&oeta,&oblock));
121 55 : PetscCall(BVSetOrthogonalization(*basis,otype,oref,oeta,oblock));
122 55 : PetscCall(PetscObjectStateIncrease((PetscObject)*basis));
123 0 : } else PetscCall(BVDuplicate(pep->V,basis));
124 55 : PetscFunctionReturn(PETSC_SUCCESS);
125 : }
126 :
127 11 : static PetscErrorCode PEPSetUp_JD(PEP pep)
128 : {
129 11 : PEP_JD *pjd = (PEP_JD*)pep->data;
130 11 : PetscBool isprecond,flg;
131 11 : PetscRandom rand;
132 11 : PetscInt i;
133 :
134 11 : PetscFunctionBegin;
135 11 : PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
136 11 : if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
137 11 : if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;
138 11 : PetscCheck(pep->which==PEP_TARGET_MAGNITUDE || pep->which==PEP_TARGET_REAL || pep->which==PEP_TARGET_IMAGINARY,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver supports only target which, see PEPSetWhichEigenpairs()");
139 :
140 11 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STPRECOND,&isprecond));
141 11 : PetscCheck(isprecond,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver only works with PRECOND spectral transformation");
142 :
143 11 : PetscCall(STGetTransform(pep->st,&flg));
144 11 : PetscCheck(!flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver requires the ST transform flag unset, see STSetTransform()");
145 11 : PEPCheckIgnored(pep,PEP_FEATURE_EXTRACT);
146 :
147 11 : if (!pjd->mmidx) pjd->mmidx = pep->nmat-1;
148 11 : pjd->mmidx = PetscMin(pjd->mmidx,pep->nmat-1);
149 11 : if (!pjd->keep) pjd->keep = 0.5;
150 11 : PetscCall(PEPBasisCoefficients(pep,pep->pbc));
151 11 : PetscCall(PEPAllocateSolution(pep,0));
152 11 : PetscCall(BVGetRandomContext(pep->V,&rand)); /* make sure the random context is available when duplicating */
153 11 : PetscCall(PEPSetWorkVecs(pep,5));
154 11 : pjd->ld = pep->nev;
155 : #if !defined (PETSC_USE_COMPLEX)
156 11 : pjd->ld++;
157 : #endif
158 11 : PetscCall(PetscMalloc2(pep->nmat,&pjd->TV,pep->nmat,&pjd->AX));
159 46 : for (i=0;i<pep->nmat;i++) PetscCall(PEPJDDuplicateBasis(pep,pjd->TV+i));
160 11 : if (pjd->ld>1) {
161 11 : PetscCall(PEPJDDuplicateBasis(pep,&pjd->V));
162 11 : PetscCall(BVSetFromOptions(pjd->V));
163 46 : for (i=0;i<pep->nmat;i++) PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->AX+i));
164 11 : PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->N));
165 11 : PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->N+1));
166 11 : pjd->X = pep->V;
167 11 : PetscCall(PetscCalloc3(pjd->ld*pjd->ld,&pjd->XpX,pep->ncv*pep->ncv,&pjd->T,pjd->ld*pjd->ld*pep->nmat,&pjd->Tj));
168 0 : } else pjd->V = pep->V;
169 11 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(PEPJDDuplicateBasis(pep,&pjd->W));
170 2 : else pjd->W = pjd->V;
171 11 : PetscCall(DSSetType(pep->ds,DSPEP));
172 11 : PetscCall(DSPEPSetDegree(pep->ds,pep->nmat-1));
173 11 : if (pep->basis!=PEP_BASIS_MONOMIAL) PetscCall(DSPEPSetCoefficients(pep->ds,pep->pbc));
174 11 : PetscCall(DSAllocate(pep->ds,pep->ncv));
175 11 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 : /*
179 : Updates columns (low to (high-1)) of TV[i]
180 : */
181 264 : static PetscErrorCode PEPJDUpdateTV(PEP pep,PetscInt low,PetscInt high,Vec *w)
182 : {
183 264 : PEP_JD *pjd = (PEP_JD*)pep->data;
184 264 : PetscInt pp,col,i,nloc,nconv;
185 264 : Vec v1,v2,t1,t2;
186 264 : PetscScalar *array1,*array2,*x2,*xx,*N,*Np,*y2=NULL,zero=0.0,sone=1.0,*pT,fact,*psc;
187 264 : PetscReal *cg,*ca,*cb;
188 264 : PetscMPIInt rk,np;
189 264 : PetscBLASInt n_,ld_,one=1;
190 264 : Mat T;
191 264 : BV pbv;
192 :
193 264 : PetscFunctionBegin;
194 264 : ca = pep->pbc; cb = ca+pep->nmat; cg = cb + pep->nmat;
195 264 : nconv = pjd->nlock;
196 264 : PetscCall(PetscMalloc5(nconv,&x2,nconv,&xx,nconv*nconv,&pT,nconv*nconv,&N,nconv*nconv,&Np));
197 264 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
198 264 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
199 264 : PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
200 264 : t1 = w[0];
201 264 : t2 = w[1];
202 264 : PetscCall(PetscBLASIntCast(pjd->nlock,&n_));
203 264 : PetscCall(PetscBLASIntCast(pjd->ld,&ld_));
204 264 : if (nconv) {
205 189 : for (i=0;i<nconv;i++) PetscCall(PetscArraycpy(pT+i*nconv,pjd->T+i*pep->ncv,nconv));
206 68 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,pT,&T));
207 : }
208 835 : for (col=low;col<high;col++) {
209 571 : PetscCall(BVGetColumn(pjd->V,col,&v1));
210 571 : PetscCall(VecGetArray(v1,&array1));
211 571 : if (nconv>0) {
212 792 : for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
213 : }
214 571 : PetscCall(VecPlaceArray(t1,array1));
215 571 : if (nconv) {
216 258 : PetscCall(BVSetActiveColumns(pjd->N[0],0,nconv));
217 258 : PetscCall(BVSetActiveColumns(pjd->N[1],0,nconv));
218 258 : PetscCall(BVDotVec(pjd->X,t1,xx));
219 : }
220 2325 : for (pp=pep->nmat-1;pp>=0;pp--) {
221 1754 : PetscCall(BVGetColumn(pjd->TV[pp],col,&v2));
222 1754 : PetscCall(VecGetArray(v2,&array2));
223 1754 : PetscCall(VecPlaceArray(t2,array2));
224 1754 : PetscCall(MatMult(pep->A[pp],t1,t2));
225 1754 : if (nconv) {
226 800 : if (pp<pep->nmat-3) {
227 26 : PetscCall(BVMult(pjd->N[0],1.0,-cg[pp+2],pjd->AX[pp+1],NULL));
228 26 : PetscCall(MatShift(T,-cb[pp+1]));
229 26 : PetscCall(BVMult(pjd->N[0],1.0/ca[pp],1.0/ca[pp],pjd->N[1],T));
230 26 : pbv = pjd->N[0]; pjd->N[0] = pjd->N[1]; pjd->N[1] = pbv;
231 26 : PetscCall(BVMultVec(pjd->N[1],1.0,1.0,t2,x2));
232 26 : PetscCall(MatShift(T,cb[pp+1]));
233 774 : } else if (pp==pep->nmat-3) {
234 258 : PetscCall(BVCopy(pjd->AX[pp+2],pjd->N[0]));
235 258 : PetscCall(BVScale(pjd->N[0],1/ca[pp+1]));
236 258 : PetscCall(BVCopy(pjd->AX[pp+1],pjd->N[1]));
237 258 : PetscCall(MatShift(T,-cb[pp+1]));
238 258 : PetscCall(BVMult(pjd->N[1],1.0/ca[pp],1.0/ca[pp],pjd->N[0],T));
239 258 : PetscCall(BVMultVec(pjd->N[1],1.0,1.0,t2,x2));
240 258 : PetscCall(MatShift(T,cb[pp+1]));
241 516 : } else if (pp==pep->nmat-2) PetscCall(BVMultVec(pjd->AX[pp+1],1.0/ca[pp],1.0,t2,x2));
242 800 : if (pp<pjd->midx) {
243 397 : y2 = array2+nloc;
244 397 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n_,&n_,&sone,pjd->Tj+pjd->ld*pjd->ld*pp,&ld_,xx,&one,&zero,y2,&one));
245 397 : if (pp<pjd->midx-2) {
246 0 : fact = -cg[pp+2];
247 0 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&fact,Np,&n_));
248 0 : fact = 1/ca[pp];
249 0 : PetscCall(MatShift(T,-cb[pp+1]));
250 0 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&fact,N,&n_,pT,&n_,&fact,Np,&n_));
251 0 : PetscCall(MatShift(T,cb[pp+1]));
252 0 : psc = Np; Np = N; N = psc;
253 0 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
254 397 : } else if (pp==pjd->midx-2) {
255 139 : fact = 1/ca[pp];
256 139 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&fact,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&zero,N,&n_));
257 139 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
258 258 : } else if (pp==pjd->midx-1) PetscCall(PetscArrayzero(Np,nconv*nconv));
259 : }
260 2428 : for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
261 : }
262 1754 : PetscCall(VecResetArray(t2));
263 1754 : PetscCall(VecRestoreArray(v2,&array2));
264 1754 : PetscCall(BVRestoreColumn(pjd->TV[pp],col,&v2));
265 : }
266 571 : PetscCall(VecResetArray(t1));
267 571 : PetscCall(VecRestoreArray(v1,&array1));
268 571 : PetscCall(BVRestoreColumn(pjd->V,col,&v1));
269 : }
270 264 : if (nconv) PetscCall(MatDestroy(&T));
271 264 : PetscCall(PetscFree5(x2,xx,pT,N,Np));
272 264 : PetscFunctionReturn(PETSC_SUCCESS);
273 : }
274 :
275 : /*
276 : RRQR of X. Xin*P=Xou*R. Rank of R is rk
277 : */
278 47 : static PetscErrorCode PEPJDOrthogonalize(PetscInt row,PetscInt col,PetscScalar *X,PetscInt ldx,PetscInt *rk,PetscInt *P,PetscScalar *R,PetscInt ldr)
279 : {
280 47 : PetscInt i,j,n,r;
281 47 : PetscBLASInt row_,col_,ldx_,*p,lwork,info,n_;
282 47 : PetscScalar *tau,*work;
283 47 : PetscReal tol,*rwork;
284 :
285 47 : PetscFunctionBegin;
286 47 : PetscCall(PetscBLASIntCast(row,&row_));
287 47 : PetscCall(PetscBLASIntCast(col,&col_));
288 47 : PetscCall(PetscBLASIntCast(ldx,&ldx_));
289 47 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
290 47 : n = PetscMin(row,col);
291 47 : PetscCall(PetscBLASIntCast(n,&n_));
292 47 : lwork = 3*col_+1;
293 47 : PetscCall(PetscMalloc4(col,&p,n,&tau,lwork,&work,2*col,&rwork));
294 442 : for (i=1;i<col;i++) p[i] = 0;
295 47 : p[0] = 1;
296 :
297 : /* rank revealing QR */
298 : #if defined(PETSC_USE_COMPLEX)
299 : PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,rwork,&info));
300 : #else
301 47 : PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,&info));
302 : #endif
303 47 : SlepcCheckLapackInfo("geqp3",info);
304 241 : if (P) for (i=0;i<col;i++) P[i] = p[i]-1;
305 :
306 : /* rank computation */
307 47 : tol = PetscMax(row,col)*PETSC_MACHINE_EPSILON*PetscAbsScalar(X[0]);
308 47 : r = 1;
309 442 : for (i=1;i<n;i++) {
310 395 : if (PetscAbsScalar(X[i+ldx*i])>tol) r++;
311 : else break;
312 : }
313 47 : if (rk) *rk=r;
314 :
315 : /* copy upper triangular matrix if requested */
316 47 : if (R) {
317 209 : for (i=0;i<r;i++) {
318 194 : PetscCall(PetscArrayzero(R+i*ldr,r));
319 1703 : for (j=0;j<=i;j++) R[i*ldr+j] = X[i*ldx+j];
320 : }
321 : }
322 47 : PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&row_,&n_,&n_,X,&ldx_,tau,work,&lwork,&info));
323 47 : SlepcCheckLapackInfo("orgqr",info);
324 47 : PetscCall(PetscFPTrapPop());
325 47 : PetscCall(PetscFree4(p,tau,work,rwork));
326 47 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
328 :
329 : /*
330 : Application of extended preconditioner
331 : */
332 4925 : static PetscErrorCode PEPJDExtendedPCApply(PC pc,Vec x,Vec y)
333 : {
334 4925 : PetscInt i,j,nloc,n,ld=0;
335 4925 : PetscMPIInt np;
336 4925 : Vec tx,ty;
337 4925 : PEP_JD_PCSHELL *ctx;
338 4925 : const PetscScalar *array1;
339 4925 : PetscScalar *x2=NULL,*t=NULL,*ps=NULL,*array2,zero=0.0,sone=1.0;
340 4925 : PetscBLASInt one=1,ld_,n_,ncv_;
341 4925 : PEP_JD *pjd=NULL;
342 :
343 4925 : PetscFunctionBegin;
344 4925 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&np));
345 4925 : PetscCall(PCShellGetContext(pc,&ctx));
346 4925 : n = ctx->n;
347 4925 : if (n) {
348 2196 : pjd = (PEP_JD*)ctx->pep->data;
349 2196 : ps = ctx->ps;
350 2196 : ld = pjd->ld;
351 2196 : PetscCall(PetscMalloc2(n,&x2,n,&t));
352 2196 : PetscCall(VecGetLocalSize(ctx->work[0],&nloc));
353 2196 : PetscCall(VecGetArrayRead(x,&array1));
354 4982 : for (i=0;i<n;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
355 2196 : PetscCall(VecRestoreArrayRead(x,&array1));
356 : }
357 :
358 : /* y = B\x apply PC */
359 4925 : tx = ctx->work[0];
360 4925 : ty = ctx->work[1];
361 4925 : PetscCall(VecGetArrayRead(x,&array1));
362 4925 : PetscCall(VecPlaceArray(tx,array1));
363 4925 : PetscCall(VecGetArray(y,&array2));
364 4925 : PetscCall(VecPlaceArray(ty,array2));
365 4925 : PetscCall(PCApply(ctx->pc,tx,ty));
366 4925 : if (n) {
367 2196 : PetscCall(PetscBLASIntCast(ld,&ld_));
368 2196 : PetscCall(PetscBLASIntCast(n,&n_));
369 4982 : for (i=0;i<n;i++) {
370 2786 : t[i] = 0.0;
371 7096 : for (j=0;j<n;j++) t[i] += ctx->M[i+j*ld]*x2[j];
372 : }
373 2196 : if (pjd->midx==1) {
374 1808 : PetscCall(PetscBLASIntCast(ctx->pep->ncv,&ncv_));
375 3664 : for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] -= ctx->theta;
376 1808 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,pjd->T,&ncv_,t,&one,&zero,x2,&one));
377 3664 : for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] += ctx->theta;
378 3664 : for (i=0;i<n;i++) array2[nloc+i] = x2[i];
379 3664 : for (i=0;i<n;i++) x2[i] = -t[i];
380 : } else {
381 1318 : for (i=0;i<n;i++) array2[nloc+i] = t[i];
382 388 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,ps,&ld_,t,&one,&zero,x2,&one));
383 : }
384 4982 : for (i=0;i<n;i++) array2[nloc+i] /= PetscSqrtReal(np);
385 2196 : PetscCall(BVSetActiveColumns(pjd->X,0,n));
386 2196 : PetscCall(BVMultVec(pjd->X,-1.0,1.0,ty,x2));
387 2196 : PetscCall(PetscFree2(x2,t));
388 : }
389 4925 : PetscCall(VecResetArray(tx));
390 4925 : PetscCall(VecResetArray(ty));
391 4925 : PetscCall(VecRestoreArrayRead(x,&array1));
392 4925 : PetscCall(VecRestoreArray(y,&array2));
393 4925 : PetscFunctionReturn(PETSC_SUCCESS);
394 : }
395 :
396 : /*
397 : Application of shell preconditioner:
398 : y = B\x - eta*B\p, with eta = (u'*B\x)/(u'*B\p)
399 : */
400 4437 : static PetscErrorCode PCShellApply_PEPJD(PC pc,Vec x,Vec y)
401 : {
402 4437 : PetscScalar rr,eta;
403 4437 : PEP_JD_PCSHELL *ctx;
404 4437 : PetscInt sz;
405 4437 : const Vec *xs,*ys;
406 : #if !defined(PETSC_USE_COMPLEX)
407 4437 : PetscScalar rx,xr,xx;
408 : #endif
409 :
410 4437 : PetscFunctionBegin;
411 4437 : PetscCall(PCShellGetContext(pc,&ctx));
412 4437 : PetscCall(VecCompGetSubVecs(x,&sz,&xs));
413 4437 : PetscCall(VecCompGetSubVecs(y,NULL,&ys));
414 : /* y = B\x apply extended PC */
415 4437 : PetscCall(PEPJDExtendedPCApply(pc,xs[0],ys[0]));
416 : #if !defined(PETSC_USE_COMPLEX)
417 4437 : if (sz==2) PetscCall(PEPJDExtendedPCApply(pc,xs[1],ys[1]));
418 : #endif
419 :
420 : /* Compute eta = u'*y / u'*Bp */
421 4437 : PetscCall(VecDot(ys[0],ctx->u[0],&rr));
422 4437 : eta = -rr*ctx->gamma[0];
423 : #if !defined(PETSC_USE_COMPLEX)
424 4437 : if (sz==2) {
425 227 : PetscCall(VecDot(ys[0],ctx->u[1],&xr));
426 227 : PetscCall(VecDot(ys[1],ctx->u[0],&rx));
427 227 : PetscCall(VecDot(ys[1],ctx->u[1],&xx));
428 227 : eta += -ctx->gamma[0]*xx-ctx->gamma[1]*(-xr+rx);
429 : }
430 : #endif
431 4437 : eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
432 :
433 : /* y = y - eta*Bp */
434 4437 : PetscCall(VecAXPY(ys[0],eta,ctx->Bp[0]));
435 : #if !defined(PETSC_USE_COMPLEX)
436 4437 : if (sz==2) {
437 227 : PetscCall(VecAXPY(ys[1],eta,ctx->Bp[1]));
438 227 : eta = -ctx->gamma[1]*(rr+xx)+ctx->gamma[0]*(-xr+rx);
439 227 : eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
440 227 : PetscCall(VecAXPY(ys[0],eta,ctx->Bp[1]));
441 227 : PetscCall(VecAXPY(ys[1],-eta,ctx->Bp[0]));
442 : }
443 : #endif
444 4437 : PetscFunctionReturn(PETSC_SUCCESS);
445 : }
446 :
447 219 : static PetscErrorCode PEPJDCopyToExtendedVec(PEP pep,Vec v,PetscScalar *a,PetscInt na,PetscInt off,Vec vex,PetscBool back)
448 : {
449 219 : PetscMPIInt np,rk,count;
450 219 : PetscScalar *array1,*array2;
451 219 : PetscInt nloc;
452 :
453 219 : PetscFunctionBegin;
454 219 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
455 219 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
456 219 : PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
457 219 : if (v) {
458 33 : PetscCall(VecGetArray(v,&array1));
459 33 : PetscCall(VecGetArray(vex,&array2));
460 33 : if (back) PetscCall(PetscArraycpy(array1,array2,nloc));
461 0 : else PetscCall(PetscArraycpy(array2,array1,nloc));
462 33 : PetscCall(VecRestoreArray(v,&array1));
463 33 : PetscCall(VecRestoreArray(vex,&array2));
464 : }
465 219 : if (a) {
466 219 : PetscCall(VecGetArray(vex,&array2));
467 219 : if (back) {
468 33 : PetscCall(PetscArraycpy(a,array2+nloc+off,na));
469 33 : PetscCall(PetscMPIIntCast(na,&count));
470 66 : PetscCallMPI(MPI_Bcast(a,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
471 : } else {
472 186 : PetscCall(PetscArraycpy(array2+nloc+off,a,na));
473 186 : PetscCall(PetscMPIIntCast(na,&count));
474 372 : PetscCallMPI(MPI_Bcast(array2+nloc+off,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
475 : }
476 219 : PetscCall(VecRestoreArray(vex,&array2));
477 : }
478 219 : PetscFunctionReturn(PETSC_SUCCESS);
479 : }
480 :
481 : /* Computes Phi^hat(lambda) times a vector or its derivative (depends on beval)
482 : if no vector is provided returns a matrix
483 : */
484 9293 : static PetscErrorCode PEPJDEvaluateHatBasis(PEP pep,PetscInt n,PetscScalar *H,PetscInt ldh,PetscScalar *beval,PetscScalar *t,PetscInt idx,PetscScalar *qpp,PetscScalar *qp,PetscScalar *q)
485 : {
486 9293 : PetscInt j,i;
487 9293 : PetscBLASInt n_,ldh_,one=1;
488 9293 : PetscReal *a,*b,*g;
489 9293 : PetscScalar sone=1.0,zero=0.0;
490 :
491 9293 : PetscFunctionBegin;
492 9293 : a = pep->pbc; b=a+pep->nmat; g=b+pep->nmat;
493 9293 : PetscCall(PetscBLASIntCast(n,&n_));
494 9293 : PetscCall(PetscBLASIntCast(ldh,&ldh_));
495 9293 : if (idx<1) PetscCall(PetscArrayzero(q,t?n:n*n));
496 6730 : else if (idx==1) {
497 6107 : if (t) {for (j=0;j<n;j++) q[j] = t[j]*beval[idx-1]/a[0];}
498 : else {
499 0 : PetscCall(PetscArrayzero(q,n*n));
500 0 : for (j=0;j<n;j++) q[(j+1)*n] = beval[idx-1]/a[0];
501 : }
502 : } else {
503 4167 : if (t) {
504 4144 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,H,&ldh_,qp,&one,&zero,q,&one));
505 9269 : for (j=0;j<n;j++) {
506 5125 : q[j] += beval[idx-1]*t[j]-b[idx-1]*qp[j]-g[idx-1]*qpp[j];
507 5125 : q[j] /= a[idx-1];
508 : }
509 : } else {
510 23 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,H,&ldh_,qp,&n_,&zero,q,&n_));
511 81 : for (j=0;j<n;j++) {
512 58 : q[j+n*j] += beval[idx-1];
513 210 : for (i=0;i<n;i++) {
514 152 : q[i+n*j] += -b[idx-1]*qp[j*n+i]-g[idx-1]*qpp[j*n+i];
515 152 : q[i+n*j] /= a[idx-1];
516 : }
517 : }
518 : }
519 : }
520 9293 : PetscFunctionReturn(PETSC_SUCCESS);
521 : }
522 :
523 705 : static PetscErrorCode PEPJDComputeResidual(PEP pep,PetscBool derivative,PetscInt sz,Vec *u,PetscScalar *theta,Vec *p,Vec *work)
524 : {
525 705 : PEP_JD *pjd = (PEP_JD*)pep->data;
526 705 : PetscMPIInt rk,np,count;
527 705 : Vec tu,tp,w;
528 705 : PetscScalar *dval,*dvali,*array1,*array2,*x2=NULL,*y2,*qj=NULL,*tt=NULL,*xx=NULL,*xxi=NULL,sone=1.0;
529 705 : PetscInt i,j,nconv,nloc;
530 705 : PetscBLASInt n,ld,one=1;
531 : #if !defined(PETSC_USE_COMPLEX)
532 705 : Vec tui=NULL,tpi=NULL;
533 705 : PetscScalar *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
534 : #endif
535 :
536 705 : PetscFunctionBegin;
537 705 : nconv = pjd->nlock;
538 705 : if (!nconv) PetscCall(PetscMalloc1(2*sz*pep->nmat,&dval));
539 : else {
540 623 : PetscCall(PetscMalloc5(2*pep->nmat,&dval,2*nconv,&xx,nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*pep->nmat,&qj));
541 316 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
542 316 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
543 316 : PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
544 316 : PetscCall(VecGetArray(u[0],&array1));
545 942 : for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]*PetscSqrtReal(np);
546 316 : PetscCall(VecRestoreArray(u[0],&array1));
547 : #if !defined(PETSC_USE_COMPLEX)
548 316 : if (sz==2) {
549 9 : x2i = x2+nconv;
550 9 : PetscCall(VecGetArray(u[1],&arrayi1));
551 27 : for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]*PetscSqrtReal(np);
552 9 : PetscCall(VecRestoreArray(u[1],&arrayi1));
553 : }
554 : #endif
555 : }
556 705 : dvali = dval+pep->nmat;
557 705 : tu = work[0];
558 705 : tp = work[1];
559 705 : w = work[2];
560 705 : PetscCall(VecGetArray(u[0],&array1));
561 705 : PetscCall(VecPlaceArray(tu,array1));
562 705 : PetscCall(VecGetArray(p[0],&array2));
563 705 : PetscCall(VecPlaceArray(tp,array2));
564 705 : PetscCall(VecSet(tp,0.0));
565 : #if !defined(PETSC_USE_COMPLEX)
566 705 : if (sz==2) {
567 57 : tui = work[3];
568 57 : tpi = work[4];
569 57 : PetscCall(VecGetArray(u[1],&arrayi1));
570 57 : PetscCall(VecPlaceArray(tui,arrayi1));
571 57 : PetscCall(VecGetArray(p[1],&arrayi2));
572 57 : PetscCall(VecPlaceArray(tpi,arrayi2));
573 57 : PetscCall(VecSet(tpi,0.0));
574 : }
575 : #endif
576 705 : if (derivative) PetscCall(PEPEvaluateBasisDerivative(pep,theta[0],theta[1],dval,dvali));
577 459 : else PetscCall(PEPEvaluateBasis(pep,theta[0],theta[1],dval,dvali));
578 2647 : for (i=derivative?1:0;i<pep->nmat;i++) {
579 1942 : PetscCall(MatMult(pep->A[i],tu,w));
580 1942 : PetscCall(VecAXPY(tp,dval[i],w));
581 : #if !defined(PETSC_USE_COMPLEX)
582 1942 : if (sz==2) {
583 147 : PetscCall(VecAXPY(tpi,dvali[i],w));
584 147 : PetscCall(MatMult(pep->A[i],tui,w));
585 147 : PetscCall(VecAXPY(tpi,dval[i],w));
586 1942 : PetscCall(VecAXPY(tp,-dvali[i],w));
587 : }
588 : #endif
589 : }
590 705 : if (nconv) {
591 1305 : for (i=0;i<pep->nmat;i++) PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dval,x2,i,i>1?qj+(i-2)*nconv:NULL,i>0?qj+(i-1)*nconv:NULL,qj+i*nconv));
592 : #if !defined(PETSC_USE_COMPLEX)
593 316 : if (sz==2) {
594 9 : qji = qj+nconv*pep->nmat;
595 9 : qq = qji+nconv*pep->nmat;
596 36 : for (i=0;i<pep->nmat;i++) PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dvali,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
597 63 : for (i=0;i<nconv*pep->nmat;i++) qj[i] -= qji[i];
598 36 : for (i=0;i<pep->nmat;i++) {
599 27 : PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dval,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
600 27 : PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dvali,x2,i,i>1?qq+(i-2)*nconv:NULL,i>0?qq+(i-1)*nconv:NULL,qq+i*nconv));
601 : }
602 63 : for (i=0;i<nconv*pep->nmat;i++) qji[i] += qq[i];
603 32 : for (i=derivative?2:1;i<pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tpi,qji+i*nconv));
604 : }
605 : #endif
606 1241 : for (i=derivative?2:1;i<pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tp,qj+i*nconv));
607 :
608 : /* extended vector part */
609 316 : PetscCall(BVSetActiveColumns(pjd->X,0,nconv));
610 316 : PetscCall(BVDotVec(pjd->X,tu,xx));
611 316 : xxi = xx+nconv;
612 : #if !defined(PETSC_USE_COMPLEX)
613 316 : if (sz==2) PetscCall(BVDotVec(pjd->X,tui,xxi));
614 : #endif
615 316 : if (sz==1) PetscCall(PetscArrayzero(xxi,nconv));
616 316 : if (rk==np-1) {
617 277 : PetscCall(PetscBLASIntCast(nconv,&n));
618 277 : PetscCall(PetscBLASIntCast(pjd->ld,&ld));
619 277 : y2 = array2+nloc;
620 277 : PetscCall(PetscArrayzero(y2,nconv));
621 619 : for (j=derivative?1:0;j<pjd->midx;j++) {
622 1081 : for (i=0;i<nconv;i++) tt[i] = dval[j]*xx[i]-dvali[j]*xxi[i];
623 342 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
624 342 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
625 : }
626 805 : for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
627 : #if !defined(PETSC_USE_COMPLEX)
628 277 : if (sz==2) {
629 9 : y2i = arrayi2+nloc;
630 9 : PetscCall(PetscArrayzero(y2i,nconv));
631 23 : for (j=derivative?1:0;j<pjd->midx;j++) {
632 42 : for (i=0;i<nconv;i++) tt[i] = dval[j]*xxi[i]+dvali[j]*xx[i];
633 14 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
634 14 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
635 : }
636 27 : for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
637 : }
638 : #endif
639 : }
640 316 : PetscCall(PetscMPIIntCast(nconv,&count));
641 632 : PetscCallMPI(MPI_Bcast(array2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
642 : #if !defined(PETSC_USE_COMPLEX)
643 325 : if (sz==2) PetscCallMPI(MPI_Bcast(arrayi2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
644 : #endif
645 : }
646 705 : if (nconv) PetscCall(PetscFree5(dval,xx,tt,x2,qj));
647 389 : else PetscCall(PetscFree(dval));
648 705 : PetscCall(VecResetArray(tu));
649 705 : PetscCall(VecRestoreArray(u[0],&array1));
650 705 : PetscCall(VecResetArray(tp));
651 705 : PetscCall(VecRestoreArray(p[0],&array2));
652 : #if !defined(PETSC_USE_COMPLEX)
653 705 : if (sz==2) {
654 57 : PetscCall(VecResetArray(tui));
655 57 : PetscCall(VecRestoreArray(u[1],&arrayi1));
656 57 : PetscCall(VecResetArray(tpi));
657 57 : PetscCall(VecRestoreArray(p[1],&arrayi2));
658 : }
659 : #endif
660 705 : PetscFunctionReturn(PETSC_SUCCESS);
661 : }
662 :
663 11 : static PetscErrorCode PEPJDProcessInitialSpace(PEP pep,Vec *w)
664 : {
665 11 : PEP_JD *pjd = (PEP_JD*)pep->data;
666 11 : PetscScalar *tt,target[2];
667 11 : Vec vg,wg;
668 11 : PetscInt i;
669 11 : PetscReal norm;
670 :
671 11 : PetscFunctionBegin;
672 11 : PetscCall(PetscMalloc1(pjd->ld-1,&tt));
673 11 : PetscCheck(pep->nini==0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for initial vectors not implemented yet");
674 11 : PetscCall(BVSetRandomColumn(pjd->V,0));
675 44 : for (i=0;i<pjd->ld-1;i++) tt[i] = 0.0;
676 11 : PetscCall(BVGetColumn(pjd->V,0,&vg));
677 11 : PetscCall(PEPJDCopyToExtendedVec(pep,NULL,tt,pjd->ld-1,0,vg,PETSC_FALSE));
678 11 : PetscCall(BVRestoreColumn(pjd->V,0,&vg));
679 11 : PetscCall(BVNormColumn(pjd->V,0,NORM_2,&norm));
680 11 : PetscCall(BVScaleColumn(pjd->V,0,1.0/norm));
681 11 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
682 9 : PetscCall(BVGetColumn(pjd->V,0,&vg));
683 9 : PetscCall(BVGetColumn(pjd->W,0,&wg));
684 9 : PetscCall(VecSet(wg,0.0));
685 9 : target[0] = pep->target; target[1] = 0.0;
686 9 : PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,1,&vg,target,&wg,w));
687 9 : PetscCall(BVRestoreColumn(pjd->W,0,&wg));
688 9 : PetscCall(BVRestoreColumn(pjd->V,0,&vg));
689 9 : PetscCall(BVNormColumn(pjd->W,0,NORM_2,&norm));
690 9 : PetscCall(BVScaleColumn(pjd->W,0,1.0/norm));
691 : }
692 11 : PetscCall(PetscFree(tt));
693 11 : PetscFunctionReturn(PETSC_SUCCESS);
694 : }
695 :
696 4200 : static PetscErrorCode MatMult_PEPJD(Mat P,Vec x,Vec y)
697 : {
698 4200 : PEP_JD_MATSHELL *matctx;
699 4200 : PEP_JD *pjd;
700 4200 : PetscInt i,j,nconv,nloc,nmat,ldt,ncv,sz;
701 4200 : Vec tx,ty;
702 4200 : const Vec *xs,*ys;
703 4200 : PetscScalar *array1,*array2,*x2=NULL,*y2,*tt=NULL,*xx=NULL,*xxi,theta[2],sone=1.0,*qj,*val,*vali=NULL;
704 4200 : PetscBLASInt n,ld,one=1;
705 4200 : PetscMPIInt np;
706 : #if !defined(PETSC_USE_COMPLEX)
707 4200 : Vec txi=NULL,tyi=NULL;
708 4200 : PetscScalar *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
709 : #endif
710 :
711 4200 : PetscFunctionBegin;
712 4200 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)P),&np));
713 4200 : PetscCall(MatShellGetContext(P,&matctx));
714 4200 : pjd = (PEP_JD*)matctx->pep->data;
715 4200 : nconv = pjd->nlock;
716 4200 : nmat = matctx->pep->nmat;
717 4200 : ncv = matctx->pep->ncv;
718 4200 : ldt = pjd->ld;
719 4200 : PetscCall(VecCompGetSubVecs(x,&sz,&xs));
720 4200 : PetscCall(VecCompGetSubVecs(y,NULL,&ys));
721 4200 : theta[0] = matctx->theta[0];
722 4200 : theta[1] = (sz==2)?matctx->theta[1]:0.0;
723 4200 : if (nconv>0) {
724 3880 : PetscCall(PetscMalloc5(nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*nmat,&qj,2*nconv,&xx,2*nmat,&val));
725 1980 : PetscCall(BVGetSizes(matctx->pep->V,&nloc,NULL,NULL));
726 1980 : PetscCall(VecGetArray(xs[0],&array1));
727 4364 : for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
728 1980 : PetscCall(VecRestoreArray(xs[0],&array1));
729 : #if !defined(PETSC_USE_COMPLEX)
730 1980 : if (sz==2) {
731 80 : x2i = x2+nconv;
732 80 : PetscCall(VecGetArray(xs[1],&arrayi1));
733 240 : for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]* PetscSqrtReal(np);
734 80 : PetscCall(VecRestoreArray(xs[1],&arrayi1));
735 : }
736 : #endif
737 1980 : vali = val+nmat;
738 : }
739 4200 : tx = matctx->work[0];
740 4200 : ty = matctx->work[1];
741 4200 : PetscCall(VecGetArray(xs[0],&array1));
742 4200 : PetscCall(VecPlaceArray(tx,array1));
743 4200 : PetscCall(VecGetArray(ys[0],&array2));
744 4200 : PetscCall(VecPlaceArray(ty,array2));
745 4200 : PetscCall(MatMult(matctx->Pr,tx,ty));
746 : #if !defined(PETSC_USE_COMPLEX)
747 4200 : if (sz==2) {
748 203 : txi = matctx->work[2];
749 203 : tyi = matctx->work[3];
750 203 : PetscCall(VecGetArray(xs[1],&arrayi1));
751 203 : PetscCall(VecPlaceArray(txi,arrayi1));
752 203 : PetscCall(VecGetArray(ys[1],&arrayi2));
753 203 : PetscCall(VecPlaceArray(tyi,arrayi2));
754 203 : PetscCall(MatMult(matctx->Pr,txi,tyi));
755 203 : if (theta[1]!=0.0) {
756 144 : PetscCall(MatMult(matctx->Pi,txi,matctx->work[4]));
757 144 : PetscCall(VecAXPY(ty,-1.0,matctx->work[4]));
758 144 : PetscCall(MatMult(matctx->Pi,tx,matctx->work[4]));
759 144 : PetscCall(VecAXPY(tyi,1.0,matctx->work[4]));
760 : }
761 : }
762 : #endif
763 4200 : if (nconv>0) {
764 1980 : PetscCall(PEPEvaluateBasis(matctx->pep,theta[0],theta[1],val,vali));
765 9460 : for (i=0;i<nmat;i++) PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,ncv,val,x2,i,i>1?qj+(i-2)*nconv:NULL,i>0?qj+(i-1)*nconv:NULL,qj+i*nconv));
766 : #if !defined(PETSC_USE_COMPLEX)
767 1980 : if (sz==2) {
768 80 : qji = qj+nconv*nmat;
769 80 : qq = qji+nconv*nmat;
770 320 : for (i=0;i<nmat;i++) PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,vali,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
771 560 : for (i=0;i<nconv*nmat;i++) qj[i] -= qji[i];
772 320 : for (i=0;i<nmat;i++) {
773 240 : PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,val,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
774 240 : PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,vali,x2,i,i>1?qq+(i-2)*nconv:NULL,i>0?qq+(i-1)*nconv:NULL,qq+i*nconv));
775 : }
776 560 : for (i=0;i<nconv*nmat;i++) qji[i] += qq[i];
777 240 : for (i=1;i<matctx->pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tyi,qji+i*nconv));
778 : }
779 : #endif
780 7480 : for (i=1;i<nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,ty,qj+i*nconv));
781 :
782 : /* extended vector part */
783 1980 : PetscCall(BVSetActiveColumns(pjd->X,0,nconv));
784 1980 : PetscCall(BVDotVec(pjd->X,tx,xx));
785 1980 : xxi = xx+nconv;
786 : #if !defined(PETSC_USE_COMPLEX)
787 1980 : if (sz==2) PetscCall(BVDotVec(pjd->X,txi,xxi));
788 : #endif
789 1980 : if (sz==1) PetscCall(PetscArrayzero(xxi,nconv));
790 1980 : PetscCall(PetscBLASIntCast(pjd->nlock,&n));
791 1980 : PetscCall(PetscBLASIntCast(ldt,&ld));
792 1980 : y2 = array2+nloc;
793 1980 : PetscCall(PetscArrayzero(y2,nconv));
794 4208 : for (j=0;j<pjd->midx;j++) {
795 5232 : for (i=0;i<nconv;i++) tt[i] = val[j]*xx[i]-vali[j]*xxi[i];
796 2228 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
797 2228 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
798 : }
799 : #if !defined(PETSC_USE_COMPLEX)
800 1980 : if (sz==2) {
801 80 : y2i = arrayi2+nloc;
802 80 : PetscCall(PetscArrayzero(y2i,nconv));
803 240 : for (j=0;j<pjd->midx;j++) {
804 480 : for (i=0;i<nconv;i++) tt[i] = val[j]*xxi[i]+vali[j]*xx[i];
805 160 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
806 160 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
807 : }
808 240 : for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
809 : }
810 : #endif
811 4364 : for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
812 1980 : PetscCall(PetscFree5(tt,x2,qj,xx,val));
813 : }
814 4200 : PetscCall(VecResetArray(tx));
815 4200 : PetscCall(VecRestoreArray(xs[0],&array1));
816 4200 : PetscCall(VecResetArray(ty));
817 4200 : PetscCall(VecRestoreArray(ys[0],&array2));
818 : #if !defined(PETSC_USE_COMPLEX)
819 4200 : if (sz==2) {
820 203 : PetscCall(VecResetArray(txi));
821 203 : PetscCall(VecRestoreArray(xs[1],&arrayi1));
822 203 : PetscCall(VecResetArray(tyi));
823 203 : PetscCall(VecRestoreArray(ys[1],&arrayi2));
824 : }
825 : #endif
826 4200 : PetscFunctionReturn(PETSC_SUCCESS);
827 : }
828 :
829 11 : static PetscErrorCode MatCreateVecs_PEPJD(Mat A,Vec *right,Vec *left)
830 : {
831 11 : PEP_JD_MATSHELL *matctx;
832 11 : PEP_JD *pjd;
833 11 : PetscInt kspsf=1,i;
834 11 : Vec v[2];
835 :
836 11 : PetscFunctionBegin;
837 11 : PetscCall(MatShellGetContext(A,&matctx));
838 11 : pjd = (PEP_JD*)matctx->pep->data;
839 : #if !defined (PETSC_USE_COMPLEX)
840 11 : kspsf = 2;
841 : #endif
842 33 : for (i=0;i<kspsf;i++) PetscCall(BVCreateVec(pjd->V,v+i));
843 11 : if (right) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,right));
844 11 : if (left) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,left));
845 33 : for (i=0;i<kspsf;i++) PetscCall(VecDestroy(&v[i]));
846 11 : PetscFunctionReturn(PETSC_SUCCESS);
847 : }
848 :
849 119 : static PetscErrorCode PEPJDUpdateExtendedPC(PEP pep,PetscScalar theta)
850 : {
851 119 : PEP_JD *pjd = (PEP_JD*)pep->data;
852 119 : PEP_JD_PCSHELL *pcctx;
853 119 : PetscInt i,j,k,n=pjd->nlock,ld=pjd->ld,deg=pep->nmat-1;
854 119 : PetscScalar *M,*ps,*work,*U,*V,*S,*Sp,*Spp,snone=-1.0,sone=1.0,zero=0.0,*val;
855 119 : PetscReal tol,maxeig=0.0,*sg,*rwork;
856 119 : PetscBLASInt n_,info,ld_,*p,lw_,rk=0;
857 :
858 119 : PetscFunctionBegin;
859 119 : if (n) {
860 48 : PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
861 48 : pcctx->theta = theta;
862 48 : pcctx->n = n;
863 48 : M = pcctx->M;
864 48 : PetscCall(PetscBLASIntCast(n,&n_));
865 48 : PetscCall(PetscBLASIntCast(ld,&ld_));
866 48 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
867 48 : if (pjd->midx==1) {
868 25 : PetscCall(PetscArraycpy(M,pjd->XpX,ld*ld));
869 25 : PetscCall(PetscCalloc2(10*n,&work,n,&p));
870 : } else {
871 23 : ps = pcctx->ps;
872 23 : PetscCall(PetscCalloc7(2*n*n,&U,3*n*n,&S,n,&sg,10*n,&work,5*n,&rwork,n,&p,deg+1,&val));
873 23 : V = U+n*n;
874 : /* pseudo-inverse */
875 81 : for (j=0;j<n;j++) {
876 210 : for (i=0;i<n;i++) S[n*j+i] = -pjd->T[pep->ncv*j+i];
877 58 : S[n*j+j] += theta;
878 : }
879 23 : lw_ = 10*n_;
880 : #if !defined (PETSC_USE_COMPLEX)
881 23 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,&info));
882 : #else
883 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,rwork,&info));
884 : #endif
885 23 : SlepcCheckLapackInfo("gesvd",info);
886 81 : for (i=0;i<n;i++) maxeig = PetscMax(maxeig,sg[i]);
887 23 : tol = 10*PETSC_MACHINE_EPSILON*n*maxeig;
888 81 : for (j=0;j<n;j++) {
889 58 : if (sg[j]>tol) {
890 210 : for (i=0;i<n;i++) U[j*n+i] /= sg[j];
891 58 : rk++;
892 : } else break;
893 : }
894 23 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&rk,&sone,U,&n_,V,&n_,&zero,ps,&ld_));
895 :
896 : /* compute M */
897 23 : PetscCall(PEPEvaluateBasis(pep,theta,0.0,val,NULL));
898 23 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&snone,pjd->XpX,&ld_,ps,&ld_,&zero,M,&ld_));
899 23 : PetscCall(PetscArrayzero(S,2*n*n));
900 23 : Sp = S+n*n;
901 81 : for (j=0;j<n;j++) S[j*(n+1)] = 1.0;
902 46 : for (k=1;k<pjd->midx;k++) {
903 233 : for (j=0;j<n;j++) for (i=0;i<n;i++) V[j*n+i] = S[j*n+i] - ps[j*ld+i]*val[k];
904 23 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pjd->XpX,&ld_,V,&n_,&zero,U,&n_));
905 23 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+k*ld*ld,&ld_,U,&n_,&sone,M,&ld_));
906 23 : Spp = Sp; Sp = S;
907 23 : PetscCall(PEPJDEvaluateHatBasis(pep,n,pjd->T,pep->ncv,val,NULL,k+1,Spp,Sp,S));
908 : }
909 : }
910 : /* inverse */
911 48 : PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,M,&ld_,p,&info));
912 48 : SlepcCheckLapackInfo("getrf",info);
913 48 : PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&n_,M,&ld_,p,work,&n_,&info));
914 48 : SlepcCheckLapackInfo("getri",info);
915 48 : PetscCall(PetscFPTrapPop());
916 48 : if (pjd->midx==1) PetscCall(PetscFree2(work,p));
917 23 : else PetscCall(PetscFree7(U,S,sg,work,rwork,p,val));
918 : }
919 119 : PetscFunctionReturn(PETSC_SUCCESS);
920 : }
921 :
922 248 : static PetscErrorCode PEPJDMatSetUp(PEP pep,PetscInt sz,PetscScalar *theta)
923 : {
924 248 : PEP_JD *pjd = (PEP_JD*)pep->data;
925 248 : PEP_JD_MATSHELL *matctx;
926 248 : PEP_JD_PCSHELL *pcctx;
927 248 : MatStructure str;
928 248 : PetscScalar *vals,*valsi;
929 248 : PetscBool skipmat=PETSC_FALSE;
930 248 : PetscInt i;
931 248 : Mat Pr=NULL;
932 :
933 248 : PetscFunctionBegin;
934 248 : if (sz==2 && theta[1]==0.0) sz = 1;
935 248 : PetscCall(MatShellGetContext(pjd->Pshell,&matctx));
936 248 : PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
937 248 : if (matctx->Pr && matctx->theta[0]==theta[0] && ((!matctx->Pi && sz==1) || (sz==2 && matctx->theta[1]==theta[1]))) {
938 127 : if (pcctx->n == pjd->nlock) PetscFunctionReturn(PETSC_SUCCESS);
939 : skipmat = PETSC_TRUE;
940 : }
941 : if (!skipmat) {
942 121 : PetscCall(PetscMalloc2(pep->nmat,&vals,pep->nmat,&valsi));
943 121 : PetscCall(STGetMatStructure(pep->st,&str));
944 121 : PetscCall(PEPEvaluateBasis(pep,theta[0],theta[1],vals,valsi));
945 121 : if (!matctx->Pr) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pr));
946 110 : else PetscCall(MatCopy(pep->A[0],matctx->Pr,str));
947 374 : for (i=1;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pr,vals[i],pep->A[i],str));
948 121 : if (!pjd->reusepc) {
949 118 : if (pcctx->PPr && sz==2) {
950 6 : PetscCall(MatCopy(matctx->Pr,pcctx->PPr,str));
951 6 : Pr = pcctx->PPr;
952 112 : } else Pr = matctx->Pr;
953 : }
954 121 : matctx->theta[0] = theta[0];
955 : #if !defined(PETSC_USE_COMPLEX)
956 121 : if (sz==2) {
957 8 : if (!matctx->Pi) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pi));
958 5 : else PetscCall(MatCopy(pep->A[1],matctx->Pi,str));
959 8 : PetscCall(MatScale(matctx->Pi,valsi[1]));
960 16 : for (i=2;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pi,valsi[i],pep->A[i],str));
961 8 : matctx->theta[1] = theta[1];
962 : }
963 : #endif
964 121 : PetscCall(PetscFree2(vals,valsi));
965 : }
966 122 : if (!pjd->reusepc) {
967 119 : if (!skipmat) {
968 118 : PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
969 118 : PetscCall(PCSetUp(pcctx->pc));
970 : }
971 119 : PetscCall(PEPJDUpdateExtendedPC(pep,theta[0]));
972 : }
973 122 : PetscFunctionReturn(PETSC_SUCCESS);
974 : }
975 :
976 11 : static PetscErrorCode PEPJDCreateShellPC(PEP pep,Vec *ww)
977 : {
978 11 : PEP_JD *pjd = (PEP_JD*)pep->data;
979 11 : PEP_JD_PCSHELL *pcctx;
980 11 : PEP_JD_MATSHELL *matctx;
981 11 : KSP ksp;
982 11 : PetscInt nloc,mloc,kspsf=1;
983 11 : Vec v[2];
984 11 : PetscScalar target[2];
985 11 : Mat Pr;
986 :
987 11 : PetscFunctionBegin;
988 : /* Create the reference vector */
989 11 : PetscCall(BVGetColumn(pjd->V,0,&v[0]));
990 11 : v[1] = v[0];
991 : #if !defined (PETSC_USE_COMPLEX)
992 11 : kspsf = 2;
993 : #endif
994 11 : PetscCall(VecCreateCompWithVecs(v,kspsf,NULL,&pjd->vtempl));
995 11 : PetscCall(BVRestoreColumn(pjd->V,0,&v[0]));
996 :
997 : /* Replace preconditioner with one containing projectors */
998 11 : PetscCall(PCCreate(PetscObjectComm((PetscObject)pep),&pjd->pcshell));
999 11 : PetscCall(PCSetType(pjd->pcshell,PCSHELL));
1000 11 : PetscCall(PCShellSetName(pjd->pcshell,"PCPEPJD"));
1001 11 : PetscCall(PCShellSetApply(pjd->pcshell,PCShellApply_PEPJD));
1002 11 : PetscCall(PetscNew(&pcctx));
1003 11 : PetscCall(PCShellSetContext(pjd->pcshell,pcctx));
1004 11 : PetscCall(STGetKSP(pep->st,&ksp));
1005 11 : PetscCall(BVCreateVec(pjd->V,&pcctx->Bp[0]));
1006 11 : PetscCall(VecDuplicate(pcctx->Bp[0],&pcctx->Bp[1]));
1007 11 : PetscCall(KSPGetPC(ksp,&pcctx->pc));
1008 11 : PetscCall(PetscObjectReference((PetscObject)pcctx->pc));
1009 11 : PetscCall(MatGetLocalSize(pep->A[0],&mloc,&nloc));
1010 11 : if (pjd->ld>1) {
1011 11 : nloc += pjd->ld-1; mloc += pjd->ld-1;
1012 : }
1013 11 : PetscCall(PetscNew(&matctx));
1014 11 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),kspsf*nloc,kspsf*mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&pjd->Pshell));
1015 11 : PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_MULT,(void(*)(void))MatMult_PEPJD));
1016 11 : PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_PEPJD));
1017 11 : matctx->pep = pep;
1018 11 : target[0] = pep->target; target[1] = 0.0;
1019 11 : PetscCall(PEPJDMatSetUp(pep,1,target));
1020 11 : Pr = matctx->Pr;
1021 11 : pcctx->PPr = NULL;
1022 : #if !defined(PETSC_USE_COMPLEX)
1023 11 : if (!pjd->reusepc) {
1024 10 : PetscCall(MatDuplicate(matctx->Pr,MAT_COPY_VALUES,&pcctx->PPr));
1025 10 : Pr = pcctx->PPr;
1026 : }
1027 : #endif
1028 11 : PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
1029 11 : PetscCall(PCSetErrorIfFailure(pcctx->pc,PETSC_TRUE));
1030 11 : PetscCall(KSPSetPC(ksp,pjd->pcshell));
1031 11 : if (pjd->reusepc) {
1032 1 : PetscCall(PCSetReusePreconditioner(pcctx->pc,PETSC_TRUE));
1033 1 : PetscCall(KSPSetReusePreconditioner(ksp,PETSC_TRUE));
1034 : }
1035 11 : PetscCall(PEP_KSPSetOperators(ksp,pjd->Pshell,pjd->Pshell));
1036 11 : PetscCall(KSPSetUp(ksp));
1037 11 : if (pjd->ld>1) {
1038 11 : PetscCall(PetscMalloc2(pjd->ld*pjd->ld,&pcctx->M,pjd->ld*pjd->ld,&pcctx->ps));
1039 11 : pcctx->pep = pep;
1040 : }
1041 11 : matctx->work = ww;
1042 11 : pcctx->work = ww;
1043 11 : PetscFunctionReturn(PETSC_SUCCESS);
1044 : }
1045 :
1046 11 : static PetscErrorCode PEPJDEigenvectors(PEP pep)
1047 : {
1048 11 : PEP_JD *pjd = (PEP_JD*)pep->data;
1049 11 : PetscBLASInt ld,nconv,info,nc;
1050 11 : PetscScalar *Z;
1051 11 : PetscReal *wr;
1052 11 : Mat U;
1053 : #if defined(PETSC_USE_COMPLEX)
1054 : PetscScalar *w;
1055 : #endif
1056 :
1057 11 : PetscFunctionBegin;
1058 11 : PetscCall(PetscBLASIntCast(pep->ncv,&ld));
1059 11 : PetscCall(PetscBLASIntCast(pep->nconv,&nconv));
1060 11 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1061 : #if !defined(PETSC_USE_COMPLEX)
1062 11 : PetscCall(PetscMalloc2(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr));
1063 11 : PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,wr,&info));
1064 : #else
1065 : PetscCall(PetscMalloc3(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr,2*pep->ncv,&w));
1066 : PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,w,wr,&info));
1067 : #endif
1068 11 : PetscCall(PetscFPTrapPop());
1069 11 : SlepcCheckLapackInfo("trevc",info);
1070 11 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,Z,&U));
1071 11 : PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv));
1072 11 : PetscCall(BVMultInPlace(pjd->X,U,0,pep->nconv));
1073 11 : PetscCall(BVNormalize(pjd->X,pep->eigi));
1074 11 : PetscCall(MatDestroy(&U));
1075 : #if !defined(PETSC_USE_COMPLEX)
1076 11 : PetscCall(PetscFree2(Z,wr));
1077 : #else
1078 : PetscCall(PetscFree3(Z,wr,w));
1079 : #endif
1080 11 : PetscFunctionReturn(PETSC_SUCCESS);
1081 : }
1082 :
1083 15 : static PetscErrorCode PEPJDLockConverged(PEP pep,PetscInt *nv,PetscInt sz)
1084 : {
1085 15 : PEP_JD *pjd = (PEP_JD*)pep->data;
1086 15 : PetscInt j,i,*P,ldds,rk=0,nvv=*nv;
1087 15 : Vec v,x,w;
1088 15 : PetscScalar *R,*r,*pX,target[2];
1089 15 : Mat X;
1090 15 : PetscBLASInt sz_,rk_,nv_,info;
1091 15 : PetscMPIInt np;
1092 :
1093 15 : PetscFunctionBegin;
1094 : /* update AX and XpX */
1095 34 : for (i=sz;i>0;i--) {
1096 19 : PetscCall(BVGetColumn(pjd->X,pjd->nlock-i,&x));
1097 78 : for (j=0;j<pep->nmat;j++) {
1098 59 : PetscCall(BVGetColumn(pjd->AX[j],pjd->nlock-i,&v));
1099 59 : PetscCall(MatMult(pep->A[j],x,v));
1100 59 : PetscCall(BVRestoreColumn(pjd->AX[j],pjd->nlock-i,&v));
1101 59 : PetscCall(BVSetActiveColumns(pjd->AX[j],0,pjd->nlock-i+1));
1102 : }
1103 19 : PetscCall(BVRestoreColumn(pjd->X,pjd->nlock-i,&x));
1104 19 : PetscCall(BVDotColumn(pjd->X,(pjd->nlock-i),pjd->XpX+(pjd->nlock-i)*pjd->ld));
1105 19 : pjd->XpX[(pjd->nlock-i)*(1+pjd->ld)] = 1.0;
1106 35 : for (j=0;j<pjd->nlock-i;j++) pjd->XpX[j*pjd->ld+pjd->nlock-i] = PetscConj(pjd->XpX[(pjd->nlock-i)*pjd->ld+j]);
1107 : }
1108 :
1109 : /* minimality index */
1110 15 : pjd->midx = PetscMin(pjd->mmidx,pjd->nlock);
1111 :
1112 : /* evaluate the polynomial basis in T */
1113 15 : PetscCall(PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat));
1114 62 : for (j=0;j<pep->nmat;j++) PetscCall(PEPEvaluateBasisMat(pep,pjd->nlock,pjd->T,pep->ncv,j,(j>1)?pjd->Tj+(j-2)*pjd->ld*pjd->ld:NULL,pjd->ld,j?pjd->Tj+(j-1)*pjd->ld*pjd->ld:NULL,pjd->ld,pjd->Tj+j*pjd->ld*pjd->ld,pjd->ld));
1115 :
1116 : /* Extend search space */
1117 15 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
1118 15 : PetscCall(PetscCalloc3(nvv,&P,nvv*nvv,&R,nvv*sz,&r));
1119 15 : PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
1120 15 : PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
1121 15 : PetscCall(PEPJDOrthogonalize(nvv,nvv,pX,ldds,&rk,P,R,nvv));
1122 34 : for (j=0;j<sz;j++) {
1123 262 : for (i=0;i<rk;i++) r[i*sz+j] = PetscConj(R[nvv*i+j]*pep->eigr[P[i]]); /* first row scaled with permuted diagonal */
1124 : }
1125 15 : PetscCall(PetscBLASIntCast(rk,&rk_));
1126 15 : PetscCall(PetscBLASIntCast(sz,&sz_));
1127 15 : PetscCall(PetscBLASIntCast(nvv,&nv_));
1128 15 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1129 15 : PetscCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&rk_,R,&nv_,&info));
1130 15 : PetscCall(PetscFPTrapPop());
1131 15 : SlepcCheckLapackInfo("trtri",info);
1132 34 : for (i=0;i<sz;i++) PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&rk_,R,&nv_,r+i,&sz_));
1133 258 : for (i=0;i<sz*rk;i++) r[i] = PetscConj(r[i])/PetscSqrtReal(np); /* revert */
1134 15 : PetscCall(BVSetActiveColumns(pjd->V,0,nvv));
1135 15 : rk -= sz;
1136 190 : for (j=0;j<rk;j++) PetscCall(PetscArraycpy(R+j*nvv,pX+(j+sz)*ldds,nvv));
1137 15 : PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
1138 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nvv,rk,R,&X));
1139 15 : PetscCall(BVMultInPlace(pjd->V,X,0,rk));
1140 15 : PetscCall(MatDestroy(&X));
1141 15 : PetscCall(BVSetActiveColumns(pjd->V,0,rk));
1142 190 : for (j=0;j<rk;j++) {
1143 175 : PetscCall(BVGetColumn(pjd->V,j,&v));
1144 175 : PetscCall(PEPJDCopyToExtendedVec(pep,NULL,r+sz*(j+sz),sz,pjd->nlock-sz,v,PETSC_FALSE));
1145 175 : PetscCall(BVRestoreColumn(pjd->V,j,&v));
1146 : }
1147 15 : PetscCall(BVOrthogonalize(pjd->V,NULL));
1148 :
1149 15 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1150 181 : for (j=0;j<rk;j++) {
1151 : /* W = P(target)*V */
1152 168 : PetscCall(BVGetColumn(pjd->W,j,&w));
1153 168 : PetscCall(BVGetColumn(pjd->V,j,&v));
1154 168 : target[0] = pep->target; target[1] = 0.0;
1155 168 : PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,1,&v,target,&w,pep->work));
1156 168 : PetscCall(BVRestoreColumn(pjd->V,j,&v));
1157 168 : PetscCall(BVRestoreColumn(pjd->W,j,&w));
1158 : }
1159 13 : PetscCall(BVSetActiveColumns(pjd->W,0,rk));
1160 13 : PetscCall(BVOrthogonalize(pjd->W,NULL));
1161 : }
1162 15 : *nv = rk;
1163 15 : PetscCall(PetscFree3(P,R,r));
1164 15 : PetscFunctionReturn(PETSC_SUCCESS);
1165 : }
1166 :
1167 237 : static PetscErrorCode PEPJDSystemSetUp(PEP pep,PetscInt sz,PetscScalar *theta,Vec *u,Vec *p,Vec *ww)
1168 : {
1169 237 : PEP_JD *pjd = (PEP_JD*)pep->data;
1170 237 : PEP_JD_PCSHELL *pcctx;
1171 : #if !defined(PETSC_USE_COMPLEX)
1172 237 : PetscScalar s[2];
1173 : #endif
1174 :
1175 237 : PetscFunctionBegin;
1176 237 : PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
1177 237 : PetscCall(PEPJDMatSetUp(pep,sz,theta));
1178 237 : pcctx->u[0] = u[0]; pcctx->u[1] = u[1];
1179 : /* Compute r'. p is a work space vector */
1180 237 : PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,sz,u,theta,p,ww));
1181 237 : PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[0],pcctx->Bp[0]));
1182 237 : PetscCall(VecDot(pcctx->Bp[0],u[0],pcctx->gamma));
1183 : #if !defined(PETSC_USE_COMPLEX)
1184 237 : if (sz==2) {
1185 24 : PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[1],pcctx->Bp[1]));
1186 24 : PetscCall(VecDot(pcctx->Bp[0],u[1],pcctx->gamma+1));
1187 24 : PetscCall(VecMDot(pcctx->Bp[1],2,u,s));
1188 24 : pcctx->gamma[0] += s[1];
1189 24 : pcctx->gamma[1] = -pcctx->gamma[1]+s[0];
1190 : }
1191 : #endif
1192 237 : if (sz==1) {
1193 213 : PetscCall(VecZeroEntries(pcctx->Bp[1]));
1194 213 : pcctx->gamma[1] = 0.0;
1195 : }
1196 237 : PetscFunctionReturn(PETSC_SUCCESS);
1197 : }
1198 :
1199 11 : static PetscErrorCode PEPSolve_JD(PEP pep)
1200 : {
1201 11 : PEP_JD *pjd = (PEP_JD*)pep->data;
1202 11 : PetscInt k,nv,nvc,ld,minv,dim,bupdated=0,sz=1,kspsf=1,idx,off,maxits,nloc;
1203 11 : PetscMPIInt np,count;
1204 11 : PetscScalar theta[2]={0.0,0.0},ritz[2]={0.0,0.0},*pX,*eig,*eigi,*array;
1205 11 : PetscReal norm,*res,tol=0.0,rtol,abstol, dtol;
1206 11 : PetscBool lindep,ini=PETSC_TRUE;
1207 11 : Vec tc,t[2]={NULL,NULL},u[2]={NULL,NULL},p[2]={NULL,NULL};
1208 11 : Vec rc,rr[2],r[2]={NULL,NULL},*ww=pep->work,v[2];
1209 11 : Mat G,X,Y;
1210 11 : KSP ksp;
1211 11 : PEP_JD_PCSHELL *pcctx;
1212 11 : PEP_JD_MATSHELL *matctx;
1213 : #if !defined(PETSC_USE_COMPLEX)
1214 11 : PetscReal norm1;
1215 : #endif
1216 :
1217 11 : PetscFunctionBegin;
1218 11 : PetscCall(PetscCitationsRegister(citation,&cited));
1219 11 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
1220 11 : PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
1221 11 : PetscCall(DSGetLeadingDimension(pep->ds,&ld));
1222 11 : PetscCall(PetscCalloc3(pep->ncv+pep->nev,&eig,pep->ncv+pep->nev,&eigi,pep->ncv+pep->nev,&res));
1223 11 : pjd->nlock = 0;
1224 11 : PetscCall(STGetKSP(pep->st,&ksp));
1225 11 : PetscCall(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
1226 : #if !defined (PETSC_USE_COMPLEX)
1227 11 : kspsf = 2;
1228 : #endif
1229 11 : PetscCall(PEPJDProcessInitialSpace(pep,ww));
1230 11 : nv = (pep->nini)?pep->nini:1;
1231 :
1232 : /* Replace preconditioner with one containing projectors */
1233 11 : PetscCall(PEPJDCreateShellPC(pep,ww));
1234 11 : PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
1235 :
1236 : /* Create auxiliary vectors */
1237 11 : PetscCall(BVCreateVec(pjd->V,&u[0]));
1238 11 : PetscCall(VecDuplicate(u[0],&p[0]));
1239 11 : PetscCall(VecDuplicate(u[0],&r[0]));
1240 : #if !defined (PETSC_USE_COMPLEX)
1241 11 : PetscCall(VecDuplicate(u[0],&u[1]));
1242 11 : PetscCall(VecDuplicate(u[0],&p[1]));
1243 11 : PetscCall(VecDuplicate(u[0],&r[1]));
1244 : #endif
1245 :
1246 : /* Restart loop */
1247 275 : while (pep->reason == PEP_CONVERGED_ITERATING) {
1248 264 : pep->its++;
1249 264 : PetscCall(DSSetDimensions(pep->ds,nv,0,0));
1250 264 : PetscCall(BVSetActiveColumns(pjd->V,bupdated,nv));
1251 264 : PetscCall(PEPJDUpdateTV(pep,bupdated,nv,ww));
1252 264 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVSetActiveColumns(pjd->W,bupdated,nv));
1253 1084 : for (k=0;k<pep->nmat;k++) {
1254 820 : PetscCall(BVSetActiveColumns(pjd->TV[k],bupdated,nv));
1255 820 : PetscCall(DSGetMat(pep->ds,DSMatExtra[k],&G));
1256 820 : PetscCall(BVMatProject(pjd->TV[k],NULL,pjd->W,G));
1257 820 : PetscCall(DSRestoreMat(pep->ds,DSMatExtra[k],&G));
1258 : }
1259 264 : PetscCall(BVSetActiveColumns(pjd->V,0,nv));
1260 264 : PetscCall(BVSetActiveColumns(pjd->W,0,nv));
1261 :
1262 : /* Solve projected problem */
1263 264 : PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
1264 264 : PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
1265 264 : PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
1266 264 : PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
1267 : idx = 0;
1268 291 : do {
1269 291 : ritz[0] = pep->eigr[idx];
1270 : #if !defined(PETSC_USE_COMPLEX)
1271 291 : ritz[1] = pep->eigi[idx];
1272 291 : sz = (ritz[1]==0.0)?1:2;
1273 : #endif
1274 : /* Compute Ritz vector u=V*X(:,1) */
1275 291 : PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
1276 291 : PetscCall(BVSetActiveColumns(pjd->V,0,nv));
1277 291 : PetscCall(BVMultVec(pjd->V,1.0,0.0,u[0],pX+idx*ld));
1278 : #if !defined(PETSC_USE_COMPLEX)
1279 291 : if (sz==2) PetscCall(BVMultVec(pjd->V,1.0,0.0,u[1],pX+(idx+1)*ld));
1280 : #endif
1281 291 : PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
1282 291 : PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,sz,u,ritz,r,ww));
1283 : /* Check convergence */
1284 291 : PetscCall(VecNorm(r[0],NORM_2,&norm));
1285 : #if !defined(PETSC_USE_COMPLEX)
1286 291 : if (sz==2) {
1287 33 : PetscCall(VecNorm(r[1],NORM_2,&norm1));
1288 33 : norm = SlepcAbs(norm,norm1);
1289 : }
1290 : #endif
1291 291 : PetscCall((*pep->converged)(pep,ritz[0],ritz[1],norm,&pep->errest[pep->nconv],pep->convergedctx));
1292 291 : if (sz==2) pep->errest[pep->nconv+1] = pep->errest[pep->nconv];
1293 291 : if (ini) {
1294 40 : tol = PetscMin(.1,pep->errest[pep->nconv]); ini = PETSC_FALSE;
1295 251 : } else tol = PetscMin(pep->errest[pep->nconv],tol/2);
1296 291 : PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,(pep->errest[pep->nconv]<pep->tol)?pep->nconv+sz:pep->nconv,pep->nev,&pep->reason,pep->stoppingctx));
1297 291 : if (pep->errest[pep->nconv]<pep->tol) {
1298 : /* Ritz pair converged */
1299 29 : ini = PETSC_TRUE;
1300 29 : minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1301 29 : if (pjd->ld>1) {
1302 29 : PetscCall(BVGetColumn(pjd->X,pep->nconv,&v[0]));
1303 29 : PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*pep->nconv,pjd->ld-1,0,u[0],PETSC_TRUE));
1304 29 : PetscCall(BVRestoreColumn(pjd->X,pep->nconv,&v[0]));
1305 29 : PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+1));
1306 29 : PetscCall(BVNormColumn(pjd->X,pep->nconv,NORM_2,&norm));
1307 29 : PetscCall(BVScaleColumn(pjd->X,pep->nconv,1.0/norm));
1308 63 : for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*pep->nconv+k] *= PetscSqrtReal(np)/norm;
1309 29 : pjd->T[(pep->ncv+1)*pep->nconv] = ritz[0];
1310 29 : eig[pep->nconv] = ritz[0];
1311 29 : idx++;
1312 : #if !defined(PETSC_USE_COMPLEX)
1313 29 : if (sz==2) {
1314 4 : PetscCall(BVGetColumn(pjd->X,pep->nconv+1,&v[0]));
1315 4 : PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*(pep->nconv+1),pjd->ld-1,0,u[1],PETSC_TRUE));
1316 4 : PetscCall(BVRestoreColumn(pjd->X,pep->nconv+1,&v[0]));
1317 4 : PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+2));
1318 4 : PetscCall(BVNormColumn(pjd->X,pep->nconv+1,NORM_2,&norm1));
1319 4 : PetscCall(BVScaleColumn(pjd->X,pep->nconv+1,1.0/norm1));
1320 6 : for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*(pep->nconv+1)+k] *= PetscSqrtReal(np)/norm1;
1321 4 : pjd->T[(pep->ncv+1)*(pep->nconv+1)] = ritz[0];
1322 4 : pjd->T[(pep->ncv+1)*pep->nconv+1] = -ritz[1]*norm1/norm;
1323 4 : pjd->T[(pep->ncv+1)*(pep->nconv+1)-1] = ritz[1]*norm/norm1;
1324 4 : eig[pep->nconv+1] = ritz[0];
1325 4 : eigi[pep->nconv] = ritz[1]; eigi[pep->nconv+1] = -ritz[1];
1326 4 : idx++;
1327 : }
1328 : #endif
1329 0 : } else PetscCall(BVInsertVec(pep->V,pep->nconv,u[0]));
1330 29 : pep->nconv += sz;
1331 : }
1332 291 : } while (pep->errest[pep->nconv]<pep->tol && pep->nconv<nv);
1333 :
1334 264 : if (pep->reason==PEP_CONVERGED_ITERATING) {
1335 253 : nvc = nv;
1336 253 : if (idx) {
1337 15 : pjd->nlock +=idx;
1338 15 : PetscCall(PEPJDLockConverged(pep,&nv,idx));
1339 : }
1340 253 : if (nv+sz>=pep->ncv-1) {
1341 : /* Basis full, force restart */
1342 16 : minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1343 16 : PetscCall(DSGetDimensions(pep->ds,&dim,NULL,NULL,NULL));
1344 16 : PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
1345 16 : PetscCall(PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld));
1346 16 : PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
1347 16 : PetscCall(DSGetArray(pep->ds,DS_MAT_Y,&pX));
1348 16 : PetscCall(PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld));
1349 16 : PetscCall(DSRestoreArray(pep->ds,DS_MAT_Y,&pX));
1350 16 : PetscCall(DSGetMat(pep->ds,DS_MAT_X,&X));
1351 16 : PetscCall(BVMultInPlace(pjd->V,X,0,minv));
1352 16 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_X,&X));
1353 16 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1354 12 : PetscCall(DSGetMat(pep->ds,DS_MAT_Y,&Y));
1355 12 : PetscCall(BVMultInPlace(pjd->W,Y,pep->nconv,minv));
1356 12 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_Y,&Y));
1357 : }
1358 16 : nv = minv;
1359 16 : bupdated = 0;
1360 : } else {
1361 237 : if (!idx && pep->errest[pep->nconv]<pjd->fix) {theta[0] = ritz[0]; theta[1] = ritz[1];}
1362 153 : else {theta[0] = pep->target; theta[1] = 0.0;}
1363 : /* Update system mat */
1364 237 : PetscCall(PEPJDSystemSetUp(pep,sz,theta,u,p,ww));
1365 : /* Solve correction equation to expand basis */
1366 237 : PetscCall(BVGetColumn(pjd->V,nv,&t[0]));
1367 237 : rr[0] = r[0];
1368 237 : if (sz==2) {
1369 24 : PetscCall(BVGetColumn(pjd->V,nv+1,&t[1]));
1370 24 : rr[1] = r[1];
1371 : } else {
1372 213 : t[1] = NULL;
1373 213 : rr[1] = NULL;
1374 : }
1375 237 : PetscCall(VecCreateCompWithVecs(t,kspsf,pjd->vtempl,&tc));
1376 237 : PetscCall(VecCreateCompWithVecs(rr,kspsf,pjd->vtempl,&rc));
1377 237 : PetscCall(VecCompSetSubVecs(pjd->vtempl,sz,NULL));
1378 237 : tol = PetscMax(rtol,tol/2);
1379 237 : PetscCall(KSPSetTolerances(ksp,tol,abstol,dtol,maxits));
1380 237 : PetscCall(KSPSolve(ksp,rc,tc));
1381 237 : PetscCall(VecDestroy(&tc));
1382 237 : PetscCall(VecDestroy(&rc));
1383 237 : PetscCall(VecGetArray(t[0],&array));
1384 237 : PetscCall(PetscMPIIntCast(pep->nconv,&count));
1385 474 : PetscCallMPI(MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
1386 237 : PetscCall(VecRestoreArray(t[0],&array));
1387 237 : PetscCall(BVRestoreColumn(pjd->V,nv,&t[0]));
1388 237 : PetscCall(BVOrthogonalizeColumn(pjd->V,nv,NULL,&norm,&lindep));
1389 237 : if (lindep || norm==0.0) {
1390 0 : PetscCheck(sz!=1,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
1391 : off = 1;
1392 : } else {
1393 237 : off = 0;
1394 237 : PetscCall(BVScaleColumn(pjd->V,nv,1.0/norm));
1395 : }
1396 : #if !defined(PETSC_USE_COMPLEX)
1397 237 : if (sz==2) {
1398 24 : PetscCall(VecGetArray(t[1],&array));
1399 48 : PetscCallMPI(MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
1400 24 : PetscCall(VecRestoreArray(t[1],&array));
1401 24 : PetscCall(BVRestoreColumn(pjd->V,nv+1,&t[1]));
1402 24 : if (off) PetscCall(BVCopyColumn(pjd->V,nv+1,nv));
1403 24 : PetscCall(BVOrthogonalizeColumn(pjd->V,nv+1-off,NULL,&norm,&lindep));
1404 24 : if (lindep || norm==0.0) {
1405 0 : PetscCheck(off==0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
1406 : off = 1;
1407 24 : } else PetscCall(BVScaleColumn(pjd->V,nv+1-off,1.0/norm));
1408 : }
1409 : #endif
1410 237 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1411 204 : PetscCall(BVInsertVec(pjd->W,nv,r[0]));
1412 204 : if (sz==2 && !off) PetscCall(BVInsertVec(pjd->W,nv+1,r[1]));
1413 204 : PetscCall(BVOrthogonalizeColumn(pjd->W,nv,NULL,&norm,&lindep));
1414 204 : PetscCheck(!lindep && norm>0.0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
1415 204 : PetscCall(BVScaleColumn(pjd->W,nv,1.0/norm));
1416 204 : if (sz==2 && !off) {
1417 18 : PetscCall(BVOrthogonalizeColumn(pjd->W,nv+1,NULL,&norm,&lindep));
1418 18 : PetscCheck(!lindep && norm>0.0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
1419 18 : PetscCall(BVScaleColumn(pjd->W,nv+1,1.0/norm));
1420 : }
1421 : }
1422 237 : bupdated = idx?0:nv;
1423 237 : nv += sz-off;
1424 : }
1425 2745 : for (k=0;k<nvc;k++) {
1426 2492 : eig[pep->nconv-idx+k] = pep->eigr[k];
1427 : #if !defined(PETSC_USE_COMPLEX)
1428 2492 : eigi[pep->nconv-idx+k] = pep->eigi[k];
1429 : #endif
1430 : }
1431 528 : PetscCall(PEPMonitor(pep,pep->its,pep->nconv,eig,eigi,pep->errest,pep->nconv+1));
1432 : }
1433 : }
1434 11 : if (pjd->ld>1) {
1435 44 : for (k=0;k<pep->nconv;k++) {
1436 33 : pep->eigr[k] = eig[k];
1437 33 : pep->eigi[k] = eigi[k];
1438 : }
1439 11 : if (pep->nconv>0) PetscCall(PEPJDEigenvectors(pep));
1440 11 : PetscCall(PetscFree2(pcctx->M,pcctx->ps));
1441 : }
1442 11 : PetscCall(VecDestroy(&u[0]));
1443 11 : PetscCall(VecDestroy(&r[0]));
1444 11 : PetscCall(VecDestroy(&p[0]));
1445 : #if !defined (PETSC_USE_COMPLEX)
1446 11 : PetscCall(VecDestroy(&u[1]));
1447 11 : PetscCall(VecDestroy(&r[1]));
1448 11 : PetscCall(VecDestroy(&p[1]));
1449 : #endif
1450 11 : PetscCall(KSPSetTolerances(ksp,rtol,abstol,dtol,maxits));
1451 11 : PetscCall(KSPSetPC(ksp,pcctx->pc));
1452 11 : PetscCall(VecDestroy(&pcctx->Bp[0]));
1453 11 : PetscCall(VecDestroy(&pcctx->Bp[1]));
1454 11 : PetscCall(MatShellGetContext(pjd->Pshell,&matctx));
1455 11 : PetscCall(MatDestroy(&matctx->Pr));
1456 11 : PetscCall(MatDestroy(&matctx->Pi));
1457 11 : PetscCall(MatDestroy(&pjd->Pshell));
1458 11 : PetscCall(MatDestroy(&pcctx->PPr));
1459 11 : PetscCall(PCDestroy(&pcctx->pc));
1460 11 : PetscCall(PetscFree(pcctx));
1461 11 : PetscCall(PetscFree(matctx));
1462 11 : PetscCall(PCDestroy(&pjd->pcshell));
1463 11 : PetscCall(PetscFree3(eig,eigi,res));
1464 11 : PetscCall(VecDestroy(&pjd->vtempl));
1465 11 : PetscFunctionReturn(PETSC_SUCCESS);
1466 : }
1467 :
1468 3 : static PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
1469 : {
1470 3 : PEP_JD *pjd = (PEP_JD*)pep->data;
1471 :
1472 3 : PetscFunctionBegin;
1473 3 : if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) pjd->keep = 0.5;
1474 : else {
1475 3 : PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
1476 3 : pjd->keep = keep;
1477 : }
1478 3 : PetscFunctionReturn(PETSC_SUCCESS);
1479 : }
1480 :
1481 : /*@
1482 : PEPJDSetRestart - Sets the restart parameter for the Jacobi-Davidson
1483 : method, in particular the proportion of basis vectors that must be kept
1484 : after restart.
1485 :
1486 : Logically Collective
1487 :
1488 : Input Parameters:
1489 : + pep - the eigenproblem solver context
1490 : - keep - the number of vectors to be kept at restart
1491 :
1492 : Options Database Key:
1493 : . -pep_jd_restart - Sets the restart parameter
1494 :
1495 : Notes:
1496 : Allowed values are in the range [0.1,0.9]. The default is 0.5.
1497 :
1498 : Level: advanced
1499 :
1500 : .seealso: PEPJDGetRestart()
1501 : @*/
1502 3 : PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
1503 : {
1504 3 : PetscFunctionBegin;
1505 3 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1506 9 : PetscValidLogicalCollectiveReal(pep,keep,2);
1507 3 : PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
1508 3 : PetscFunctionReturn(PETSC_SUCCESS);
1509 : }
1510 :
1511 2 : static PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
1512 : {
1513 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1514 :
1515 2 : PetscFunctionBegin;
1516 2 : *keep = pjd->keep;
1517 2 : PetscFunctionReturn(PETSC_SUCCESS);
1518 : }
1519 :
1520 : /*@
1521 : PEPJDGetRestart - Gets the restart parameter used in the Jacobi-Davidson method.
1522 :
1523 : Not Collective
1524 :
1525 : Input Parameter:
1526 : . pep - the eigenproblem solver context
1527 :
1528 : Output Parameter:
1529 : . keep - the restart parameter
1530 :
1531 : Level: advanced
1532 :
1533 : .seealso: PEPJDSetRestart()
1534 : @*/
1535 2 : PetscErrorCode PEPJDGetRestart(PEP pep,PetscReal *keep)
1536 : {
1537 2 : PetscFunctionBegin;
1538 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1539 2 : PetscAssertPointer(keep,2);
1540 2 : PetscUseMethod(pep,"PEPJDGetRestart_C",(PEP,PetscReal*),(pep,keep));
1541 2 : PetscFunctionReturn(PETSC_SUCCESS);
1542 : }
1543 :
1544 2 : static PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
1545 : {
1546 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1547 :
1548 2 : PetscFunctionBegin;
1549 2 : if (fix == (PetscReal)PETSC_DEFAULT || fix == (PetscReal)PETSC_DECIDE) pjd->fix = 0.01;
1550 : else {
1551 2 : PetscCheck(fix>=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid fix value, must be >0");
1552 2 : pjd->fix = fix;
1553 : }
1554 2 : PetscFunctionReturn(PETSC_SUCCESS);
1555 : }
1556 :
1557 : /*@
1558 : PEPJDSetFix - Sets the threshold for changing the target in the correction
1559 : equation.
1560 :
1561 : Logically Collective
1562 :
1563 : Input Parameters:
1564 : + pep - the eigenproblem solver context
1565 : - fix - threshold for changing the target
1566 :
1567 : Options Database Key:
1568 : . -pep_jd_fix - the fix value
1569 :
1570 : Note:
1571 : The target in the correction equation is fixed at the first iterations.
1572 : When the norm of the residual vector is lower than the fix value,
1573 : the target is set to the corresponding eigenvalue.
1574 :
1575 : Level: advanced
1576 :
1577 : .seealso: PEPJDGetFix()
1578 : @*/
1579 2 : PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
1580 : {
1581 2 : PetscFunctionBegin;
1582 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1583 6 : PetscValidLogicalCollectiveReal(pep,fix,2);
1584 2 : PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
1585 2 : PetscFunctionReturn(PETSC_SUCCESS);
1586 : }
1587 :
1588 2 : static PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
1589 : {
1590 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1591 :
1592 2 : PetscFunctionBegin;
1593 2 : *fix = pjd->fix;
1594 2 : PetscFunctionReturn(PETSC_SUCCESS);
1595 : }
1596 :
1597 : /*@
1598 : PEPJDGetFix - Returns the threshold for changing the target in the correction
1599 : equation.
1600 :
1601 : Not Collective
1602 :
1603 : Input Parameter:
1604 : . pep - the eigenproblem solver context
1605 :
1606 : Output Parameter:
1607 : . fix - threshold for changing the target
1608 :
1609 : Note:
1610 : The target in the correction equation is fixed at the first iterations.
1611 : When the norm of the residual vector is lower than the fix value,
1612 : the target is set to the corresponding eigenvalue.
1613 :
1614 : Level: advanced
1615 :
1616 : .seealso: PEPJDSetFix()
1617 : @*/
1618 2 : PetscErrorCode PEPJDGetFix(PEP pep,PetscReal *fix)
1619 : {
1620 2 : PetscFunctionBegin;
1621 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1622 2 : PetscAssertPointer(fix,2);
1623 2 : PetscUseMethod(pep,"PEPJDGetFix_C",(PEP,PetscReal*),(pep,fix));
1624 2 : PetscFunctionReturn(PETSC_SUCCESS);
1625 : }
1626 :
1627 1 : static PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
1628 : {
1629 1 : PEP_JD *pjd = (PEP_JD*)pep->data;
1630 :
1631 1 : PetscFunctionBegin;
1632 1 : pjd->reusepc = reusepc;
1633 1 : PetscFunctionReturn(PETSC_SUCCESS);
1634 : }
1635 :
1636 : /*@
1637 : PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
1638 : must be reused or not.
1639 :
1640 : Logically Collective
1641 :
1642 : Input Parameters:
1643 : + pep - the eigenproblem solver context
1644 : - reusepc - the reuse flag
1645 :
1646 : Options Database Key:
1647 : . -pep_jd_reuse_preconditioner - the reuse flag
1648 :
1649 : Note:
1650 : The default value is False. If set to True, the preconditioner is built
1651 : only at the beginning, using the target value. Otherwise, it may be rebuilt
1652 : (depending on the fix parameter) at each iteration from the Ritz value.
1653 :
1654 : Level: advanced
1655 :
1656 : .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
1657 : @*/
1658 1 : PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
1659 : {
1660 1 : PetscFunctionBegin;
1661 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1662 3 : PetscValidLogicalCollectiveBool(pep,reusepc,2);
1663 1 : PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
1664 1 : PetscFunctionReturn(PETSC_SUCCESS);
1665 : }
1666 :
1667 2 : static PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
1668 : {
1669 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1670 :
1671 2 : PetscFunctionBegin;
1672 2 : *reusepc = pjd->reusepc;
1673 2 : PetscFunctionReturn(PETSC_SUCCESS);
1674 : }
1675 :
1676 : /*@
1677 : PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.
1678 :
1679 : Not Collective
1680 :
1681 : Input Parameter:
1682 : . pep - the eigenproblem solver context
1683 :
1684 : Output Parameter:
1685 : . reusepc - the reuse flag
1686 :
1687 : Level: advanced
1688 :
1689 : .seealso: PEPJDSetReusePreconditioner()
1690 : @*/
1691 2 : PetscErrorCode PEPJDGetReusePreconditioner(PEP pep,PetscBool *reusepc)
1692 : {
1693 2 : PetscFunctionBegin;
1694 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1695 2 : PetscAssertPointer(reusepc,2);
1696 2 : PetscUseMethod(pep,"PEPJDGetReusePreconditioner_C",(PEP,PetscBool*),(pep,reusepc));
1697 2 : PetscFunctionReturn(PETSC_SUCCESS);
1698 : }
1699 :
1700 2 : static PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
1701 : {
1702 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1703 :
1704 2 : PetscFunctionBegin;
1705 2 : if (mmidx == PETSC_DEFAULT || mmidx == PETSC_DECIDE) {
1706 0 : if (pjd->mmidx != 1) pep->state = PEP_STATE_INITIAL;
1707 0 : pjd->mmidx = 1;
1708 : } else {
1709 2 : PetscCheck(mmidx>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mmidx value, should be >0");
1710 2 : if (pjd->mmidx != mmidx) pep->state = PEP_STATE_INITIAL;
1711 2 : pjd->mmidx = mmidx;
1712 : }
1713 2 : PetscFunctionReturn(PETSC_SUCCESS);
1714 : }
1715 :
1716 : /*@
1717 : PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.
1718 :
1719 : Logically Collective
1720 :
1721 : Input Parameters:
1722 : + pep - the eigenproblem solver context
1723 : - mmidx - maximum minimality index
1724 :
1725 : Options Database Key:
1726 : . -pep_jd_minimality_index - the minimality index value
1727 :
1728 : Note:
1729 : The default value is equal to the degree of the polynomial. A smaller value
1730 : can be used if the wanted eigenvectors are known to be linearly independent.
1731 :
1732 : Level: advanced
1733 :
1734 : .seealso: PEPJDGetMinimalityIndex()
1735 : @*/
1736 2 : PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
1737 : {
1738 2 : PetscFunctionBegin;
1739 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1740 6 : PetscValidLogicalCollectiveInt(pep,mmidx,2);
1741 2 : PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
1742 2 : PetscFunctionReturn(PETSC_SUCCESS);
1743 : }
1744 :
1745 2 : static PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
1746 : {
1747 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1748 :
1749 2 : PetscFunctionBegin;
1750 2 : *mmidx = pjd->mmidx;
1751 2 : PetscFunctionReturn(PETSC_SUCCESS);
1752 : }
1753 :
1754 : /*@
1755 : PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
1756 : index.
1757 :
1758 : Not Collective
1759 :
1760 : Input Parameter:
1761 : . pep - the eigenproblem solver context
1762 :
1763 : Output Parameter:
1764 : . mmidx - minimality index
1765 :
1766 : Level: advanced
1767 :
1768 : .seealso: PEPJDSetMinimalityIndex()
1769 : @*/
1770 2 : PetscErrorCode PEPJDGetMinimalityIndex(PEP pep,PetscInt *mmidx)
1771 : {
1772 2 : PetscFunctionBegin;
1773 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1774 2 : PetscAssertPointer(mmidx,2);
1775 2 : PetscUseMethod(pep,"PEPJDGetMinimalityIndex_C",(PEP,PetscInt*),(pep,mmidx));
1776 2 : PetscFunctionReturn(PETSC_SUCCESS);
1777 : }
1778 :
1779 2 : static PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
1780 : {
1781 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1782 :
1783 2 : PetscFunctionBegin;
1784 2 : switch (proj) {
1785 2 : case PEP_JD_PROJECTION_HARMONIC:
1786 : case PEP_JD_PROJECTION_ORTHOGONAL:
1787 2 : if (pjd->proj != proj) {
1788 2 : pep->state = PEP_STATE_INITIAL;
1789 2 : pjd->proj = proj;
1790 : }
1791 2 : break;
1792 0 : default:
1793 0 : SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'proj' value");
1794 : }
1795 2 : PetscFunctionReturn(PETSC_SUCCESS);
1796 : }
1797 :
1798 : /*@
1799 : PEPJDSetProjection - Sets the type of projection to be used in the Jacobi-Davidson solver.
1800 :
1801 : Logically Collective
1802 :
1803 : Input Parameters:
1804 : + pep - the eigenproblem solver context
1805 : - proj - the type of projection
1806 :
1807 : Options Database Key:
1808 : . -pep_jd_projection - the projection type, either orthogonal or harmonic
1809 :
1810 : Level: advanced
1811 :
1812 : .seealso: PEPJDGetProjection()
1813 : @*/
1814 2 : PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
1815 : {
1816 2 : PetscFunctionBegin;
1817 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1818 6 : PetscValidLogicalCollectiveEnum(pep,proj,2);
1819 2 : PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
1820 2 : PetscFunctionReturn(PETSC_SUCCESS);
1821 : }
1822 :
1823 2 : static PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
1824 : {
1825 2 : PEP_JD *pjd = (PEP_JD*)pep->data;
1826 :
1827 2 : PetscFunctionBegin;
1828 2 : *proj = pjd->proj;
1829 2 : PetscFunctionReturn(PETSC_SUCCESS);
1830 : }
1831 :
1832 : /*@
1833 : PEPJDGetProjection - Returns the type of projection used by the Jacobi-Davidson solver.
1834 :
1835 : Not Collective
1836 :
1837 : Input Parameter:
1838 : . pep - the eigenproblem solver context
1839 :
1840 : Output Parameter:
1841 : . proj - the type of projection
1842 :
1843 : Level: advanced
1844 :
1845 : .seealso: PEPJDSetProjection()
1846 : @*/
1847 2 : PetscErrorCode PEPJDGetProjection(PEP pep,PEPJDProjection *proj)
1848 : {
1849 2 : PetscFunctionBegin;
1850 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1851 2 : PetscAssertPointer(proj,2);
1852 2 : PetscUseMethod(pep,"PEPJDGetProjection_C",(PEP,PEPJDProjection*),(pep,proj));
1853 2 : PetscFunctionReturn(PETSC_SUCCESS);
1854 : }
1855 :
1856 10 : static PetscErrorCode PEPSetFromOptions_JD(PEP pep,PetscOptionItems *PetscOptionsObject)
1857 : {
1858 10 : PetscBool flg,b1;
1859 10 : PetscReal r1;
1860 10 : PetscInt i1;
1861 10 : PEPJDProjection proj;
1862 :
1863 10 : PetscFunctionBegin;
1864 10 : PetscOptionsHeadBegin(PetscOptionsObject,"PEP JD Options");
1865 :
1866 10 : PetscCall(PetscOptionsReal("-pep_jd_restart","Proportion of vectors kept after restart","PEPJDSetRestart",0.5,&r1,&flg));
1867 10 : if (flg) PetscCall(PEPJDSetRestart(pep,r1));
1868 :
1869 10 : PetscCall(PetscOptionsReal("-pep_jd_fix","Tolerance for changing the target in the correction equation","PEPJDSetFix",0.01,&r1,&flg));
1870 10 : if (flg) PetscCall(PEPJDSetFix(pep,r1));
1871 :
1872 10 : PetscCall(PetscOptionsBool("-pep_jd_reuse_preconditioner","Whether to reuse the preconditioner","PEPJDSetReusePreconditoiner",PETSC_FALSE,&b1,&flg));
1873 10 : if (flg) PetscCall(PEPJDSetReusePreconditioner(pep,b1));
1874 :
1875 10 : PetscCall(PetscOptionsInt("-pep_jd_minimality_index","Maximum allowed minimality index","PEPJDSetMinimalityIndex",1,&i1,&flg));
1876 10 : if (flg) PetscCall(PEPJDSetMinimalityIndex(pep,i1));
1877 :
1878 10 : PetscCall(PetscOptionsEnum("-pep_jd_projection","Type of projection","PEPJDSetProjection",PEPJDProjectionTypes,(PetscEnum)PEP_JD_PROJECTION_HARMONIC,(PetscEnum*)&proj,&flg));
1879 10 : if (flg) PetscCall(PEPJDSetProjection(pep,proj));
1880 :
1881 10 : PetscOptionsHeadEnd();
1882 10 : PetscFunctionReturn(PETSC_SUCCESS);
1883 : }
1884 :
1885 0 : static PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
1886 : {
1887 0 : PEP_JD *pjd = (PEP_JD*)pep->data;
1888 0 : PetscBool isascii;
1889 :
1890 0 : PetscFunctionBegin;
1891 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1892 0 : if (isascii) {
1893 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*pjd->keep)));
1894 0 : PetscCall(PetscViewerASCIIPrintf(viewer," threshold for changing the target in the correction equation (fix): %g\n",(double)pjd->fix));
1895 0 : PetscCall(PetscViewerASCIIPrintf(viewer," projection type: %s\n",PEPJDProjectionTypes[pjd->proj]));
1896 0 : PetscCall(PetscViewerASCIIPrintf(viewer," maximum allowed minimality index: %" PetscInt_FMT "\n",pjd->mmidx));
1897 0 : if (pjd->reusepc) PetscCall(PetscViewerASCIIPrintf(viewer," reusing the preconditioner\n"));
1898 : }
1899 0 : PetscFunctionReturn(PETSC_SUCCESS);
1900 : }
1901 :
1902 21 : static PetscErrorCode PEPSetDefaultST_JD(PEP pep)
1903 : {
1904 21 : KSP ksp;
1905 :
1906 21 : PetscFunctionBegin;
1907 21 : if (!((PetscObject)pep->st)->type_name) {
1908 10 : PetscCall(STSetType(pep->st,STPRECOND));
1909 10 : PetscCall(STPrecondSetKSPHasMat(pep->st,PETSC_TRUE));
1910 : }
1911 21 : PetscCall(STSetTransform(pep->st,PETSC_FALSE));
1912 21 : PetscCall(STGetKSP(pep->st,&ksp));
1913 21 : if (!((PetscObject)ksp)->type_name) {
1914 10 : PetscCall(KSPSetType(ksp,KSPBCGSL));
1915 10 : PetscCall(KSPSetTolerances(ksp,1e-5,PETSC_CURRENT,PETSC_CURRENT,100));
1916 : }
1917 21 : PetscFunctionReturn(PETSC_SUCCESS);
1918 : }
1919 :
1920 12 : static PetscErrorCode PEPReset_JD(PEP pep)
1921 : {
1922 12 : PEP_JD *pjd = (PEP_JD*)pep->data;
1923 12 : PetscInt i;
1924 :
1925 12 : PetscFunctionBegin;
1926 47 : for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->TV+i));
1927 12 : if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVDestroy(&pjd->W));
1928 12 : if (pjd->ld>1) {
1929 12 : PetscCall(BVDestroy(&pjd->V));
1930 47 : for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->AX+i));
1931 12 : PetscCall(BVDestroy(&pjd->N[0]));
1932 12 : PetscCall(BVDestroy(&pjd->N[1]));
1933 12 : PetscCall(PetscFree3(pjd->XpX,pjd->T,pjd->Tj));
1934 : }
1935 12 : PetscCall(PetscFree2(pjd->TV,pjd->AX));
1936 12 : PetscFunctionReturn(PETSC_SUCCESS);
1937 : }
1938 :
1939 10 : static PetscErrorCode PEPDestroy_JD(PEP pep)
1940 : {
1941 10 : PetscFunctionBegin;
1942 10 : PetscCall(PetscFree(pep->data));
1943 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",NULL));
1944 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",NULL));
1945 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",NULL));
1946 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",NULL));
1947 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",NULL));
1948 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",NULL));
1949 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",NULL));
1950 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",NULL));
1951 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",NULL));
1952 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",NULL));
1953 10 : PetscFunctionReturn(PETSC_SUCCESS);
1954 : }
1955 :
1956 10 : SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
1957 : {
1958 10 : PEP_JD *pjd;
1959 :
1960 10 : PetscFunctionBegin;
1961 10 : PetscCall(PetscNew(&pjd));
1962 10 : pep->data = (void*)pjd;
1963 :
1964 10 : pep->lineariz = PETSC_FALSE;
1965 10 : pjd->fix = 0.01;
1966 10 : pjd->mmidx = 0;
1967 :
1968 10 : pep->ops->solve = PEPSolve_JD;
1969 10 : pep->ops->setup = PEPSetUp_JD;
1970 10 : pep->ops->setfromoptions = PEPSetFromOptions_JD;
1971 10 : pep->ops->destroy = PEPDestroy_JD;
1972 10 : pep->ops->reset = PEPReset_JD;
1973 10 : pep->ops->view = PEPView_JD;
1974 10 : pep->ops->setdefaultst = PEPSetDefaultST_JD;
1975 :
1976 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",PEPJDSetRestart_JD));
1977 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",PEPJDGetRestart_JD));
1978 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",PEPJDSetFix_JD));
1979 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",PEPJDGetFix_JD));
1980 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",PEPJDSetReusePreconditioner_JD));
1981 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",PEPJDGetReusePreconditioner_JD));
1982 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",PEPJDSetMinimalityIndex_JD));
1983 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",PEPJDGetMinimalityIndex_JD));
1984 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",PEPJDSetProjection_JD));
1985 10 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",PEPJDGetProjection_JD));
1986 10 : PetscFunctionReturn(PETSC_SUCCESS);
1987 : }
|
