Actual source code: pjd.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc polynomial eigensolver: "jd"

 13:    Method: Jacobi-Davidson

 15:    Algorithm:

 17:        Jacobi-Davidson for polynomial eigenvalue problems.

 19:    References:

 21:        [1] C. Campos and J.E. Roman, "A polynomial Jacobi-Davidson solver
 22:            with support for non-monomial bases and deflation", BIT Numer.
 23:            Math. 60:295-318, 2020.

 25:        [2] G.L.G. Sleijpen et al., "Jacobi-Davidson type methods for
 26:            generalized eigenproblems and polynomial eigenproblems", BIT
 27:            36(3):595-633, 1996.

 29:        [3] Feng-Nan Hwang, Zih-Hao Wei, Tsung-Ming Huang, Weichung Wang,
 30:            "A Parallel Additive Schwarz Preconditioned Jacobi-Davidson
 31:            Algorithm for Polynomial Eigenvalue Problems in Quantum Dot
 32:            Simulation", J. Comput. Phys. 229(8):2932-2947, 2010.
 33: */

 35: #include <slepc/private/pepimpl.h>
 36: #include <slepcblaslapack.h>

 38: static PetscBool  cited = PETSC_FALSE;
 39: static const char citation[] =
 40:   "@Article{slepc-slice-qep,\n"
 41:   "   author = \"C. Campos and J. E. Roman\",\n"
 42:   "   title = \"A polynomial {Jacobi-Davidson} solver with support for non-monomial bases and deflation\",\n"
 43:   "   journal = \"{BIT} Numer. Math.\",\n"
 44:   "   volume = \"60\",\n"
 45:   "   pages = \"295--318\",\n"
 46:   "   year = \"2020,\"\n"
 47:   "   doi = \"https://doi.org/10.1007/s10543-019-00778-z\"\n"
 48:   "}\n";

 50: typedef struct {
 51:   PetscReal   keep;          /* restart parameter */
 52:   PetscReal   fix;           /* fix parameter */
 53:   PetscBool   reusepc;       /* flag indicating whether pc is rebuilt or not */
 54:   BV          V;             /* work basis vectors to store the search space */
 55:   BV          W;             /* work basis vectors to store the test space */
 56:   BV          *TV;           /* work basis vectors to store T*V (each TV[i] is the coefficient for \lambda^i of T*V for the extended T) */
 57:   BV          *AX;           /* work basis vectors to store A_i*X for locked eigenvectors */
 58:   BV          N[2];          /* auxiliary work BVs */
 59:   BV          X;             /* locked eigenvectors */
 60:   PetscScalar *T;            /* matrix of the invariant pair */
 61:   PetscScalar *Tj;           /* matrix containing the powers of the invariant pair matrix */
 62:   PetscScalar *XpX;          /* X^H*X */
 63:   PetscInt    ld;            /* leading dimension for Tj and XpX */
 64:   PC          pcshell;       /* preconditioner including basic precond+projector */
 65:   Mat         Pshell;        /* auxiliary shell matrix */
 66:   PetscInt    nlock;         /* number of locked vectors in the invariant pair */
 67:   Vec         vtempl;        /* reference nested vector */
 68:   PetscInt    midx;          /* minimality index */
 69:   PetscInt    mmidx;         /* maximum allowed minimality index */
 70:   PEPJDProjection proj;      /* projection type (orthogonal, harmonic) */
 71: } PEP_JD;

 73: typedef struct {
 74:   PEP         pep;
 75:   PC          pc;            /* basic preconditioner */
 76:   Vec         Bp[2];         /* preconditioned residual of derivative polynomial, B\p */
 77:   Vec         u[2];          /* Ritz vector */
 78:   PetscScalar gamma[2];      /* precomputed scalar u'*B\p */
 79:   PetscScalar theta;
 80:   PetscScalar *M;
 81:   PetscScalar *ps;
 82:   PetscInt    ld;
 83:   Vec         *work;
 84:   Mat         PPr;
 85:   BV          X;
 86:   PetscInt    n;
 87: } PEP_JD_PCSHELL;

 89: typedef struct {
 90:   Mat         Pr,Pi;         /* matrix polynomial evaluated at theta */
 91:   PEP         pep;
 92:   Vec         *work;
 93:   PetscScalar theta[2];
 94: } PEP_JD_MATSHELL;

 96: /*
 97:    Duplicate and resize auxiliary basis
 98: */
 99: static PetscErrorCode PEPJDDuplicateBasis(PEP pep,BV *basis)
100: {
101:   PEP_JD             *pjd = (PEP_JD*)pep->data;
102:   PetscInt           nloc,m;
103:   BVType             type;
104:   BVOrthogType       otype;
105:   BVOrthogRefineType oref;
106:   PetscReal          oeta;
107:   BVOrthogBlockType  oblock;

109:   if (pjd->ld>1) {
110:     BVCreate(PetscObjectComm((PetscObject)pep),basis);
111:     BVGetSizes(pep->V,&nloc,NULL,&m);
112:     nloc += pjd->ld-1;
113:     BVSetSizes(*basis,nloc,PETSC_DECIDE,m);
114:     BVGetType(pep->V,&type);
115:     BVSetType(*basis,type);
116:     BVGetOrthogonalization(pep->V,&otype,&oref,&oeta,&oblock);
117:     BVSetOrthogonalization(*basis,otype,oref,oeta,oblock);
118:     PetscObjectStateIncrease((PetscObject)*basis);
119:   } else BVDuplicate(pep->V,basis);
120:   return 0;
121: }

123: PetscErrorCode PEPSetUp_JD(PEP pep)
124: {
125:   PEP_JD         *pjd = (PEP_JD*)pep->data;
126:   PetscBool      isprecond,flg;
127:   PetscRandom    rand;
128:   PetscInt       i;

130:   PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
131:   if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
132:   if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;

135:   PetscObjectTypeCompare((PetscObject)pep->st,STPRECOND,&isprecond);

138:   STGetTransform(pep->st,&flg);
140:   PEPCheckIgnored(pep,PEP_FEATURE_EXTRACT);

142:   if (!pjd->mmidx) pjd->mmidx = pep->nmat-1;
143:   pjd->mmidx = PetscMin(pjd->mmidx,pep->nmat-1);
144:   if (!pjd->keep) pjd->keep = 0.5;
145:   PEPBasisCoefficients(pep,pep->pbc);
146:   PEPAllocateSolution(pep,0);
147:   BVGetRandomContext(pep->V,&rand);  /* make sure the random context is available when duplicating */
148:   PEPSetWorkVecs(pep,5);
149:   pjd->ld = pep->nev;
150: #if !defined (PETSC_USE_COMPLEX)
151:   pjd->ld++;
152: #endif
153:   PetscMalloc2(pep->nmat,&pjd->TV,pep->nmat,&pjd->AX);
154:   for (i=0;i<pep->nmat;i++) PEPJDDuplicateBasis(pep,pjd->TV+i);
155:   if (pjd->ld>1) {
156:     PEPJDDuplicateBasis(pep,&pjd->V);
157:     BVSetFromOptions(pjd->V);
158:     for (i=0;i<pep->nmat;i++) BVDuplicateResize(pep->V,pjd->ld-1,pjd->AX+i);
159:     BVDuplicateResize(pep->V,pjd->ld-1,pjd->N);
160:     BVDuplicateResize(pep->V,pjd->ld-1,pjd->N+1);
161:     pjd->X = pep->V;
162:     PetscCalloc3((pjd->ld)*(pjd->ld),&pjd->XpX,pep->ncv*pep->ncv,&pjd->T,pjd->ld*pjd->ld*pep->nmat,&pjd->Tj);
163:   } else pjd->V = pep->V;
164:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PEPJDDuplicateBasis(pep,&pjd->W);
165:   else pjd->W = pjd->V;
166:   DSSetType(pep->ds,DSPEP);
167:   DSPEPSetDegree(pep->ds,pep->nmat-1);
168:   if (pep->basis!=PEP_BASIS_MONOMIAL) DSPEPSetCoefficients(pep->ds,pep->pbc);
169:   DSAllocate(pep->ds,pep->ncv);
170:   return 0;
171: }

173: /*
174:    Updates columns (low to (high-1)) of TV[i]
175: */
176: static PetscErrorCode PEPJDUpdateTV(PEP pep,PetscInt low,PetscInt high,Vec *w)
177: {
178:   PEP_JD         *pjd = (PEP_JD*)pep->data;
179:   PetscInt       pp,col,i,nloc,nconv;
180:   Vec            v1,v2,t1,t2;
181:   PetscScalar    *array1,*array2,*x2,*xx,*N,*Np,*y2=NULL,zero=0.0,sone=1.0,*pT,fact,*psc;
182:   PetscReal      *cg,*ca,*cb;
183:   PetscMPIInt    rk,np;
184:   PetscBLASInt   n_,ld_,one=1;
185:   Mat            T;
186:   BV             pbv;

188:   ca = pep->pbc; cb = ca+pep->nmat; cg = cb + pep->nmat;
189:   nconv = pjd->nlock;
190:   PetscMalloc5(nconv,&x2,nconv,&xx,nconv*nconv,&pT,nconv*nconv,&N,nconv*nconv,&Np);
191:   MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk);
192:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
193:   BVGetSizes(pep->V,&nloc,NULL,NULL);
194:   t1 = w[0];
195:   t2 = w[1];
196:   PetscBLASIntCast(pjd->nlock,&n_);
197:   PetscBLASIntCast(pjd->ld,&ld_);
198:   if (nconv) {
199:     for (i=0;i<nconv;i++) PetscArraycpy(pT+i*nconv,pjd->T+i*pep->ncv,nconv);
200:     MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,pT,&T);
201:   }
202:   for (col=low;col<high;col++) {
203:     BVGetColumn(pjd->V,col,&v1);
204:     VecGetArray(v1,&array1);
205:     if (nconv>0) {
206:       for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
207:     }
208:     VecPlaceArray(t1,array1);
209:     if (nconv) {
210:       BVSetActiveColumns(pjd->N[0],0,nconv);
211:       BVSetActiveColumns(pjd->N[1],0,nconv);
212:       BVDotVec(pjd->X,t1,xx);
213:     }
214:     for (pp=pep->nmat-1;pp>=0;pp--) {
215:       BVGetColumn(pjd->TV[pp],col,&v2);
216:       VecGetArray(v2,&array2);
217:       VecPlaceArray(t2,array2);
218:       MatMult(pep->A[pp],t1,t2);
219:       if (nconv) {
220:         if (pp<pep->nmat-3) {
221:           BVMult(pjd->N[0],1.0,-cg[pp+2],pjd->AX[pp+1],NULL);
222:           MatShift(T,-cb[pp+1]);
223:           BVMult(pjd->N[0],1.0/ca[pp],1.0/ca[pp],pjd->N[1],T);
224:           pbv = pjd->N[0]; pjd->N[0] = pjd->N[1]; pjd->N[1] = pbv;
225:           BVMultVec(pjd->N[1],1.0,1.0,t2,x2);
226:           MatShift(T,cb[pp+1]);
227:         } else if (pp==pep->nmat-3) {
228:           BVCopy(pjd->AX[pp+2],pjd->N[0]);
229:           BVScale(pjd->N[0],1/ca[pp+1]);
230:           BVCopy(pjd->AX[pp+1],pjd->N[1]);
231:           MatShift(T,-cb[pp+1]);
232:           BVMult(pjd->N[1],1.0/ca[pp],1.0/ca[pp],pjd->N[0],T);
233:           BVMultVec(pjd->N[1],1.0,1.0,t2,x2);
234:           MatShift(T,cb[pp+1]);
235:         } else if (pp==pep->nmat-2) BVMultVec(pjd->AX[pp+1],1.0/ca[pp],1.0,t2,x2);
236:         if (pp<pjd->midx) {
237:           y2 = array2+nloc;
238:           PetscCallBLAS("BLASgemv",BLASgemv_("C",&n_,&n_,&sone,pjd->Tj+pjd->ld*pjd->ld*pp,&ld_,xx,&one,&zero,y2,&one));
239:           if (pp<pjd->midx-2) {
240:             fact = -cg[pp+2];
241:             PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&fact,Np,&n_));
242:             fact = 1/ca[pp];
243:             MatShift(T,-cb[pp+1]);
244:             PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&fact,N,&n_,pT,&n_,&fact,Np,&n_));
245:             MatShift(T,cb[pp+1]);
246:             psc = Np; Np = N; N = psc;
247:             PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
248:           } else if (pp==pjd->midx-2) {
249:             fact = 1/ca[pp];
250:             PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&fact,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&zero,N,&n_));
251:             PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
252:           } else if (pp==pjd->midx-1) PetscArrayzero(Np,nconv*nconv);
253:         }
254:         for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
255:       }
256:       VecResetArray(t2);
257:       VecRestoreArray(v2,&array2);
258:       BVRestoreColumn(pjd->TV[pp],col,&v2);
259:     }
260:     VecResetArray(t1);
261:     VecRestoreArray(v1,&array1);
262:     BVRestoreColumn(pjd->V,col,&v1);
263:   }
264:   if (nconv) MatDestroy(&T);
265:   PetscFree5(x2,xx,pT,N,Np);
266:   return 0;
267: }

269: /*
270:    RRQR of X. Xin*P=Xou*R. Rank of R is rk
271: */
272: static PetscErrorCode PEPJDOrthogonalize(PetscInt row,PetscInt col,PetscScalar *X,PetscInt ldx,PetscInt *rk,PetscInt *P,PetscScalar *R,PetscInt ldr)
273: {
274:   PetscInt       i,j,n,r;
275:   PetscBLASInt   row_,col_,ldx_,*p,lwork,info,n_;
276:   PetscScalar    *tau,*work;
277:   PetscReal      tol,*rwork;

279:   PetscBLASIntCast(row,&row_);
280:   PetscBLASIntCast(col,&col_);
281:   PetscBLASIntCast(ldx,&ldx_);
282:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
283:   n = PetscMin(row,col);
284:   PetscBLASIntCast(n,&n_);
285:   lwork = 3*col_+1;
286:   PetscMalloc4(col,&p,n,&tau,lwork,&work,2*col,&rwork);
287:   for (i=1;i<col;i++) p[i] = 0;
288:   p[0] = 1;

290:   /* rank revealing QR */
291: #if defined(PETSC_USE_COMPLEX)
292:   PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,rwork,&info));
293: #else
294:   PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,&info));
295: #endif
296:   SlepcCheckLapackInfo("geqp3",info);
297:   if (P) for (i=0;i<col;i++) P[i] = p[i]-1;

299:   /* rank computation */
300:   tol = PetscMax(row,col)*PETSC_MACHINE_EPSILON*PetscAbsScalar(X[0]);
301:   r = 1;
302:   for (i=1;i<n;i++) {
303:     if (PetscAbsScalar(X[i+ldx*i])>tol) r++;
304:     else break;
305:   }
306:   if (rk) *rk=r;

308:   /* copy upper triangular matrix if requested */
309:   if (R) {
310:      for (i=0;i<r;i++) {
311:        PetscArrayzero(R+i*ldr,r);
312:        for (j=0;j<=i;j++) R[i*ldr+j] = X[i*ldx+j];
313:      }
314:   }
315:   PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&row_,&n_,&n_,X,&ldx_,tau,work,&lwork,&info));
316:   SlepcCheckLapackInfo("orgqr",info);
317:   PetscFPTrapPop();
318:   PetscFree4(p,tau,work,rwork);
319:   return 0;
320: }

322: /*
323:    Application of extended preconditioner
324: */
325: static PetscErrorCode PEPJDExtendedPCApply(PC pc,Vec x,Vec y)
326: {
327:   PetscInt          i,j,nloc,n,ld=0;
328:   PetscMPIInt       np;
329:   Vec               tx,ty;
330:   PEP_JD_PCSHELL    *ctx;
331:   const PetscScalar *array1;
332:   PetscScalar       *x2=NULL,*t=NULL,*ps=NULL,*array2,zero=0.0,sone=1.0;
333:   PetscBLASInt      one=1,ld_,n_,ncv_;
334:   PEP_JD            *pjd=NULL;

336:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&np);
337:   PCShellGetContext(pc,&ctx);
338:   n  = ctx->n;
339:   if (n) {
340:     pjd = (PEP_JD*)ctx->pep->data;
341:     ps = ctx->ps;
342:     ld = pjd->ld;
343:     PetscMalloc2(n,&x2,n,&t);
344:     VecGetLocalSize(ctx->work[0],&nloc);
345:     VecGetArrayRead(x,&array1);
346:     for (i=0;i<n;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
347:     VecRestoreArrayRead(x,&array1);
348:   }

350:   /* y = B\x apply PC */
351:   tx = ctx->work[0];
352:   ty = ctx->work[1];
353:   VecGetArrayRead(x,&array1);
354:   VecPlaceArray(tx,array1);
355:   VecGetArray(y,&array2);
356:   VecPlaceArray(ty,array2);
357:   PCApply(ctx->pc,tx,ty);
358:   if (n) {
359:     PetscBLASIntCast(ld,&ld_);
360:     PetscBLASIntCast(n,&n_);
361:     for (i=0;i<n;i++) {
362:       t[i] = 0.0;
363:       for (j=0;j<n;j++) t[i] += ctx->M[i+j*ld]*x2[j];
364:     }
365:     if (pjd->midx==1) {
366:       PetscBLASIntCast(ctx->pep->ncv,&ncv_);
367:       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] -= ctx->theta;
368:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,pjd->T,&ncv_,t,&one,&zero,x2,&one));
369:       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] += ctx->theta;
370:       for (i=0;i<n;i++) array2[nloc+i] = x2[i];
371:       for (i=0;i<n;i++) x2[i] = -t[i];
372:     } else {
373:       for (i=0;i<n;i++) array2[nloc+i] = t[i];
374:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,ps,&ld_,t,&one,&zero,x2,&one));
375:     }
376:     for (i=0;i<n;i++) array2[nloc+i] /= PetscSqrtReal(np);
377:     BVSetActiveColumns(pjd->X,0,n);
378:     BVMultVec(pjd->X,-1.0,1.0,ty,x2);
379:     PetscFree2(x2,t);
380:   }
381:   VecResetArray(tx);
382:   VecResetArray(ty);
383:   VecRestoreArrayRead(x,&array1);
384:   VecRestoreArray(y,&array2);
385:   return 0;
386: }

388: /*
389:    Application of shell preconditioner:
390:       y = B\x - eta*B\p,  with eta = (u'*B\x)/(u'*B\p)
391: */
392: static PetscErrorCode PCShellApply_PEPJD(PC pc,Vec x,Vec y)
393: {
394:   PetscScalar    rr,eta;
395:   PEP_JD_PCSHELL *ctx;
396:   PetscInt       sz;
397:   const Vec      *xs,*ys;
398: #if !defined(PETSC_USE_COMPLEX)
399:   PetscScalar    rx,xr,xx;
400: #endif

402:   PCShellGetContext(pc,&ctx);
403:   VecCompGetSubVecs(x,&sz,&xs);
404:   VecCompGetSubVecs(y,NULL,&ys);
405:   /* y = B\x apply extended PC */
406:   PEPJDExtendedPCApply(pc,xs[0],ys[0]);
407: #if !defined(PETSC_USE_COMPLEX)
408:   if (sz==2) PEPJDExtendedPCApply(pc,xs[1],ys[1]);
409: #endif

411:   /* Compute eta = u'*y / u'*Bp */
412:   VecDot(ys[0],ctx->u[0],&rr);
413:   eta  = -rr*ctx->gamma[0];
414: #if !defined(PETSC_USE_COMPLEX)
415:   if (sz==2) {
416:     VecDot(ys[0],ctx->u[1],&xr);
417:     VecDot(ys[1],ctx->u[0],&rx);
418:     VecDot(ys[1],ctx->u[1],&xx);
419:     eta += -ctx->gamma[0]*xx-ctx->gamma[1]*(-xr+rx);
420:   }
421: #endif
422:   eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];

424:   /* y = y - eta*Bp */
425:   VecAXPY(ys[0],eta,ctx->Bp[0]);
426: #if !defined(PETSC_USE_COMPLEX)
427:   if (sz==2) {
428:     VecAXPY(ys[1],eta,ctx->Bp[1]);
429:     eta = -ctx->gamma[1]*(rr+xx)+ctx->gamma[0]*(-xr+rx);
430:     eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
431:     VecAXPY(ys[0],eta,ctx->Bp[1]);
432:     VecAXPY(ys[1],-eta,ctx->Bp[0]);
433:   }
434: #endif
435:   return 0;
436: }

438: static PetscErrorCode PEPJDCopyToExtendedVec(PEP pep,Vec v,PetscScalar *a,PetscInt na,PetscInt off,Vec vex,PetscBool back)
439: {
440:   PetscMPIInt    np,rk,count;
441:   PetscScalar    *array1,*array2;
442:   PetscInt       nloc;

444:   MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk);
445:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
446:   BVGetSizes(pep->V,&nloc,NULL,NULL);
447:   if (v) {
448:     VecGetArray(v,&array1);
449:     VecGetArray(vex,&array2);
450:     if (back) PetscArraycpy(array1,array2,nloc);
451:     else PetscArraycpy(array2,array1,nloc);
452:     VecRestoreArray(v,&array1);
453:     VecRestoreArray(vex,&array2);
454:   }
455:   if (a) {
456:     VecGetArray(vex,&array2);
457:     if (back) {
458:       PetscArraycpy(a,array2+nloc+off,na);
459:       PetscMPIIntCast(na,&count);
460:       MPI_Bcast(a,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
461:     } else {
462:       PetscArraycpy(array2+nloc+off,a,na);
463:       PetscMPIIntCast(na,&count);
464:       MPI_Bcast(array2+nloc+off,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
465:     }
466:     VecRestoreArray(vex,&array2);
467:   }
468:   return 0;
469: }

471: /* Computes Phi^hat(lambda) times a vector or its derivative (depends on beval)
472:      if no vector is provided returns a matrix
473:  */
474: static PetscErrorCode PEPJDEvaluateHatBasis(PEP pep,PetscInt n,PetscScalar *H,PetscInt ldh,PetscScalar *beval,PetscScalar *t,PetscInt idx,PetscScalar *qpp,PetscScalar *qp,PetscScalar *q)
475: {
476:   PetscInt       j,i;
477:   PetscBLASInt   n_,ldh_,one=1;
478:   PetscReal      *a,*b,*g;
479:   PetscScalar    sone=1.0,zero=0.0;

481:   a = pep->pbc; b=a+pep->nmat; g=b+pep->nmat;
482:   PetscBLASIntCast(n,&n_);
483:   PetscBLASIntCast(ldh,&ldh_);
484:   if (idx<1) PetscArrayzero(q,t?n:n*n);
485:   else if (idx==1) {
486:     if (t) {for (j=0;j<n;j++) q[j] = t[j]*beval[idx-1]/a[0];}
487:     else {
488:       PetscArrayzero(q,n*n);
489:       for (j=0;j<n;j++) q[(j+1)*n] = beval[idx-1]/a[0];
490:     }
491:   } else {
492:     if (t) {
493:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,H,&ldh_,qp,&one,&zero,q,&one));
494:       for (j=0;j<n;j++) {
495:         q[j] += beval[idx-1]*t[j]-b[idx-1]*qp[j]-g[idx-1]*qpp[j];
496:         q[j] /= a[idx-1];
497:       }
498:     } else {
499:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,H,&ldh_,qp,&n_,&zero,q,&n_));
500:       for (j=0;j<n;j++) {
501:         q[j+n*j] += beval[idx-1];
502:         for (i=0;i<n;i++) {
503:           q[i+n*j] += -b[idx-1]*qp[j*n+i]-g[idx-1]*qpp[j*n+i];
504:           q[i+n*j] /= a[idx-1];
505:         }
506:       }
507:     }
508:   }
509:   return 0;
510: }

512: static PetscErrorCode PEPJDComputeResidual(PEP pep,PetscBool derivative,PetscInt sz,Vec *u,PetscScalar *theta,Vec *p,Vec *work)
513: {
514:   PEP_JD         *pjd = (PEP_JD*)pep->data;
515:   PetscMPIInt    rk,np,count;
516:   Vec            tu,tp,w;
517:   PetscScalar    *dval,*dvali,*array1,*array2,*x2=NULL,*y2,*qj=NULL,*tt=NULL,*xx=NULL,*xxi=NULL,sone=1.0;
518:   PetscInt       i,j,nconv,nloc;
519:   PetscBLASInt   n,ld,one=1;
520: #if !defined(PETSC_USE_COMPLEX)
521:   Vec            tui=NULL,tpi=NULL;
522:   PetscScalar    *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
523: #endif

525:   nconv = pjd->nlock;
526:   if (!nconv) PetscMalloc1(2*sz*pep->nmat,&dval);
527:   else {
528:     PetscMalloc5(2*pep->nmat,&dval,2*nconv,&xx,nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*pep->nmat,&qj);
529:     MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk);
530:     MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
531:     BVGetSizes(pep->V,&nloc,NULL,NULL);
532:     VecGetArray(u[0],&array1);
533:     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]*PetscSqrtReal(np);
534:     VecRestoreArray(u[0],&array1);
535: #if !defined(PETSC_USE_COMPLEX)
536:     if (sz==2) {
537:       x2i = x2+nconv;
538:       VecGetArray(u[1],&arrayi1);
539:       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]*PetscSqrtReal(np);
540:       VecRestoreArray(u[1],&arrayi1);
541:     }
542: #endif
543:   }
544:   dvali = dval+pep->nmat;
545:   tu = work[0];
546:   tp = work[1];
547:   w  = work[2];
548:   VecGetArray(u[0],&array1);
549:   VecPlaceArray(tu,array1);
550:   VecGetArray(p[0],&array2);
551:   VecPlaceArray(tp,array2);
552:   VecSet(tp,0.0);
553: #if !defined(PETSC_USE_COMPLEX)
554:   if (sz==2) {
555:     tui = work[3];
556:     tpi = work[4];
557:     VecGetArray(u[1],&arrayi1);
558:     VecPlaceArray(tui,arrayi1);
559:     VecGetArray(p[1],&arrayi2);
560:     VecPlaceArray(tpi,arrayi2);
561:     VecSet(tpi,0.0);
562:   }
563: #endif
564:   if (derivative) PEPEvaluateBasisDerivative(pep,theta[0],theta[1],dval,dvali);
565:   else PEPEvaluateBasis(pep,theta[0],theta[1],dval,dvali);
566:   for (i=derivative?1:0;i<pep->nmat;i++) {
567:     MatMult(pep->A[i],tu,w);
568:     VecAXPY(tp,dval[i],w);
569: #if !defined(PETSC_USE_COMPLEX)
570:     if (sz==2) {
571:       VecAXPY(tpi,dvali[i],w);
572:       MatMult(pep->A[i],tui,w);
573:       VecAXPY(tpi,dval[i],w);
574:       VecAXPY(tp,-dvali[i],w);
575:     }
576: #endif
577:   }
578:   if (nconv) {
579:     for (i=0;i<pep->nmat;i++) 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);
580: #if !defined(PETSC_USE_COMPLEX)
581:     if (sz==2) {
582:       qji = qj+nconv*pep->nmat;
583:       qq = qji+nconv*pep->nmat;
584:       for (i=0;i<pep->nmat;i++) 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);
585:       for (i=0;i<nconv*pep->nmat;i++) qj[i] -= qji[i];
586:       for (i=0;i<pep->nmat;i++) {
587:         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);
588:         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);
589:       }
590:       for (i=0;i<nconv*pep->nmat;i++) qji[i] += qq[i];
591:       for (i=derivative?2:1;i<pep->nmat;i++) BVMultVec(pjd->AX[i],1.0,1.0,tpi,qji+i*nconv);
592:     }
593: #endif
594:     for (i=derivative?2:1;i<pep->nmat;i++) BVMultVec(pjd->AX[i],1.0,1.0,tp,qj+i*nconv);

596:     /* extended vector part */
597:     BVSetActiveColumns(pjd->X,0,nconv);
598:     BVDotVec(pjd->X,tu,xx);
599:     xxi = xx+nconv;
600: #if !defined(PETSC_USE_COMPLEX)
601:     if (sz==2) BVDotVec(pjd->X,tui,xxi);
602: #endif
603:     if (sz==1) PetscArrayzero(xxi,nconv);
604:     if (rk==np-1) {
605:       PetscBLASIntCast(nconv,&n);
606:       PetscBLASIntCast(pjd->ld,&ld);
607:       y2  = array2+nloc;
608:       PetscArrayzero(y2,nconv);
609:       for (j=derivative?1:0;j<pjd->midx;j++) {
610:         for (i=0;i<nconv;i++) tt[i] = dval[j]*xx[i]-dvali[j]*xxi[i];
611:         PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
612:         PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
613:       }
614:       for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
615: #if !defined(PETSC_USE_COMPLEX)
616:       if (sz==2) {
617:         y2i = arrayi2+nloc;
618:         PetscArrayzero(y2i,nconv);
619:         for (j=derivative?1:0;j<pjd->midx;j++) {
620:           for (i=0;i<nconv;i++) tt[i] = dval[j]*xxi[i]+dvali[j]*xx[i];
621:           PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
622:           PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
623:         }
624:         for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
625:       }
626: #endif
627:     }
628:     PetscMPIIntCast(nconv,&count);
629:     MPI_Bcast(array2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
630: #if !defined(PETSC_USE_COMPLEX)
631:     if (sz==2) MPI_Bcast(arrayi2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
632: #endif
633:   }
634:   if (nconv) PetscFree5(dval,xx,tt,x2,qj);
635:   else PetscFree(dval);
636:   VecResetArray(tu);
637:   VecRestoreArray(u[0],&array1);
638:   VecResetArray(tp);
639:   VecRestoreArray(p[0],&array2);
640: #if !defined(PETSC_USE_COMPLEX)
641:   if (sz==2) {
642:     VecResetArray(tui);
643:     VecRestoreArray(u[1],&arrayi1);
644:     VecResetArray(tpi);
645:     VecRestoreArray(p[1],&arrayi2);
646:   }
647: #endif
648:   return 0;
649: }

651: static PetscErrorCode PEPJDProcessInitialSpace(PEP pep,Vec *w)
652: {
653:   PEP_JD         *pjd = (PEP_JD*)pep->data;
654:   PetscScalar    *tt,target[2];
655:   Vec            vg,wg;
656:   PetscInt       i;
657:   PetscReal      norm;

659:   PetscMalloc1(pjd->ld-1,&tt);
661:   BVSetRandomColumn(pjd->V,0);
662:   for (i=0;i<pjd->ld-1;i++) tt[i] = 0.0;
663:   BVGetColumn(pjd->V,0,&vg);
664:   PEPJDCopyToExtendedVec(pep,NULL,tt,pjd->ld-1,0,vg,PETSC_FALSE);
665:   BVRestoreColumn(pjd->V,0,&vg);
666:   BVNormColumn(pjd->V,0,NORM_2,&norm);
667:   BVScaleColumn(pjd->V,0,1.0/norm);
668:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
669:     BVGetColumn(pjd->V,0,&vg);
670:     BVGetColumn(pjd->W,0,&wg);
671:     VecSet(wg,0.0);
672:     target[0] = pep->target; target[1] = 0.0;
673:     PEPJDComputeResidual(pep,PETSC_TRUE,1,&vg,target,&wg,w);
674:     BVRestoreColumn(pjd->W,0,&wg);
675:     BVRestoreColumn(pjd->V,0,&vg);
676:     BVNormColumn(pjd->W,0,NORM_2,&norm);
677:     BVScaleColumn(pjd->W,0,1.0/norm);
678:   }
679:   PetscFree(tt);
680:   return 0;
681: }

683: static PetscErrorCode MatMult_PEPJD(Mat P,Vec x,Vec y)
684: {
685:   PEP_JD_MATSHELL *matctx;
686:   PEP_JD          *pjd;
687:   PetscInt        i,j,nconv,nloc,nmat,ldt,ncv,sz;
688:   Vec             tx,ty;
689:   const Vec       *xs,*ys;
690:   PetscScalar     *array1,*array2,*x2=NULL,*y2,*tt=NULL,*xx=NULL,*xxi,theta[2],sone=1.0,*qj,*val,*vali=NULL;
691:   PetscBLASInt    n,ld,one=1;
692:   PetscMPIInt     np;
693: #if !defined(PETSC_USE_COMPLEX)
694:   Vec             txi=NULL,tyi=NULL;
695:   PetscScalar     *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
696: #endif

698:   MPI_Comm_size(PetscObjectComm((PetscObject)P),&np);
699:   MatShellGetContext(P,&matctx);
700:   pjd   = (PEP_JD*)(matctx->pep->data);
701:   nconv = pjd->nlock;
702:   nmat  = matctx->pep->nmat;
703:   ncv   = matctx->pep->ncv;
704:   ldt   = pjd->ld;
705:   VecCompGetSubVecs(x,&sz,&xs);
706:   VecCompGetSubVecs(y,NULL,&ys);
707:   theta[0] = matctx->theta[0];
708:   theta[1] = (sz==2)?matctx->theta[1]:0.0;
709:   if (nconv>0) {
710:     PetscMalloc5(nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*nmat,&qj,2*nconv,&xx,2*nmat,&val);
711:     BVGetSizes(matctx->pep->V,&nloc,NULL,NULL);
712:     VecGetArray(xs[0],&array1);
713:     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
714:     VecRestoreArray(xs[0],&array1);
715: #if !defined(PETSC_USE_COMPLEX)
716:     if (sz==2) {
717:       x2i = x2+nconv;
718:       VecGetArray(xs[1],&arrayi1);
719:       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]* PetscSqrtReal(np);
720:       VecRestoreArray(xs[1],&arrayi1);
721:     }
722: #endif
723:     vali = val+nmat;
724:   }
725:   tx = matctx->work[0];
726:   ty = matctx->work[1];
727:   VecGetArray(xs[0],&array1);
728:   VecPlaceArray(tx,array1);
729:   VecGetArray(ys[0],&array2);
730:   VecPlaceArray(ty,array2);
731:   MatMult(matctx->Pr,tx,ty);
732: #if !defined(PETSC_USE_COMPLEX)
733:   if (sz==2) {
734:     txi = matctx->work[2];
735:     tyi = matctx->work[3];
736:     VecGetArray(xs[1],&arrayi1);
737:     VecPlaceArray(txi,arrayi1);
738:     VecGetArray(ys[1],&arrayi2);
739:     VecPlaceArray(tyi,arrayi2);
740:     MatMult(matctx->Pr,txi,tyi);
741:     if (theta[1]!=0.0) {
742:       MatMult(matctx->Pi,txi,matctx->work[4]);
743:       VecAXPY(ty,-1.0,matctx->work[4]);
744:       MatMult(matctx->Pi,tx,matctx->work[4]);
745:       VecAXPY(tyi,1.0,matctx->work[4]);
746:     }
747:   }
748: #endif
749:   if (nconv>0) {
750:     PEPEvaluateBasis(matctx->pep,theta[0],theta[1],val,vali);
751:     for (i=0;i<nmat;i++) 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);
752: #if !defined(PETSC_USE_COMPLEX)
753:     if (sz==2) {
754:       qji = qj+nconv*nmat;
755:       qq = qji+nconv*nmat;
756:       for (i=0;i<nmat;i++) 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);
757:       for (i=0;i<nconv*nmat;i++) qj[i] -= qji[i];
758:       for (i=0;i<nmat;i++) {
759:         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);
760:         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);
761:       }
762:       for (i=0;i<nconv*nmat;i++) qji[i] += qq[i];
763:       for (i=1;i<matctx->pep->nmat;i++) BVMultVec(pjd->AX[i],1.0,1.0,tyi,qji+i*nconv);
764:     }
765: #endif
766:     for (i=1;i<nmat;i++) BVMultVec(pjd->AX[i],1.0,1.0,ty,qj+i*nconv);

768:     /* extended vector part */
769:     BVSetActiveColumns(pjd->X,0,nconv);
770:     BVDotVec(pjd->X,tx,xx);
771:     xxi = xx+nconv;
772: #if !defined(PETSC_USE_COMPLEX)
773:     if (sz==2) BVDotVec(pjd->X,txi,xxi);
774: #endif
775:     if (sz==1) PetscArrayzero(xxi,nconv);
776:     PetscBLASIntCast(pjd->nlock,&n);
777:     PetscBLASIntCast(ldt,&ld);
778:     y2 = array2+nloc;
779:     PetscArrayzero(y2,nconv);
780:     for (j=0;j<pjd->midx;j++) {
781:       for (i=0;i<nconv;i++) tt[i] = val[j]*xx[i]-vali[j]*xxi[i];
782:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
783:       PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
784:     }
785: #if !defined(PETSC_USE_COMPLEX)
786:     if (sz==2) {
787:       y2i = arrayi2+nloc;
788:       PetscArrayzero(y2i,nconv);
789:       for (j=0;j<pjd->midx;j++) {
790:         for (i=0;i<nconv;i++) tt[i] = val[j]*xxi[i]+vali[j]*xx[i];
791:         PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
792:         PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
793:       }
794:       for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
795:     }
796: #endif
797:     for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
798:     PetscFree5(tt,x2,qj,xx,val);
799:   }
800:   VecResetArray(tx);
801:   VecRestoreArray(xs[0],&array1);
802:   VecResetArray(ty);
803:   VecRestoreArray(ys[0],&array2);
804: #if !defined(PETSC_USE_COMPLEX)
805:   if (sz==2) {
806:     VecResetArray(txi);
807:     VecRestoreArray(xs[1],&arrayi1);
808:     VecResetArray(tyi);
809:     VecRestoreArray(ys[1],&arrayi2);
810:   }
811: #endif
812:   return 0;
813: }

815: static PetscErrorCode MatCreateVecs_PEPJD(Mat A,Vec *right,Vec *left)
816: {
817:   PEP_JD_MATSHELL *matctx;
818:   PEP_JD          *pjd;
819:   PetscInt        kspsf=1,i;
820:   Vec             v[2];

822:   MatShellGetContext(A,&matctx);
823:   pjd   = (PEP_JD*)(matctx->pep->data);
824: #if !defined (PETSC_USE_COMPLEX)
825:   kspsf = 2;
826: #endif
827:   for (i=0;i<kspsf;i++) BVCreateVec(pjd->V,v+i);
828:   if (right) VecCreateCompWithVecs(v,kspsf,pjd->vtempl,right);
829:   if (left) VecCreateCompWithVecs(v,kspsf,pjd->vtempl,left);
830:   for (i=0;i<kspsf;i++) VecDestroy(&v[i]);
831:   return 0;
832: }

834: static PetscErrorCode PEPJDUpdateExtendedPC(PEP pep,PetscScalar theta)
835: {
836:   PEP_JD         *pjd = (PEP_JD*)pep->data;
837:   PEP_JD_PCSHELL *pcctx;
838:   PetscInt       i,j,k,n=pjd->nlock,ld=pjd->ld,deg=pep->nmat-1;
839:   PetscScalar    *M,*ps,*work,*U,*V,*S,*Sp,*Spp,snone=-1.0,sone=1.0,zero=0.0,*val;
840:   PetscReal      tol,maxeig=0.0,*sg,*rwork;
841:   PetscBLASInt   n_,info,ld_,*p,lw_,rk=0;

843:   if (n) {
844:     PCShellGetContext(pjd->pcshell,&pcctx);
845:     pcctx->theta = theta;
846:     pcctx->n = n;
847:     M  = pcctx->M;
848:     PetscBLASIntCast(n,&n_);
849:     PetscBLASIntCast(ld,&ld_);
850:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
851:     if (pjd->midx==1) {
852:       PetscArraycpy(M,pjd->XpX,ld*ld);
853:       PetscCalloc2(10*n,&work,n,&p);
854:     } else {
855:       ps = pcctx->ps;
856:       PetscCalloc7(2*n*n,&U,3*n*n,&S,n,&sg,10*n,&work,5*n,&rwork,n,&p,deg+1,&val);
857:       V = U+n*n;
858:       /* pseudo-inverse */
859:       for (j=0;j<n;j++) {
860:         for (i=0;i<n;i++) S[n*j+i] = -pjd->T[pep->ncv*j+i];
861:         S[n*j+j] += theta;
862:       }
863:       lw_ = 10*n_;
864: #if !defined (PETSC_USE_COMPLEX)
865:       PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,&info));
866: #else
867:       PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,rwork,&info));
868: #endif
869:       SlepcCheckLapackInfo("gesvd",info);
870:       for (i=0;i<n;i++) maxeig = PetscMax(maxeig,sg[i]);
871:       tol = 10*PETSC_MACHINE_EPSILON*n*maxeig;
872:       for (j=0;j<n;j++) {
873:         if (sg[j]>tol) {
874:           for (i=0;i<n;i++) U[j*n+i] /= sg[j];
875:           rk++;
876:         } else break;
877:       }
878:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&rk,&sone,U,&n_,V,&n_,&zero,ps,&ld_));

880:       /* compute M */
881:       PEPEvaluateBasis(pep,theta,0.0,val,NULL);
882:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&snone,pjd->XpX,&ld_,ps,&ld_,&zero,M,&ld_));
883:       PetscArrayzero(S,2*n*n);
884:       Sp = S+n*n;
885:       for (j=0;j<n;j++) S[j*(n+1)] = 1.0;
886:       for (k=1;k<pjd->midx;k++) {
887:         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];
888:         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pjd->XpX,&ld_,V,&n_,&zero,U,&n_));
889:         PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+k*ld*ld,&ld_,U,&n_,&sone,M,&ld_));
890:         Spp = Sp; Sp = S;
891:         PEPJDEvaluateHatBasis(pep,n,pjd->T,pep->ncv,val,NULL,k+1,Spp,Sp,S);
892:       }
893:     }
894:     /* inverse */
895:     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,M,&ld_,p,&info));
896:     SlepcCheckLapackInfo("getrf",info);
897:     PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&n_,M,&ld_,p,work,&n_,&info));
898:     SlepcCheckLapackInfo("getri",info);
899:     PetscFPTrapPop();
900:     if (pjd->midx==1) PetscFree2(work,p);
901:     else PetscFree7(U,S,sg,work,rwork,p,val);
902:   }
903:   return 0;
904: }

906: static PetscErrorCode PEPJDMatSetUp(PEP pep,PetscInt sz,PetscScalar *theta)
907: {
908:   PEP_JD          *pjd = (PEP_JD*)pep->data;
909:   PEP_JD_MATSHELL *matctx;
910:   PEP_JD_PCSHELL  *pcctx;
911:   MatStructure    str;
912:   PetscScalar     *vals,*valsi;
913:   PetscBool       skipmat=PETSC_FALSE;
914:   PetscInt        i;
915:   Mat             Pr=NULL;

917:   if (sz==2 && theta[1]==0.0) sz = 1;
918:   MatShellGetContext(pjd->Pshell,&matctx);
919:   PCShellGetContext(pjd->pcshell,&pcctx);
920:   if (matctx->Pr && matctx->theta[0]==theta[0] && ((!matctx->Pi && sz==1) || (sz==2 && matctx->theta[1]==theta[1]))) {
921:     if (pcctx->n == pjd->nlock) return 0;
922:     skipmat = PETSC_TRUE;
923:   }
924:   if (!skipmat) {
925:     PetscMalloc2(pep->nmat,&vals,pep->nmat,&valsi);
926:     STGetMatStructure(pep->st,&str);
927:     PEPEvaluateBasis(pep,theta[0],theta[1],vals,valsi);
928:     if (!matctx->Pr) MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pr);
929:     else MatCopy(pep->A[0],matctx->Pr,str);
930:     for (i=1;i<pep->nmat;i++) MatAXPY(matctx->Pr,vals[i],pep->A[i],str);
931:     if (!pjd->reusepc) {
932:       if (pcctx->PPr && sz==2) {
933:         MatCopy(matctx->Pr,pcctx->PPr,str);
934:         Pr = pcctx->PPr;
935:       } else Pr = matctx->Pr;
936:     }
937:     matctx->theta[0] = theta[0];
938: #if !defined(PETSC_USE_COMPLEX)
939:     if (sz==2) {
940:       if (!matctx->Pi) MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pi);
941:       else MatCopy(pep->A[1],matctx->Pi,str);
942:       MatScale(matctx->Pi,valsi[1]);
943:       for (i=2;i<pep->nmat;i++) MatAXPY(matctx->Pi,valsi[i],pep->A[i],str);
944:       matctx->theta[1] = theta[1];
945:     }
946: #endif
947:     PetscFree2(vals,valsi);
948:   }
949:   if (!pjd->reusepc) {
950:     if (!skipmat) {
951:       PCSetOperators(pcctx->pc,Pr,Pr);
952:       PCSetUp(pcctx->pc);
953:     }
954:     PEPJDUpdateExtendedPC(pep,theta[0]);
955:   }
956:   return 0;
957: }

959: static PetscErrorCode PEPJDCreateShellPC(PEP pep,Vec *ww)
960: {
961:   PEP_JD          *pjd = (PEP_JD*)pep->data;
962:   PEP_JD_PCSHELL  *pcctx;
963:   PEP_JD_MATSHELL *matctx;
964:   KSP             ksp;
965:   PetscInt        nloc,mloc,kspsf=1;
966:   Vec             v[2];
967:   PetscScalar     target[2];
968:   Mat             Pr;

970:   /* Create the reference vector */
971:   BVGetColumn(pjd->V,0,&v[0]);
972:   v[1] = v[0];
973: #if !defined (PETSC_USE_COMPLEX)
974:   kspsf = 2;
975: #endif
976:   VecCreateCompWithVecs(v,kspsf,NULL,&pjd->vtempl);
977:   BVRestoreColumn(pjd->V,0,&v[0]);

979:   /* Replace preconditioner with one containing projectors */
980:   PCCreate(PetscObjectComm((PetscObject)pep),&pjd->pcshell);
981:   PCSetType(pjd->pcshell,PCSHELL);
982:   PCShellSetName(pjd->pcshell,"PCPEPJD");
983:   PCShellSetApply(pjd->pcshell,PCShellApply_PEPJD);
984:   PetscNew(&pcctx);
985:   PCShellSetContext(pjd->pcshell,pcctx);
986:   STGetKSP(pep->st,&ksp);
987:   BVCreateVec(pjd->V,&pcctx->Bp[0]);
988:   VecDuplicate(pcctx->Bp[0],&pcctx->Bp[1]);
989:   KSPGetPC(ksp,&pcctx->pc);
990:   PetscObjectReference((PetscObject)pcctx->pc);
991:   MatGetLocalSize(pep->A[0],&mloc,&nloc);
992:   if (pjd->ld>1) {
993:     nloc += pjd->ld-1; mloc += pjd->ld-1;
994:   }
995:   PetscNew(&matctx);
996:   MatCreateShell(PetscObjectComm((PetscObject)pep),kspsf*nloc,kspsf*mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&pjd->Pshell);
997:   MatShellSetOperation(pjd->Pshell,MATOP_MULT,(void(*)(void))MatMult_PEPJD);
998:   MatShellSetOperation(pjd->Pshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_PEPJD);
999:   matctx->pep = pep;
1000:   target[0] = pep->target; target[1] = 0.0;
1001:   PEPJDMatSetUp(pep,1,target);
1002:   Pr = matctx->Pr;
1003:   pcctx->PPr = NULL;
1004: #if !defined(PETSC_USE_COMPLEX)
1005:   if (!pjd->reusepc) {
1006:     MatDuplicate(matctx->Pr,MAT_COPY_VALUES,&pcctx->PPr);
1007:     Pr = pcctx->PPr;
1008:   }
1009: #endif
1010:   PCSetOperators(pcctx->pc,Pr,Pr);
1011:   PCSetErrorIfFailure(pcctx->pc,PETSC_TRUE);
1012:   KSPSetPC(ksp,pjd->pcshell);
1013:   if (pjd->reusepc) {
1014:     PCSetReusePreconditioner(pcctx->pc,PETSC_TRUE);
1015:     KSPSetReusePreconditioner(ksp,PETSC_TRUE);
1016:   }
1017:   PEP_KSPSetOperators(ksp,pjd->Pshell,pjd->Pshell);
1018:   KSPSetUp(ksp);
1019:   if (pjd->ld>1) {
1020:     PetscMalloc2(pjd->ld*pjd->ld,&pcctx->M,pjd->ld*pjd->ld,&pcctx->ps);
1021:     pcctx->pep = pep;
1022:   }
1023:   matctx->work = ww;
1024:   pcctx->work  = ww;
1025:   return 0;
1026: }

1028: static PetscErrorCode PEPJDEigenvectors(PEP pep)
1029: {
1030:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1031:   PetscBLASInt   ld,nconv,info,nc;
1032:   PetscScalar    *Z;
1033:   PetscReal      *wr;
1034:   Mat            U;
1035: #if defined(PETSC_USE_COMPLEX)
1036:   PetscScalar    *w;
1037: #endif

1039:   PetscBLASIntCast(pep->ncv,&ld);
1040:   PetscBLASIntCast(pep->nconv,&nconv);
1041:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1042: #if !defined(PETSC_USE_COMPLEX)
1043:   PetscMalloc2(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr);
1044:   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,wr,&info));
1045: #else
1046:   PetscMalloc3(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr,2*pep->ncv,&w);
1047:   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,w,wr,&info));
1048: #endif
1049:   PetscFPTrapPop();
1050:   SlepcCheckLapackInfo("trevc",info);
1051:   MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,Z,&U);
1052:   BVSetActiveColumns(pjd->X,0,pep->nconv);
1053:   BVMultInPlace(pjd->X,U,0,pep->nconv);
1054:   BVNormalize(pjd->X,pep->eigi);
1055:   MatDestroy(&U);
1056: #if !defined(PETSC_USE_COMPLEX)
1057:   PetscFree2(Z,wr);
1058: #else
1059:   PetscFree3(Z,wr,w);
1060: #endif
1061:   return 0;
1062: }

1064: static PetscErrorCode PEPJDLockConverged(PEP pep,PetscInt *nv,PetscInt sz)
1065: {
1066:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1067:   PetscInt       j,i,*P,ldds,rk=0,nvv=*nv;
1068:   Vec            v,x,w;
1069:   PetscScalar    *R,*r,*pX,target[2];
1070:   Mat            X;
1071:   PetscBLASInt   sz_,rk_,nv_,info;
1072:   PetscMPIInt    np;

1074:   /* update AX and XpX */
1075:   for (i=sz;i>0;i--) {
1076:     BVGetColumn(pjd->X,pjd->nlock-i,&x);
1077:     for (j=0;j<pep->nmat;j++) {
1078:       BVGetColumn(pjd->AX[j],pjd->nlock-i,&v);
1079:       MatMult(pep->A[j],x,v);
1080:       BVRestoreColumn(pjd->AX[j],pjd->nlock-i,&v);
1081:       BVSetActiveColumns(pjd->AX[j],0,pjd->nlock-i+1);
1082:     }
1083:     BVRestoreColumn(pjd->X,pjd->nlock-i,&x);
1084:     BVDotColumn(pjd->X,(pjd->nlock-i),pjd->XpX+(pjd->nlock-i)*(pjd->ld));
1085:     pjd->XpX[(pjd->nlock-i)*(1+pjd->ld)] = 1.0;
1086:     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]);
1087:   }

1089:   /* minimality index */
1090:   pjd->midx = PetscMin(pjd->mmidx,pjd->nlock);

1092:   /* evaluate the polynomial basis in T */
1093:   PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat);
1094:   for (j=0;j<pep->nmat;j++) 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);

1096:   /* Extend search space */
1097:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
1098:   PetscCalloc3(nvv,&P,nvv*nvv,&R,nvv*sz,&r);
1099:   DSGetLeadingDimension(pep->ds,&ldds);
1100:   DSGetArray(pep->ds,DS_MAT_X,&pX);
1101:   PEPJDOrthogonalize(nvv,nvv,pX,ldds,&rk,P,R,nvv);
1102:   for (j=0;j<sz;j++) {
1103:     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 */
1104:   }
1105:   PetscBLASIntCast(rk,&rk_);
1106:   PetscBLASIntCast(sz,&sz_);
1107:   PetscBLASIntCast(nvv,&nv_);
1108:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1109:   PetscCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&rk_,R,&nv_,&info));
1110:   PetscFPTrapPop();
1111:   SlepcCheckLapackInfo("trtri",info);
1112:   for (i=0;i<sz;i++) PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&rk_,R,&nv_,r+i,&sz_));
1113:   for (i=0;i<sz*rk;i++) r[i] = PetscConj(r[i])/PetscSqrtReal(np); /* revert */
1114:   BVSetActiveColumns(pjd->V,0,nvv);
1115:   rk -= sz;
1116:   for (j=0;j<rk;j++) PetscArraycpy(R+j*nvv,pX+(j+sz)*ldds,nvv);
1117:   DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1118:   MatCreateSeqDense(PETSC_COMM_SELF,nvv,rk,R,&X);
1119:   BVMultInPlace(pjd->V,X,0,rk);
1120:   MatDestroy(&X);
1121:   BVSetActiveColumns(pjd->V,0,rk);
1122:   for (j=0;j<rk;j++) {
1123:     BVGetColumn(pjd->V,j,&v);
1124:     PEPJDCopyToExtendedVec(pep,NULL,r+sz*(j+sz),sz,pjd->nlock-sz,v,PETSC_FALSE);
1125:     BVRestoreColumn(pjd->V,j,&v);
1126:   }
1127:   BVOrthogonalize(pjd->V,NULL);

1129:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1130:     for (j=0;j<rk;j++) {
1131:       /* W = P(target)*V */
1132:       BVGetColumn(pjd->W,j,&w);
1133:       BVGetColumn(pjd->V,j,&v);
1134:       target[0] = pep->target; target[1] = 0.0;
1135:       PEPJDComputeResidual(pep,PETSC_FALSE,1,&v,target,&w,pep->work);
1136:       BVRestoreColumn(pjd->V,j,&v);
1137:       BVRestoreColumn(pjd->W,j,&w);
1138:     }
1139:     BVSetActiveColumns(pjd->W,0,rk);
1140:     BVOrthogonalize(pjd->W,NULL);
1141:   }
1142:   *nv = rk;
1143:   PetscFree3(P,R,r);
1144:   return 0;
1145: }

1147: PetscErrorCode PEPJDSystemSetUp(PEP pep,PetscInt sz,PetscScalar *theta,Vec *u,Vec *p,Vec *ww)
1148: {
1149:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1150:   PEP_JD_PCSHELL *pcctx;
1151: #if !defined(PETSC_USE_COMPLEX)
1152:   PetscScalar    s[2];
1153: #endif

1155:   PCShellGetContext(pjd->pcshell,&pcctx);
1156:   PEPJDMatSetUp(pep,sz,theta);
1157:   pcctx->u[0] = u[0]; pcctx->u[1] = u[1];
1158:   /* Compute r'. p is a work space vector */
1159:   PEPJDComputeResidual(pep,PETSC_TRUE,sz,u,theta,p,ww);
1160:   PEPJDExtendedPCApply(pjd->pcshell,p[0],pcctx->Bp[0]);
1161:   VecDot(pcctx->Bp[0],u[0],pcctx->gamma);
1162: #if !defined(PETSC_USE_COMPLEX)
1163:   if (sz==2) {
1164:     PEPJDExtendedPCApply(pjd->pcshell,p[1],pcctx->Bp[1]);
1165:     VecDot(pcctx->Bp[0],u[1],pcctx->gamma+1);
1166:     VecMDot(pcctx->Bp[1],2,u,s);
1167:     pcctx->gamma[0] += s[1];
1168:     pcctx->gamma[1] = -pcctx->gamma[1]+s[0];
1169:   }
1170: #endif
1171:   if (sz==1) {
1172:     VecZeroEntries(pcctx->Bp[1]);
1173:     pcctx->gamma[1] = 0.0;
1174:   }
1175:   return 0;
1176: }

1178: PetscErrorCode PEPSolve_JD(PEP pep)
1179: {
1180:   PEP_JD          *pjd = (PEP_JD*)pep->data;
1181:   PetscInt        k,nv,nvc,ld,minv,dim,bupdated=0,sz=1,kspsf=1,idx,off,maxits,nloc;
1182:   PetscMPIInt     np,count;
1183:   PetscScalar     theta[2]={0.0,0.0},ritz[2]={0.0,0.0},*pX,*eig,*eigi,*array;
1184:   PetscReal       norm,*res,tol=0.0,rtol,abstol, dtol;
1185:   PetscBool       lindep,ini=PETSC_TRUE;
1186:   Vec             tc,t[2]={NULL,NULL},u[2]={NULL,NULL},p[2]={NULL,NULL};
1187:   Vec             rc,rr[2],r[2]={NULL,NULL},*ww=pep->work,v[2];
1188:   Mat             G,X,Y;
1189:   KSP             ksp;
1190:   PEP_JD_PCSHELL  *pcctx;
1191:   PEP_JD_MATSHELL *matctx;
1192: #if !defined(PETSC_USE_COMPLEX)
1193:   PetscReal       norm1;
1194: #endif

1196:   PetscCitationsRegister(citation,&cited);
1197:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
1198:   BVGetSizes(pep->V,&nloc,NULL,NULL);
1199:   DSGetLeadingDimension(pep->ds,&ld);
1200:   PetscCalloc3(pep->ncv+pep->nev,&eig,pep->ncv+pep->nev,&eigi,pep->ncv+pep->nev,&res);
1201:   pjd->nlock = 0;
1202:   STGetKSP(pep->st,&ksp);
1203:   KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits);
1204: #if !defined (PETSC_USE_COMPLEX)
1205:   kspsf = 2;
1206: #endif
1207:   PEPJDProcessInitialSpace(pep,ww);
1208:   nv = (pep->nini)?pep->nini:1;

1210:   /* Replace preconditioner with one containing projectors */
1211:   PEPJDCreateShellPC(pep,ww);
1212:   PCShellGetContext(pjd->pcshell,&pcctx);

1214:   /* Create auxiliary vectors */
1215:   BVCreateVec(pjd->V,&u[0]);
1216:   VecDuplicate(u[0],&p[0]);
1217:   VecDuplicate(u[0],&r[0]);
1218: #if !defined (PETSC_USE_COMPLEX)
1219:   VecDuplicate(u[0],&u[1]);
1220:   VecDuplicate(u[0],&p[1]);
1221:   VecDuplicate(u[0],&r[1]);
1222: #endif

1224:   /* Restart loop */
1225:   while (pep->reason == PEP_CONVERGED_ITERATING) {
1226:     pep->its++;
1227:     DSSetDimensions(pep->ds,nv,0,0);
1228:     BVSetActiveColumns(pjd->V,bupdated,nv);
1229:     PEPJDUpdateTV(pep,bupdated,nv,ww);
1230:     if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) BVSetActiveColumns(pjd->W,bupdated,nv);
1231:     for (k=0;k<pep->nmat;k++) {
1232:       BVSetActiveColumns(pjd->TV[k],bupdated,nv);
1233:       DSGetMat(pep->ds,DSMatExtra[k],&G);
1234:       BVMatProject(pjd->TV[k],NULL,pjd->W,G);
1235:       DSRestoreMat(pep->ds,DSMatExtra[k],&G);
1236:     }
1237:     BVSetActiveColumns(pjd->V,0,nv);
1238:     BVSetActiveColumns(pjd->W,0,nv);

1240:     /* Solve projected problem */
1241:     DSSetState(pep->ds,DS_STATE_RAW);
1242:     DSSolve(pep->ds,pep->eigr,pep->eigi);
1243:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
1244:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);
1245:     idx = 0;
1246:     do {
1247:       ritz[0] = pep->eigr[idx];
1248: #if !defined(PETSC_USE_COMPLEX)
1249:       ritz[1] = pep->eigi[idx];
1250:       sz = (ritz[1]==0.0)?1:2;
1251: #endif
1252:       /* Compute Ritz vector u=V*X(:,1) */
1253:       DSGetArray(pep->ds,DS_MAT_X,&pX);
1254:       BVSetActiveColumns(pjd->V,0,nv);
1255:       BVMultVec(pjd->V,1.0,0.0,u[0],pX+idx*ld);
1256: #if !defined(PETSC_USE_COMPLEX)
1257:       if (sz==2) BVMultVec(pjd->V,1.0,0.0,u[1],pX+(idx+1)*ld);
1258: #endif
1259:       DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1260:       PEPJDComputeResidual(pep,PETSC_FALSE,sz,u,ritz,r,ww);
1261:       /* Check convergence */
1262:       VecNorm(r[0],NORM_2,&norm);
1263: #if !defined(PETSC_USE_COMPLEX)
1264:       if (sz==2) {
1265:         VecNorm(r[1],NORM_2,&norm1);
1266:         norm = SlepcAbs(norm,norm1);
1267:       }
1268: #endif
1269:       (*pep->converged)(pep,ritz[0],ritz[1],norm,&pep->errest[pep->nconv],pep->convergedctx);
1270:       if (sz==2) pep->errest[pep->nconv+1] = pep->errest[pep->nconv];
1271:       if (ini) {
1272:         tol = PetscMin(.1,pep->errest[pep->nconv]); ini = PETSC_FALSE;
1273:       } else tol = PetscMin(pep->errest[pep->nconv],tol/2);
1274:       (*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);
1275:       if (pep->errest[pep->nconv]<pep->tol) {
1276:         /* Ritz pair converged */
1277:         ini = PETSC_TRUE;
1278:         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1279:         if (pjd->ld>1) {
1280:           BVGetColumn(pjd->X,pep->nconv,&v[0]);
1281:           PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*pep->nconv,pjd->ld-1,0,u[0],PETSC_TRUE);
1282:           BVRestoreColumn(pjd->X,pep->nconv,&v[0]);
1283:           BVSetActiveColumns(pjd->X,0,pep->nconv+1);
1284:           BVNormColumn(pjd->X,pep->nconv,NORM_2,&norm);
1285:           BVScaleColumn(pjd->X,pep->nconv,1.0/norm);
1286:           for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*pep->nconv+k] *= PetscSqrtReal(np)/norm;
1287:           pjd->T[(pep->ncv+1)*pep->nconv] = ritz[0];
1288:           eig[pep->nconv] = ritz[0];
1289:           idx++;
1290: #if !defined(PETSC_USE_COMPLEX)
1291:           if (sz==2) {
1292:             BVGetColumn(pjd->X,pep->nconv+1,&v[0]);
1293:             PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*(pep->nconv+1),pjd->ld-1,0,u[1],PETSC_TRUE);
1294:             BVRestoreColumn(pjd->X,pep->nconv+1,&v[0]);
1295:             BVSetActiveColumns(pjd->X,0,pep->nconv+2);
1296:             BVNormColumn(pjd->X,pep->nconv+1,NORM_2,&norm1);
1297:             BVScaleColumn(pjd->X,pep->nconv+1,1.0/norm1);
1298:             for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*(pep->nconv+1)+k] *= PetscSqrtReal(np)/norm1;
1299:             pjd->T[(pep->ncv+1)*(pep->nconv+1)] = ritz[0];
1300:             pjd->T[(pep->ncv+1)*pep->nconv+1] = -ritz[1]*norm1/norm;
1301:             pjd->T[(pep->ncv+1)*(pep->nconv+1)-1] = ritz[1]*norm/norm1;
1302:             eig[pep->nconv+1] = ritz[0];
1303:             eigi[pep->nconv] = ritz[1]; eigi[pep->nconv+1] = -ritz[1];
1304:             idx++;
1305:           }
1306: #endif
1307:         } else BVInsertVec(pep->V,pep->nconv,u[0]);
1308:         pep->nconv += sz;
1309:       }
1310:     } while (pep->errest[pep->nconv]<pep->tol && pep->nconv<nv);

1312:     if (pep->reason==PEP_CONVERGED_ITERATING) {
1313:       nvc = nv;
1314:       if (idx) {
1315:         pjd->nlock +=idx;
1316:         PEPJDLockConverged(pep,&nv,idx);
1317:       }
1318:       if (nv+sz>=pep->ncv-1) {
1319:         /* Basis full, force restart */
1320:         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1321:         DSGetDimensions(pep->ds,&dim,NULL,NULL,NULL);
1322:         DSGetArray(pep->ds,DS_MAT_X,&pX);
1323:         PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld);
1324:         DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1325:         DSGetArray(pep->ds,DS_MAT_Y,&pX);
1326:         PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld);
1327:         DSRestoreArray(pep->ds,DS_MAT_Y,&pX);
1328:         DSGetMat(pep->ds,DS_MAT_X,&X);
1329:         BVMultInPlace(pjd->V,X,0,minv);
1330:         DSRestoreMat(pep->ds,DS_MAT_X,&X);
1331:         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1332:          DSGetMat(pep->ds,DS_MAT_Y,&Y);
1333:          BVMultInPlace(pjd->W,Y,pep->nconv,minv);
1334:          DSRestoreMat(pep->ds,DS_MAT_Y,&Y);
1335:         }
1336:         nv = minv;
1337:         bupdated = 0;
1338:       } else {
1339:         if (!idx && pep->errest[pep->nconv]<pjd->fix) {theta[0] = ritz[0]; theta[1] = ritz[1];}
1340:         else {theta[0] = pep->target; theta[1] = 0.0;}
1341:         /* Update system mat */
1342:         PEPJDSystemSetUp(pep,sz,theta,u,p,ww);
1343:         /* Solve correction equation to expand basis */
1344:         BVGetColumn(pjd->V,nv,&t[0]);
1345:         rr[0] = r[0];
1346:         if (sz==2) {
1347:           BVGetColumn(pjd->V,nv+1,&t[1]);
1348:           rr[1] = r[1];
1349:         } else {
1350:           t[1] = NULL;
1351:           rr[1] = NULL;
1352:         }
1353:         VecCreateCompWithVecs(t,kspsf,pjd->vtempl,&tc);
1354:         VecCreateCompWithVecs(rr,kspsf,pjd->vtempl,&rc);
1355:         VecCompSetSubVecs(pjd->vtempl,sz,NULL);
1356:         tol  = PetscMax(rtol,tol/2);
1357:         KSPSetTolerances(ksp,tol,abstol,dtol,maxits);
1358:         KSPSolve(ksp,rc,tc);
1359:         VecDestroy(&tc);
1360:         VecDestroy(&rc);
1361:         VecGetArray(t[0],&array);
1362:         PetscMPIIntCast(pep->nconv,&count);
1363:         MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
1364:         VecRestoreArray(t[0],&array);
1365:         BVRestoreColumn(pjd->V,nv,&t[0]);
1366:         BVOrthogonalizeColumn(pjd->V,nv,NULL,&norm,&lindep);
1367:         if (lindep || norm==0.0) {
1369:           off = 1;
1370:         } else {
1371:           off = 0;
1372:           BVScaleColumn(pjd->V,nv,1.0/norm);
1373:         }
1374: #if !defined(PETSC_USE_COMPLEX)
1375:         if (sz==2) {
1376:           VecGetArray(t[1],&array);
1377:           MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
1378:           VecRestoreArray(t[1],&array);
1379:           BVRestoreColumn(pjd->V,nv+1,&t[1]);
1380:           if (off) BVCopyColumn(pjd->V,nv+1,nv);
1381:           BVOrthogonalizeColumn(pjd->V,nv+1-off,NULL,&norm,&lindep);
1382:           if (lindep || norm==0.0) {
1384:             off = 1;
1385:           } else BVScaleColumn(pjd->V,nv+1-off,1.0/norm);
1386:         }
1387: #endif
1388:         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1389:           BVInsertVec(pjd->W,nv,r[0]);
1390:           if (sz==2 && !off) BVInsertVec(pjd->W,nv+1,r[1]);
1391:           BVOrthogonalizeColumn(pjd->W,nv,NULL,&norm,&lindep);
1393:           BVScaleColumn(pjd->W,nv,1.0/norm);
1394:           if (sz==2 && !off) {
1395:             BVOrthogonalizeColumn(pjd->W,nv+1,NULL,&norm,&lindep);
1397:             BVScaleColumn(pjd->W,nv+1,1.0/norm);
1398:           }
1399:         }
1400:         bupdated = idx?0:nv;
1401:         nv += sz-off;
1402:       }
1403:       for (k=0;k<nvc;k++) {
1404:         eig[pep->nconv-idx+k] = pep->eigr[k];
1405: #if !defined(PETSC_USE_COMPLEX)
1406:         eigi[pep->nconv-idx+k] = pep->eigi[k];
1407: #endif
1408:       }
1409:       PEPMonitor(pep,pep->its,pep->nconv,eig,eigi,pep->errest,pep->nconv+1);
1410:     }
1411:   }
1412:   if (pjd->ld>1) {
1413:     for (k=0;k<pep->nconv;k++) {
1414:       pep->eigr[k] = eig[k];
1415:       pep->eigi[k] = eigi[k];
1416:     }
1417:     if (pep->nconv>0) PEPJDEigenvectors(pep);
1418:     PetscFree2(pcctx->M,pcctx->ps);
1419:   }
1420:   VecDestroy(&u[0]);
1421:   VecDestroy(&r[0]);
1422:   VecDestroy(&p[0]);
1423: #if !defined (PETSC_USE_COMPLEX)
1424:   VecDestroy(&u[1]);
1425:   VecDestroy(&r[1]);
1426:   VecDestroy(&p[1]);
1427: #endif
1428:   KSPSetTolerances(ksp,rtol,abstol,dtol,maxits);
1429:   KSPSetPC(ksp,pcctx->pc);
1430:   VecDestroy(&pcctx->Bp[0]);
1431:   VecDestroy(&pcctx->Bp[1]);
1432:   MatShellGetContext(pjd->Pshell,&matctx);
1433:   MatDestroy(&matctx->Pr);
1434:   MatDestroy(&matctx->Pi);
1435:   MatDestroy(&pjd->Pshell);
1436:   MatDestroy(&pcctx->PPr);
1437:   PCDestroy(&pcctx->pc);
1438:   PetscFree(pcctx);
1439:   PetscFree(matctx);
1440:   PCDestroy(&pjd->pcshell);
1441:   PetscFree3(eig,eigi,res);
1442:   VecDestroy(&pjd->vtempl);
1443:   return 0;
1444: }

1446: PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
1447: {
1448:   PEP_JD *pjd = (PEP_JD*)pep->data;

1450:   if (keep==PETSC_DEFAULT) pjd->keep = 0.5;
1451:   else {
1453:     pjd->keep = keep;
1454:   }
1455:   return 0;
1456: }

1458: /*@
1459:    PEPJDSetRestart - Sets the restart parameter for the Jacobi-Davidson
1460:    method, in particular the proportion of basis vectors that must be kept
1461:    after restart.

1463:    Logically Collective on pep

1465:    Input Parameters:
1466: +  pep  - the eigenproblem solver context
1467: -  keep - the number of vectors to be kept at restart

1469:    Options Database Key:
1470: .  -pep_jd_restart - Sets the restart parameter

1472:    Notes:
1473:    Allowed values are in the range [0.1,0.9]. The default is 0.5.

1475:    Level: advanced

1477: .seealso: PEPJDGetRestart()
1478: @*/
1479: PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
1480: {
1483:   PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
1484:   return 0;
1485: }

1487: PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
1488: {
1489:   PEP_JD *pjd = (PEP_JD*)pep->data;

1491:   *keep = pjd->keep;
1492:   return 0;
1493: }

1495: /*@
1496:    PEPJDGetRestart - Gets the restart parameter used in the Jacobi-Davidson method.

1498:    Not Collective

1500:    Input Parameter:
1501: .  pep - the eigenproblem solver context

1503:    Output Parameter:
1504: .  keep - the restart parameter

1506:    Level: advanced

1508: .seealso: PEPJDSetRestart()
1509: @*/
1510: PetscErrorCode PEPJDGetRestart(PEP pep,PetscReal *keep)
1511: {
1514:   PetscUseMethod(pep,"PEPJDGetRestart_C",(PEP,PetscReal*),(pep,keep));
1515:   return 0;
1516: }

1518: PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
1519: {
1520:   PEP_JD *pjd = (PEP_JD*)pep->data;

1522:   if (fix == PETSC_DEFAULT || fix == PETSC_DECIDE) pjd->fix = 0.01;
1523:   else {
1525:     pjd->fix = fix;
1526:   }
1527:   return 0;
1528: }

1530: /*@
1531:    PEPJDSetFix - Sets the threshold for changing the target in the correction
1532:    equation.

1534:    Logically Collective on pep

1536:    Input Parameters:
1537: +  pep - the eigenproblem solver context
1538: -  fix - threshold for changing the target

1540:    Options Database Key:
1541: .  -pep_jd_fix - the fix value

1543:    Note:
1544:    The target in the correction equation is fixed at the first iterations.
1545:    When the norm of the residual vector is lower than the fix value,
1546:    the target is set to the corresponding eigenvalue.

1548:    Level: advanced

1550: .seealso: PEPJDGetFix()
1551: @*/
1552: PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
1553: {
1556:   PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
1557:   return 0;
1558: }

1560: PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
1561: {
1562:   PEP_JD *pjd = (PEP_JD*)pep->data;

1564:   *fix = pjd->fix;
1565:   return 0;
1566: }

1568: /*@
1569:    PEPJDGetFix - Returns the threshold for changing the target in the correction
1570:    equation.

1572:    Not Collective

1574:    Input Parameter:
1575: .  pep - the eigenproblem solver context

1577:    Output Parameter:
1578: .  fix - threshold for changing the target

1580:    Note:
1581:    The target in the correction equation is fixed at the first iterations.
1582:    When the norm of the residual vector is lower than the fix value,
1583:    the target is set to the corresponding eigenvalue.

1585:    Level: advanced

1587: .seealso: PEPJDSetFix()
1588: @*/
1589: PetscErrorCode PEPJDGetFix(PEP pep,PetscReal *fix)
1590: {
1593:   PetscUseMethod(pep,"PEPJDGetFix_C",(PEP,PetscReal*),(pep,fix));
1594:   return 0;
1595: }

1597: PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
1598: {
1599:   PEP_JD *pjd = (PEP_JD*)pep->data;

1601:   pjd->reusepc = reusepc;
1602:   return 0;
1603: }

1605: /*@
1606:    PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
1607:    must be reused or not.

1609:    Logically Collective on pep

1611:    Input Parameters:
1612: +  pep     - the eigenproblem solver context
1613: -  reusepc - the reuse flag

1615:    Options Database Key:
1616: .  -pep_jd_reuse_preconditioner - the reuse flag

1618:    Note:
1619:    The default value is False. If set to True, the preconditioner is built
1620:    only at the beginning, using the target value. Otherwise, it may be rebuilt
1621:    (depending on the fix parameter) at each iteration from the Ritz value.

1623:    Level: advanced

1625: .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
1626: @*/
1627: PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
1628: {
1631:   PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
1632:   return 0;
1633: }

1635: PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
1636: {
1637:   PEP_JD *pjd = (PEP_JD*)pep->data;

1639:   *reusepc = pjd->reusepc;
1640:   return 0;
1641: }

1643: /*@
1644:    PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.

1646:    Not Collective

1648:    Input Parameter:
1649: .  pep - the eigenproblem solver context

1651:    Output Parameter:
1652: .  reusepc - the reuse flag

1654:    Level: advanced

1656: .seealso: PEPJDSetReusePreconditioner()
1657: @*/
1658: PetscErrorCode PEPJDGetReusePreconditioner(PEP pep,PetscBool *reusepc)
1659: {
1662:   PetscUseMethod(pep,"PEPJDGetReusePreconditioner_C",(PEP,PetscBool*),(pep,reusepc));
1663:   return 0;
1664: }

1666: PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
1667: {
1668:   PEP_JD *pjd = (PEP_JD*)pep->data;

1670:   if (mmidx == PETSC_DEFAULT || mmidx == PETSC_DECIDE) {
1671:     if (pjd->mmidx != 1) pep->state = PEP_STATE_INITIAL;
1672:     pjd->mmidx = 1;
1673:   } else {
1675:     if (pjd->mmidx != mmidx) pep->state = PEP_STATE_INITIAL;
1676:     pjd->mmidx = mmidx;
1677:   }
1678:   return 0;
1679: }

1681: /*@
1682:    PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.

1684:    Logically Collective on pep

1686:    Input Parameters:
1687: +  pep   - the eigenproblem solver context
1688: -  mmidx - maximum minimality index

1690:    Options Database Key:
1691: .  -pep_jd_minimality_index - the minimality index value

1693:    Note:
1694:    The default value is equal to the degree of the polynomial. A smaller value
1695:    can be used if the wanted eigenvectors are known to be linearly independent.

1697:    Level: advanced

1699: .seealso: PEPJDGetMinimalityIndex()
1700: @*/
1701: PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
1702: {
1705:   PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
1706:   return 0;
1707: }

1709: PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
1710: {
1711:   PEP_JD *pjd = (PEP_JD*)pep->data;

1713:   *mmidx = pjd->mmidx;
1714:   return 0;
1715: }

1717: /*@
1718:    PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
1719:    index.

1721:    Not Collective

1723:    Input Parameter:
1724: .  pep - the eigenproblem solver context

1726:    Output Parameter:
1727: .  mmidx - minimality index

1729:    Level: advanced

1731: .seealso: PEPJDSetMinimalityIndex()
1732: @*/
1733: PetscErrorCode PEPJDGetMinimalityIndex(PEP pep,PetscInt *mmidx)
1734: {
1737:   PetscUseMethod(pep,"PEPJDGetMinimalityIndex_C",(PEP,PetscInt*),(pep,mmidx));
1738:   return 0;
1739: }

1741: PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
1742: {
1743:   PEP_JD *pjd = (PEP_JD*)pep->data;

1745:   switch (proj) {
1746:     case PEP_JD_PROJECTION_HARMONIC:
1747:     case PEP_JD_PROJECTION_ORTHOGONAL:
1748:       if (pjd->proj != proj) {
1749:         pep->state = PEP_STATE_INITIAL;
1750:         pjd->proj = proj;
1751:       }
1752:       break;
1753:     default:
1754:       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'proj' value");
1755:   }
1756:   return 0;
1757: }

1759: /*@
1760:    PEPJDSetProjection - Sets the type of projection to be used in the Jacobi-Davidson solver.

1762:    Logically Collective on pep

1764:    Input Parameters:
1765: +  pep  - the eigenproblem solver context
1766: -  proj - the type of projection

1768:    Options Database Key:
1769: .  -pep_jd_projection - the projection type, either orthogonal or harmonic

1771:    Level: advanced

1773: .seealso: PEPJDGetProjection()
1774: @*/
1775: PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
1776: {
1779:   PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
1780:   return 0;
1781: }

1783: PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
1784: {
1785:   PEP_JD *pjd = (PEP_JD*)pep->data;

1787:   *proj = pjd->proj;
1788:   return 0;
1789: }

1791: /*@
1792:    PEPJDGetProjection - Returns the type of projection used by the Jacobi-Davidson solver.

1794:    Not Collective

1796:    Input Parameter:
1797: .  pep - the eigenproblem solver context

1799:    Output Parameter:
1800: .  proj - the type of projection

1802:    Level: advanced

1804: .seealso: PEPJDSetProjection()
1805: @*/
1806: PetscErrorCode PEPJDGetProjection(PEP pep,PEPJDProjection *proj)
1807: {
1810:   PetscUseMethod(pep,"PEPJDGetProjection_C",(PEP,PEPJDProjection*),(pep,proj));
1811:   return 0;
1812: }

1814: PetscErrorCode PEPSetFromOptions_JD(PEP pep,PetscOptionItems *PetscOptionsObject)
1815: {
1816:   PetscBool       flg,b1;
1817:   PetscReal       r1;
1818:   PetscInt        i1;
1819:   PEPJDProjection proj;

1821:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP JD Options");

1823:     PetscOptionsReal("-pep_jd_restart","Proportion of vectors kept after restart","PEPJDSetRestart",0.5,&r1,&flg);
1824:     if (flg) PEPJDSetRestart(pep,r1);

1826:     PetscOptionsReal("-pep_jd_fix","Tolerance for changing the target in the correction equation","PEPJDSetFix",0.01,&r1,&flg);
1827:     if (flg) PEPJDSetFix(pep,r1);

1829:     PetscOptionsBool("-pep_jd_reuse_preconditioner","Whether to reuse the preconditioner","PEPJDSetReusePreconditoiner",PETSC_FALSE,&b1,&flg);
1830:     if (flg) PEPJDSetReusePreconditioner(pep,b1);

1832:     PetscOptionsInt("-pep_jd_minimality_index","Maximum allowed minimality index","PEPJDSetMinimalityIndex",1,&i1,&flg);
1833:     if (flg) PEPJDSetMinimalityIndex(pep,i1);

1835:     PetscOptionsEnum("-pep_jd_projection","Type of projection","PEPJDSetProjection",PEPJDProjectionTypes,(PetscEnum)PEP_JD_PROJECTION_HARMONIC,(PetscEnum*)&proj,&flg);
1836:     if (flg) PEPJDSetProjection(pep,proj);

1838:   PetscOptionsHeadEnd();
1839:   return 0;
1840: }

1842: PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
1843: {
1844:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1845:   PetscBool      isascii;

1847:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1848:   if (isascii) {
1849:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*pjd->keep));
1850:     PetscViewerASCIIPrintf(viewer,"  threshold for changing the target in the correction equation (fix): %g\n",(double)pjd->fix);
1851:     PetscViewerASCIIPrintf(viewer,"  projection type: %s\n",PEPJDProjectionTypes[pjd->proj]);
1852:     PetscViewerASCIIPrintf(viewer,"  maximum allowed minimality index: %" PetscInt_FMT "\n",pjd->mmidx);
1853:     if (pjd->reusepc) PetscViewerASCIIPrintf(viewer,"  reusing the preconditioner\n");
1854:   }
1855:   return 0;
1856: }

1858: PetscErrorCode PEPSetDefaultST_JD(PEP pep)
1859: {
1860:   KSP            ksp;

1862:   if (!((PetscObject)pep->st)->type_name) {
1863:     STSetType(pep->st,STPRECOND);
1864:     STPrecondSetKSPHasMat(pep->st,PETSC_TRUE);
1865:   }
1866:   STSetTransform(pep->st,PETSC_FALSE);
1867:   STGetKSP(pep->st,&ksp);
1868:   if (!((PetscObject)ksp)->type_name) {
1869:     KSPSetType(ksp,KSPBCGSL);
1870:     KSPSetTolerances(ksp,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,100);
1871:   }
1872:   return 0;
1873: }

1875: PetscErrorCode PEPReset_JD(PEP pep)
1876: {
1877:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1878:   PetscInt       i;

1880:   for (i=0;i<pep->nmat;i++) BVDestroy(pjd->TV+i);
1881:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) BVDestroy(&pjd->W);
1882:   if (pjd->ld>1) {
1883:     BVDestroy(&pjd->V);
1884:     for (i=0;i<pep->nmat;i++) BVDestroy(pjd->AX+i);
1885:     BVDestroy(&pjd->N[0]);
1886:     BVDestroy(&pjd->N[1]);
1887:     PetscFree3(pjd->XpX,pjd->T,pjd->Tj);
1888:   }
1889:   PetscFree2(pjd->TV,pjd->AX);
1890:   return 0;
1891: }

1893: PetscErrorCode PEPDestroy_JD(PEP pep)
1894: {
1895:   PetscFree(pep->data);
1896:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",NULL);
1897:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",NULL);
1898:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",NULL);
1899:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",NULL);
1900:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",NULL);
1901:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",NULL);
1902:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",NULL);
1903:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",NULL);
1904:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",NULL);
1905:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",NULL);
1906:   return 0;
1907: }

1909: SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
1910: {
1911:   PEP_JD         *pjd;

1913:   PetscNew(&pjd);
1914:   pep->data = (void*)pjd;

1916:   pep->lineariz = PETSC_FALSE;
1917:   pjd->fix      = 0.01;
1918:   pjd->mmidx    = 0;

1920:   pep->ops->solve          = PEPSolve_JD;
1921:   pep->ops->setup          = PEPSetUp_JD;
1922:   pep->ops->setfromoptions = PEPSetFromOptions_JD;
1923:   pep->ops->destroy        = PEPDestroy_JD;
1924:   pep->ops->reset          = PEPReset_JD;
1925:   pep->ops->view           = PEPView_JD;
1926:   pep->ops->setdefaultst   = PEPSetDefaultST_JD;

1928:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",PEPJDSetRestart_JD);
1929:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",PEPJDGetRestart_JD);
1930:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",PEPJDSetFix_JD);
1931:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",PEPJDGetFix_JD);
1932:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",PEPJDSetReusePreconditioner_JD);
1933:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",PEPJDGetReusePreconditioner_JD);
1934:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",PEPJDSetMinimalityIndex_JD);
1935:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",PEPJDGetMinimalityIndex_JD);
1936:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",PEPJDSetProjection_JD);
1937:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",PEPJDGetProjection_JD);
1938:   return 0;
1939: }