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