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