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