Actual source code: pjd.c

slepc-3.12.2 2020-01-13
Report Typos and Errors
  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.

 19:    References:

 22:            with support for non-monomial bases and deflation", BIT (in
 23:            press), 2019.

 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>    /*I "slepcpep.h" I*/
 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"
 43:   "   journal = \"{BIT} Numer. Math.\",\n"
 44:   "   volume = \"IP\",\n"
 45:   "   number = \"-\",\n"
 46:   "   pages = \"1--24\",\n"
 47:   "   year = \"2019,\"\n"
 48:   "   doi = \"https://doi.org/10.1007/s10543-019-00778-z\"\n"
 49:   "}\n";

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

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

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

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

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

128: PetscErrorCode PEPSetUp_JD(PEP pep)
129: {
131:   PEP_JD         *pjd = (PEP_JD*)pep->data;
132:   PetscBool      isprecond,flg;
133:   PetscInt       i;

136:   pep->lineariz = PETSC_FALSE;
137:   PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
138:   if (!pep->max_it) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
139:   if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;
140:   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");

142:   PetscObjectTypeCompare((PetscObject)pep->st,STPRECOND,&isprecond);
143:   if (!isprecond) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"JD only works with PRECOND spectral transformation");

145:   STGetTransform(pep->st,&flg);
146:   if (flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag unset, see STSetTransform()");

148:   if (!pjd->mmidx) pjd->mmidx = pep->nmat-1;
149:   pjd->mmidx = PetscMin(pjd->mmidx,pep->nmat-1);
150:   if (!pjd->keep) pjd->keep = 0.5;
151:   PEPBasisCoefficients(pep,pep->pbc);
152:   PEPAllocateSolution(pep,0);
153:   PEPSetWorkVecs(pep,5);
154:   pjd->ld = pep->nev;
155: #if !defined (PETSC_USE_COMPLEX)
156:   pjd->ld++;
157: #endif
158:   PetscMalloc2(pep->nmat,&pjd->TV,pep->nmat,&pjd->AX);
159:   for (i=0;i<pep->nmat;i++) {
160:     PEPJDDuplicateBasis(pep,pjd->TV+i);
161:   }
162:   if (pjd->ld>1) {
163:     PEPJDDuplicateBasis(pep,&pjd->V);
164:     BVSetFromOptions(pjd->V);
165:     for (i=0;i<pep->nmat;i++) {
166:       BVDuplicateResize(pep->V,pjd->ld-1,pjd->AX+i);
167:     }
168:     BVDuplicateResize(pep->V,pjd->ld-1,pjd->N);
169:     BVDuplicateResize(pep->V,pjd->ld-1,pjd->N+1);
170:     pjd->X = pep->V;
171:     PetscCalloc3((pjd->ld)*(pjd->ld),&pjd->XpX,pep->ncv*pep->ncv,&pjd->T,pjd->ld*pjd->ld*pep->nmat,&pjd->Tj);
172:   } else pjd->V = pep->V;
173:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) { PEPJDDuplicateBasis(pep,&pjd->W); }
174:   else pjd->W = pjd->V;
175:   DSSetType(pep->ds,DSPEP);
176:   DSPEPSetDegree(pep->ds,pep->nmat-1);
177:   if (pep->basis!=PEP_BASIS_MONOMIAL) {
178:     DSPEPSetCoefficients(pep->ds,pep->pbc);
179:   }
180:   DSAllocate(pep->ds,pep->ncv);
181:   return(0);
182: }

184: /*
185:    Updates columns (low to (high-1)) of TV[i]
186: */
187: static PetscErrorCode PEPJDUpdateTV(PEP pep,PetscInt low,PetscInt high,Vec *w)
188: {
190:   PEP_JD         *pjd = (PEP_JD*)pep->data;
191:   PetscInt       pp,col,i,nloc,nconv;
192:   Vec            v1,v2,t1,t2;
193:   PetscScalar    *array1,*array2,*x2,*xx,*N,*Np,*y2=NULL,zero=0.0,sone=1.0,*pT,fact,*psc;
194:   PetscReal      *cg,*ca,*cb;
195:   PetscMPIInt    rk,np;
196:   PetscBLASInt   n_,ld_,one=1;
197:   Mat            T;
198:   BV             pbv;

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

288: /*
289:    RRQR of X. Xin*P=Xou*R. Rank of R is rk
290: */
291: static PetscErrorCode PEPJDOrthogonalize(PetscInt row,PetscInt col,PetscScalar *X,PetscInt ldx,PetscInt *rk,PetscInt *P,PetscScalar *R,PetscInt ldr)
292: {
293: #if defined(SLEPC_MISSING_LAPACK_GEQP3) || defined(PETSC_MISSING_LAPACK_ORGQR)
295:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEQP3/QRGQR - Lapack routines are unavailable");
296: #else
298:   PetscInt       i,j,n,r;
299:   PetscBLASInt   row_,col_,ldx_,*p,lwork,info,n_;
300:   PetscScalar    *tau,*work;
301:   PetscReal      tol,*rwork;

304:   PetscBLASIntCast(row,&row_);
305:   PetscBLASIntCast(col,&col_);
306:   PetscBLASIntCast(ldx,&ldx_);
307:   n = PetscMin(row,col);
308:   PetscBLASIntCast(n,&n_);
309:   lwork = 3*col_+1;
310:   PetscMalloc4(col,&p,n,&tau,lwork,&work,2*col,&rwork);
311:   for (i=1;i<col;i++) p[i] = 0;
312:   p[0] = 1;

314:   /* rank revealing QR */
315: #if defined(PETSC_USE_COMPLEX)
316:   PetscStackCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,rwork,&info));
317: #else
318:   PetscStackCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,&info));
319: #endif
320:   SlepcCheckLapackInfo("geqp3",info);
321:   if (P) for (i=0;i<col;i++) P[i] = p[i]-1;

323:   /* rank computation */
324:   tol = PetscMax(row,col)*PETSC_MACHINE_EPSILON*PetscAbsScalar(X[0]);
325:   r = 1;
326:   for (i=1;i<n;i++) {
327:     if (PetscAbsScalar(X[i+ldx*i])>tol) r++;
328:     else break;
329:   }
330:   if (rk) *rk=r;

332:   /* copy upper triangular matrix if requested */
333:   if (R) {
334:      for (i=0;i<r;i++) {
335:        PetscArrayzero(R+i*ldr,r);
336:        for (j=0;j<=i;j++) R[i*ldr+j] = X[i*ldx+j];
337:      }
338:   }
339:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&row_,&n_,&n_,X,&ldx_,tau,work,&lwork,&info));
340:   SlepcCheckLapackInfo("orgqr",info);
341:   PetscFree4(p,tau,work,rwork);
342:   return(0);
343: #endif
344: }

346: /*
347:    Application of extended preconditioner
348: */
349: static PetscErrorCode PEPJDExtendedPCApply(PC pc,Vec x,Vec y)
350: {
351:   PetscInt          i,j,nloc,n,ld=0;
352:   PetscMPIInt       np;
353:   Vec               tx,ty;
354:   PEP_JD_PCSHELL    *ctx;
355:   PetscErrorCode    ierr;
356:   const PetscScalar *array1;
357:   PetscScalar       *x2=NULL,*t=NULL,*ps=NULL,*array2,zero=0.0,sone=1.0;
358:   PetscBLASInt      one=1.0,ld_,n_,ncv_;
359:   PEP_JD            *pjd=NULL;

362:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&np);
363:   PCShellGetContext(pc,(void**)&ctx);
364:   n  = ctx->n;
365:   if (n) {
366:     pjd = (PEP_JD*)ctx->pep->data;
367:     ps = ctx->ps;
368:     ld = pjd->ld;
369:     PetscMalloc2(n,&x2,n,&t);
370:     VecGetLocalSize(ctx->work[0],&nloc);
371:     VecGetArrayRead(x,&array1);
372:     for (i=0;i<n;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
373:     VecRestoreArrayRead(x,&array1);
374:   }

376:   /* y = B\x apply PC */
377:   tx = ctx->work[0];
378:   ty = ctx->work[1];
379:   VecGetArrayRead(x,&array1);
380:   VecPlaceArray(tx,array1);
381:   VecGetArray(y,&array2);
382:   VecPlaceArray(ty,array2);
383:   PCApply(ctx->pc,tx,ty);
384:   if (n) {
385:     PetscBLASIntCast(ld,&ld_);
386:     PetscBLASIntCast(n,&n_);
387:     for (i=0;i<n;i++) {
388:       t[i] = 0.0;
389:       for (j=0;j<n;j++) t[i] += ctx->M[i+j*ld]*x2[j];
390:     }
391:     if (pjd->midx==1) {
392:       PetscBLASIntCast(ctx->pep->ncv,&ncv_);
393:       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] -= ctx->theta;
394:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,pjd->T,&ncv_,t,&one,&zero,x2,&one));
395:       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] += ctx->theta;
396:       for (i=0;i<n;i++) array2[nloc+i] = x2[i];
397:       for (i=0;i<n;i++) x2[i] = -t[i];
398:     } else {
399:       for (i=0;i<n;i++) array2[nloc+i] = t[i];
400:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,ps,&ld_,t,&one,&zero,x2,&one));
401:     }
402:     for (i=0;i<n;i++) array2[nloc+i] /= PetscSqrtReal(np);
403:     BVSetActiveColumns(pjd->X,0,n);
404:     BVMultVec(pjd->X,-1.0,1.0,ty,x2);
405:     PetscFree2(x2,t);
406:   }
407:   VecResetArray(tx);
408:   VecResetArray(ty);
409:   VecRestoreArrayRead(x,&array1);
410:   VecRestoreArray(y,&array2);
411:   return(0);
412: }

414: /*
415:    Application of shell preconditioner:
416:       y = B\x - eta*B\p,  with eta = (u'*B\x)/(u'*B\p)
417: */
418: static PetscErrorCode PCShellApply_PEPJD(PC pc,Vec x,Vec y)
419: {
421:   PetscScalar    rr,eta;
422:   PEP_JD_PCSHELL *ctx;
423:   PetscInt       sz;
424:   const Vec      *xs,*ys;
425: #if !defined(PETSC_USE_COMPLEX)
426:   PetscScalar    rx,xr,xx;
427: #endif

430:   PCShellGetContext(pc,(void**)&ctx);
431:   VecCompGetSubVecs(x,&sz,&xs);
432:   VecCompGetSubVecs(y,NULL,&ys);
433:   /* y = B\x apply extended PC */
434:   PEPJDExtendedPCApply(pc,xs[0],ys[0]);
435: #if !defined(PETSC_USE_COMPLEX)
436:   if (sz==2) {
437:     PEPJDExtendedPCApply(pc,xs[1],ys[1]);
438:   }
439: #endif

441:   /* Compute eta = u'*y / u'*Bp */
442:   VecDot(ys[0],ctx->u[0],&rr);
443:   eta  = -rr*ctx->gamma[0];
444: #if !defined(PETSC_USE_COMPLEX)
445:   if (sz==2) {
446:     VecDot(ys[0],ctx->u[1],&xr);
447:     VecDot(ys[1],ctx->u[0],&rx);
448:     VecDot(ys[1],ctx->u[1],&xx);
449:     eta += -ctx->gamma[0]*xx-ctx->gamma[1]*(-xr+rx);
450:   }
451: #endif
452:   eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];

454:   /* y = y - eta*Bp */
455:   VecAXPY(ys[0],eta,ctx->Bp[0]);
456: #if !defined(PETSC_USE_COMPLEX)
457:   if (sz==2) {
458:     VecAXPY(ys[1],eta,ctx->Bp[1]);
459:     eta = -ctx->gamma[1]*(rr+xx)+ctx->gamma[0]*(-xr+rx);
460:     eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
461:     VecAXPY(ys[0],eta,ctx->Bp[1]);
462:     VecAXPY(ys[1],-eta,ctx->Bp[0]);
463:   }
464: #endif
465:   return(0);
466: }

468: static PetscErrorCode PEPJDCopyToExtendedVec(PEP pep,Vec v,PetscScalar *a,PetscInt na,PetscInt off,Vec vex,PetscBool back)
469: {
471:   PetscMPIInt    np,rk,count;
472:   PetscScalar    *array1,*array2;
473:   PetscInt       nloc;

476:   MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk);
477:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
478:   BVGetSizes(pep->V,&nloc,NULL,NULL);
479:   if (v) {
480:     VecGetArray(v,&array1);
481:     VecGetArray(vex,&array2);
482:     if (back) {
483:       PetscArraycpy(array1,array2,nloc);
484:     } else {
485:       PetscArraycpy(array2,array1,nloc);
486:     }
487:     VecRestoreArray(v,&array1);
488:     VecRestoreArray(vex,&array2);
489:   }
490:   if (a) {
491:     VecGetArray(vex,&array2);
492:     if (back) {
493:       PetscArraycpy(a,array2+nloc+off,na);
494:       PetscMPIIntCast(na,&count);
495:       MPI_Bcast(a,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
496:     } else {
497:       PetscArraycpy(array2+nloc+off,a,na);
498:       PetscMPIIntCast(na,&count);
499:       MPI_Bcast(array2+nloc+off,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
500:     }
501:     VecRestoreArray(vex,&array2);
502:   }
503:   return(0);
504: }

506: /* Computes Phi^hat(lambda) times a vector or its derivative (depends on beval)
507:      if no vector is provided returns a matrix
508:  */
509: static PetscErrorCode PEPJDEvaluateHatBasis(PEP pep,PetscInt n,PetscScalar *H,PetscInt ldh,PetscScalar *beval,PetscScalar *t,PetscInt idx,PetscScalar *qpp,PetscScalar *qp,PetscScalar *q)
510: {
512:   PetscInt       j,i;
513:   PetscBLASInt   n_,ldh_,one=1;
514:   PetscReal      *a,*b,*g;
515:   PetscScalar    sone=1.0,zero=0.0;

518:   a = pep->pbc; b=a+pep->nmat; g=b+pep->nmat;
519:   PetscBLASIntCast(n,&n_);
520:   PetscBLASIntCast(ldh,&ldh_);
521:   if (idx<1) {
522:     PetscArrayzero(q,t?n:n*n);
523:   } else if (idx==1) {
524:     if (t) {for (j=0;j<n;j++) q[j] = t[j]*beval[idx-1]/a[0];}
525:     else {
526:       PetscArrayzero(q,n*n);
527:       for (j=0;j<n;j++) q[(j+1)*n] = beval[idx-1]/a[0];
528:     }
529:   } else {
530:     if (t) {
531:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,H,&ldh_,qp,&one,&zero,q,&one));
532:       for (j=0;j<n;j++) {
533:         q[j] += beval[idx-1]*t[j]-b[idx-1]*qp[j]-g[idx-1]*qpp[j];
534:         q[j] /= a[idx-1];
535:       }
536:     } else {
537:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,H,&ldh_,qp,&n_,&zero,q,&n_));
538:       for (j=0;j<n;j++) {
539:         q[j+n*j] += beval[idx-1];
540:         for (i=0;i<n;i++) {
541:           q[i+n*j] += -b[idx-1]*qp[j*n+i]-g[idx-1]*qpp[j*n+i];
542:           q[i+n*j] /= a[idx-1];
543:         }
544:       }
545:     }
546:   }
547:   return(0);
548: }

550: static PetscErrorCode PEPJDComputeResidual(PEP pep,PetscBool derivative,PetscInt sz,Vec *u,PetscScalar *theta,Vec *p,Vec *work)
551: {
552:   PEP_JD         *pjd = (PEP_JD*)pep->data;
554:   PetscMPIInt    rk,np,count;
555:   Vec            tu,tp,w;
556:   PetscScalar    *dval,*dvali,*array1,*array2,*x2=NULL,*y2,*qj=NULL,*tt=NULL,*xx=NULL,*xxi=NULL,sone=1.0;
557:   PetscInt       i,j,nconv,nloc;
558:   PetscBLASInt   n,ld,one=1;
559: #if !defined(PETSC_USE_COMPLEX)
560:   Vec            tui=NULL,tpi=NULL;
561:   PetscScalar    *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
562: #endif

565:   nconv = pjd->nlock;
566:   if (!nconv) {
567:     PetscMalloc1(2*sz*pep->nmat,&dval);
568:   } else {
569:     PetscMalloc5(2*pep->nmat,&dval,2*nconv,&xx,nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*pep->nmat,&qj);
570:     MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk);
571:     MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
572:     BVGetSizes(pep->V,&nloc,NULL,NULL);
573:     VecGetArray(u[0],&array1);
574:     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]*PetscSqrtReal(np);
575:     VecRestoreArray(u[0],&array1);
576: #if !defined(PETSC_USE_COMPLEX)
577:     if (sz==2) {
578:       x2i = x2+nconv;
579:       VecGetArray(u[1],&arrayi1);
580:       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]*PetscSqrtReal(np);
581:       VecRestoreArray(u[1],&arrayi1);
582:     }
583: #endif
584:   }
585:   dvali = dval+pep->nmat;
586:   tu = work[0];
587:   tp = work[1];
588:   w  = work[2];
589:   VecGetArray(u[0],&array1);
590:   VecPlaceArray(tu,array1);
591:   VecGetArray(p[0],&array2);
592:   VecPlaceArray(tp,array2);
593:   VecSet(tp,0.0);
594: #if !defined(PETSC_USE_COMPLEX)
595:   if (sz==2) {
596:     tui = work[3];
597:     tpi = work[4];
598:     VecGetArray(u[1],&arrayi1);
599:     VecPlaceArray(tui,arrayi1);
600:     VecGetArray(p[1],&arrayi2);
601:     VecPlaceArray(tpi,arrayi2);
602:     VecSet(tpi,0.0);
603:   }
604: #endif
605:   if (derivative) {
606:     PEPEvaluateBasisDerivative(pep,theta[0],theta[1],dval,dvali);
607:   } else {
608:     PEPEvaluateBasis(pep,theta[0],theta[1],dval,dvali);
609:   }
610:   for (i=derivative?1:0;i<pep->nmat;i++) {
611:     MatMult(pep->A[i],tu,w);
612:     VecAXPY(tp,dval[i],w);
613: #if !defined(PETSC_USE_COMPLEX)
614:     if (sz==2) {
615:       VecAXPY(tpi,dvali[i],w);
616:       MatMult(pep->A[i],tui,w);
617:       VecAXPY(tpi,dval[i],w);
618:       VecAXPY(tp,-dvali[i],w);
619:     }
620: #endif
621:   }
622:   if (nconv) {
623:     for (i=0;i<pep->nmat;i++) {
624:       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);
625:     }
626: #if !defined(PETSC_USE_COMPLEX)
627:     if (sz==2) {
628:       qji = qj+nconv*pep->nmat;
629:       qq = qji+nconv*pep->nmat;
630:       for (i=0;i<pep->nmat;i++) {
631:         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);
632:       }
633:       for (i=0;i<nconv*pep->nmat;i++) qj[i] -= qji[i];
634:       for (i=0;i<pep->nmat;i++) {
635:         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);
636:         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);
637:       }
638:       for (i=0;i<nconv*pep->nmat;i++) qji[i] += qq[i];
639:       for (i=derivative?2:1;i<pep->nmat;i++) {
640:         BVMultVec(pjd->AX[i],1.0,1.0,tpi,qji+i*nconv);
641:       }
642:     }
643: #endif
644:     for (i=derivative?2:1;i<pep->nmat;i++) {
645:       BVMultVec(pjd->AX[i],1.0,1.0,tp,qj+i*nconv);
646:     }

648:     /* extended vector part */
649:     BVSetActiveColumns(pjd->X,0,nconv);
650:     BVDotVec(pjd->X,tu,xx);
651:     xxi = xx+nconv;
652: #if !defined(PETSC_USE_COMPLEX)
653:     if (sz==2) { BVDotVec(pjd->X,tui,xxi); }
654: #endif
655:     if (sz==1) { PetscArrayzero(xxi,nconv); }
656:     if (rk==np-1) {
657:       PetscBLASIntCast(nconv,&n);
658:       PetscBLASIntCast(pjd->ld,&ld);
659:       y2  = array2+nloc;
660:       PetscArrayzero(y2,nconv);
661:       for (j=derivative?1:0;j<pjd->midx;j++) {
662:         for (i=0;i<nconv;i++) tt[i] = dval[j]*xx[i]-dvali[j]*xxi[i];
663:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
664:         PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
665:       }
666:       for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
667: #if !defined(PETSC_USE_COMPLEX)
668:       if (sz==2) {
669:         y2i = arrayi2+nloc;
670:         PetscArrayzero(y2i,nconv);
671:         for (j=derivative?1:0;j<pjd->midx;j++) {
672:           for (i=0;i<nconv;i++) tt[i] = dval[j]*xxi[i]+dvali[j]*xx[i];
673:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
674:           PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
675:         }
676:         for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
677:       }
678: #endif
679:     }
680:     PetscMPIIntCast(nconv,&count);
681:     MPI_Bcast(array2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
682: #if !defined(PETSC_USE_COMPLEX)
683:     if (sz==2) {
684:       MPI_Bcast(arrayi2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
685:     }
686: #endif
687:   }
688:   if (nconv) {
689:     PetscFree5(dval,xx,tt,x2,qj);
690:   } else {
691:     PetscFree(dval);
692:   }
693:   VecResetArray(tu);
694:   VecRestoreArray(u[0],&array1);
695:   VecResetArray(tp);
696:   VecRestoreArray(p[0],&array2);
697: #if !defined(PETSC_USE_COMPLEX)
698:   if (sz==2) {
699:     VecResetArray(tui);
700:     VecRestoreArray(u[1],&arrayi1);
701:     VecResetArray(tpi);
702:     VecRestoreArray(p[1],&arrayi2);
703:   }
704: #endif
705:   return(0);
706: }

708: static PetscErrorCode PEPJDProcessInitialSpace(PEP pep,Vec *w)
709: {
710:   PEP_JD         *pjd = (PEP_JD*)pep->data;
712:   PetscScalar    *tt,target[2];
713:   Vec            vg,wg;
714:   PetscInt       i;
715:   PetscReal      norm;

718:   PetscMalloc1(pjd->ld-1,&tt);
719:   if (pep->nini==0) {
720:     BVSetRandomColumn(pjd->V,0);
721:     for (i=0;i<pjd->ld-1;i++) tt[i] = 0.0;
722:     BVGetColumn(pjd->V,0,&vg);
723:     PEPJDCopyToExtendedVec(pep,NULL,tt,pjd->ld-1,0,vg,PETSC_FALSE);
724:     BVRestoreColumn(pjd->V,0,&vg);
725:     BVNormColumn(pjd->V,0,NORM_2,&norm);
726:     BVScaleColumn(pjd->V,0,1.0/norm);
727:     if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
728:       BVGetColumn(pjd->V,0,&vg);
729:       BVGetColumn(pjd->W,0,&wg);
730:       VecSet(wg,0.0);
731:       target[0] = pep->target; target[1] = 0.0;
732:       PEPJDComputeResidual(pep,PETSC_TRUE,1,&vg,target,&wg,w);
733:       BVRestoreColumn(pjd->W,0,&wg);
734:       BVRestoreColumn(pjd->V,0,&vg);
735:       BVNormColumn(pjd->W,0,NORM_2,&norm);
736:       BVScaleColumn(pjd->W,0,1.0/norm);
737:     }
738:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for initial vectors not implemented yet");
739:   PetscFree(tt);
740:   return(0);
741: }

743: static PetscErrorCode PEPJDShellMatMult(Mat P,Vec x,Vec y)
744: {
745:   PetscErrorCode  ierr;
746:   PEP_JD_MATSHELL *matctx;
747:   PEP_JD          *pjd;
748:   PetscInt        i,j,nconv,nloc,nmat,ldt,ncv,sz;
749:   Vec             tx,ty;
750:   const Vec       *xs,*ys;
751:   PetscScalar     *array1,*array2,*x2=NULL,*y2,*tt=NULL,*xx=NULL,*xxi,theta[2],sone=1.0,*qj,*val,*vali=NULL;
752:   PetscBLASInt    n,ld,one=1;
753:   PetscMPIInt     np;
754: #if !defined(PETSC_USE_COMPLEX)
755:   Vec             txi=NULL,tyi=NULL;
756:   PetscScalar     *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
757: #endif

760:   MPI_Comm_size(PetscObjectComm((PetscObject)P),&np);
761:   MatShellGetContext(P,(void**)&matctx);
762:   pjd   = (PEP_JD*)(matctx->pep->data);
763:   nconv = pjd->nlock;
764:   nmat  = matctx->pep->nmat;
765:   ncv   = matctx->pep->ncv;
766:   ldt   = pjd->ld;
767:   VecCompGetSubVecs(x,&sz,&xs);
768:   VecCompGetSubVecs(y,NULL,&ys);
769:   theta[0] = matctx->theta[0];
770:   theta[1] = (sz==2)?matctx->theta[1]:0.0;
771:   if (nconv>0) {
772:     PetscMalloc5(nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*nmat,&qj,2*nconv,&xx,2*nmat,&val);
773:     BVGetSizes(matctx->pep->V,&nloc,NULL,NULL);
774:     VecGetArray(xs[0],&array1);
775:     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
776:     VecRestoreArray(xs[0],&array1);
777: #if !defined(PETSC_USE_COMPLEX)
778:     if (sz==2) {
779:       x2i = x2+nconv;
780:       VecGetArray(xs[1],&arrayi1);
781:       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]* PetscSqrtReal(np);
782:       VecRestoreArray(xs[1],&arrayi1);
783:     }
784: #endif
785:     vali = val+nmat;
786:   }
787:   tx = matctx->work[0];
788:   ty = matctx->work[1];
789:   VecGetArray(xs[0],&array1);
790:   VecPlaceArray(tx,array1);
791:   VecGetArray(ys[0],&array2);
792:   VecPlaceArray(ty,array2);
793:   MatMult(matctx->Pr,tx,ty);
794: #if !defined(PETSC_USE_COMPLEX)
795:   if (sz==2) {
796:     txi = matctx->work[2];
797:     tyi = matctx->work[3];
798:     VecGetArray(xs[1],&arrayi1);
799:     VecPlaceArray(txi,arrayi1);
800:     VecGetArray(ys[1],&arrayi2);
801:     VecPlaceArray(tyi,arrayi2);
802:     MatMult(matctx->Pr,txi,tyi);
803:     if (theta[1]!=0.0) {
804:       MatMult(matctx->Pi,txi,matctx->work[4]);
805:       VecAXPY(ty,-1.0,matctx->work[4]);
806:       MatMult(matctx->Pi,tx,matctx->work[4]);
807:       VecAXPY(tyi,1.0,matctx->work[4]);
808:     }
809:   }
810: #endif
811:   if (nconv>0) {
812:     PEPEvaluateBasis(matctx->pep,theta[0],theta[1],val,vali);
813:     for (i=0;i<nmat;i++) {
814:       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);
815:     }
816: #if !defined(PETSC_USE_COMPLEX)
817:     if (sz==2) {
818:       qji = qj+nconv*nmat;
819:       qq = qji+nconv*nmat;
820:       for (i=0;i<nmat;i++) {
821:         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);
822:       }
823:       for (i=0;i<nconv*nmat;i++) qj[i] -= qji[i];
824:       for (i=0;i<nmat;i++) {
825:         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);
826:         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);
827:       }
828:       for (i=0;i<nconv*nmat;i++) qji[i] += qq[i];
829:       for (i=1;i<matctx->pep->nmat;i++) {
830:         BVMultVec(pjd->AX[i],1.0,1.0,tyi,qji+i*nconv);
831:       }
832:     }
833: #endif
834:     for (i=1;i<nmat;i++) {
835:       BVMultVec(pjd->AX[i],1.0,1.0,ty,qj+i*nconv);
836:     }

838:     /* extended vector part */
839:     BVSetActiveColumns(pjd->X,0,nconv);
840:     BVDotVec(pjd->X,tx,xx);
841:     xxi = xx+nconv;
842: #if !defined(PETSC_USE_COMPLEX)
843:     if (sz==2) { BVDotVec(pjd->X,txi,xxi); }
844: #endif
845:     if (sz==1) { PetscArrayzero(xxi,nconv); }
846:     PetscBLASIntCast(pjd->nlock,&n);
847:     PetscBLASIntCast(ldt,&ld);
848:     y2 = array2+nloc;
849:     PetscArrayzero(y2,nconv);
850:     for (j=0;j<pjd->midx;j++) {
851:       for (i=0;i<nconv;i++) tt[i] = val[j]*xx[i]-vali[j]*xxi[i];
852:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
853:       PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
854:     }
855: #if !defined(PETSC_USE_COMPLEX)
856:     if (sz==2) {
857:       y2i = arrayi2+nloc;
858:       PetscArrayzero(y2i,nconv);
859:       for (j=0;j<pjd->midx;j++) {
860:         for (i=0;i<nconv;i++) tt[i] = val[j]*xxi[i]+vali[j]*xx[i];
861:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
862:         PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
863:       }
864:       for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
865:     }
866: #endif
867:     for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
868:     PetscFree5(tt,x2,qj,xx,val);
869:   }
870:   VecResetArray(tx);
871:   VecRestoreArray(xs[0],&array1);
872:   VecResetArray(ty);
873:   VecRestoreArray(ys[0],&array2);
874: #if !defined(PETSC_USE_COMPLEX)
875:   if (sz==2) {
876:     VecResetArray(txi);
877:     VecRestoreArray(xs[1],&arrayi1);
878:     VecResetArray(tyi);
879:     VecRestoreArray(ys[1],&arrayi2);
880:   }
881: #endif
882:   return(0);
883: }

885: static PetscErrorCode PEPJDSellMatCreateVecs(Mat A,Vec *right,Vec *left)
886: {
887:   PetscErrorCode  ierr;
888:   PEP_JD_MATSHELL *matctx;
889:   PEP_JD          *pjd;
890:   PetscInt        kspsf=1,i;
891:   Vec             v[2];

894:   MatShellGetContext(A,(void**)&matctx);
895:   pjd   = (PEP_JD*)(matctx->pep->data);
896: #if !defined (PETSC_USE_COMPLEX)
897:   kspsf = 2;
898: #endif
899:   for (i=0;i<kspsf;i++){
900:     BVCreateVec(pjd->V,v+i);
901:   }
902:   if (right) {
903:     VecCreateCompWithVecs(v,kspsf,pjd->vtempl,right);
904:   }
905:   if (left) {
906:     VecCreateCompWithVecs(v,kspsf,pjd->vtempl,left);
907:   }
908:   for (i=0;i<kspsf;i++) {
909:     VecDestroy(&v[i]);
910:   }
911:   return(0);
912: }

914: static PetscErrorCode PEPJDUpdateExtendedPC(PEP pep,PetscScalar theta)
915: {
916: #if defined(PETSC_MISSING_LAPACK_GESVD) || defined(PETSC_MISSING_LAPACK_GETRI) || defined(PETSC_MISSING_LAPACK_GETRF)
918:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD/GETRI/GETRF - Lapack routines are unavailable");
919: #else
921:   PEP_JD         *pjd = (PEP_JD*)pep->data;
922:   PEP_JD_PCSHELL *pcctx;
923:   PetscInt       i,j,k,n=pjd->nlock,ld=pjd->ld,deg=pep->nmat-1;
924:   PetscScalar    *M,*ps,*work,*U,*V,*S,*Sp,*Spp,snone=-1.0,sone=1.0,zero=0.0,*val;
925:   PetscReal      tol,maxeig=0.0,*sg,*rwork;
926:   PetscBLASInt   n_,info,ld_,*p,lw_,rk=0;

929:   if (n) {
930:     PCShellGetContext(pjd->pcshell,(void**)&pcctx);
931:     pcctx->theta = theta;
932:     pcctx->n = n;
933:     M  = pcctx->M;
934:     PetscBLASIntCast(n,&n_);
935:     PetscBLASIntCast(ld,&ld_);
936:     if (pjd->midx==1) {
937:       PetscArraycpy(M,pjd->XpX,ld*ld);
938:       PetscCalloc2(10*n,&work,n,&p);
939:     } else {
940:       ps = pcctx->ps;
941:       PetscCalloc7(2*n*n,&U,3*n*n,&S,n,&sg,10*n,&work,5*n,&rwork,n,&p,deg+1,&val);
942:       V = U+n*n;
943:       /* pseudo-inverse */
944:       for (j=0;j<n;j++) {
945:         for (i=0;i<n;i++) S[n*j+i] = -pjd->T[pep->ncv*j+i];
946:         S[n*j+j] += theta;
947:       }
948:       lw_ = 10*n_;
949: #if !defined (PETSC_USE_COMPLEX)
950:       PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,&info));
951: #else
952:       PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,rwork,&info));
953: #endif
954:       SlepcCheckLapackInfo("gesvd",info);
955:       for (i=0;i<n;i++) maxeig = PetscMax(maxeig,sg[i]);
956:       tol = 10*PETSC_MACHINE_EPSILON*n*maxeig;
957:       for (j=0;j<n;j++) {
958:         if (sg[j]>tol) {
959:           for (i=0;i<n;i++) U[j*n+i] /= sg[j];
960:           rk++;
961:         } else break;
962:       }
963:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&rk,&sone,U,&n_,V,&n_,&zero,ps,&ld_));

965:       /* compute M */
966:       PEPEvaluateBasis(pep,theta,0.0,val,NULL);
967:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&snone,pjd->XpX,&ld_,ps,&ld_,&zero,M,&ld_));
968:       PetscArrayzero(S,2*n*n);
969:       Sp = S+n*n;
970:       for (j=0;j<n;j++) S[j*(n+1)] = 1.0;
971:       for (k=1;k<pjd->midx;k++) {
972:         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];
973:         PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pjd->XpX,&ld_,V,&n_,&zero,U,&n_));
974:         PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+k*ld*ld,&ld_,U,&n_,&sone,M,&ld_));
975:         Spp = Sp; Sp = S;
976:         PEPJDEvaluateHatBasis(pep,n,pjd->T,pep->ncv,val,NULL,k+1,Spp,Sp,S);
977:       }
978:     }
979:     /* inverse */
980:     PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,M,&ld_,p,&info));
981:     SlepcCheckLapackInfo("getrf",info);
982:     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n_,M,&ld_,p,work,&n_,&info));
983:     SlepcCheckLapackInfo("getri",info);
984:     if (pjd->midx==1) {
985:       PetscFree2(work,p);
986:     } else {
987:       PetscFree7(U,S,sg,work,rwork,p,val);
988:     }
989:   }
990:   return(0);
991: #endif
992: }

994: static PetscErrorCode PEPJDMatSetUp(PEP pep,PetscInt sz,PetscScalar *theta)
995: {
996:   PetscErrorCode  ierr;
997:   PEP_JD          *pjd = (PEP_JD*)pep->data;
998:   PEP_JD_MATSHELL *matctx;
999:   PEP_JD_PCSHELL  *pcctx;
1000:   MatStructure    str;
1001:   PetscScalar     *vals,*valsi;
1002:   PetscBool       skipmat=PETSC_FALSE;
1003:   PetscInt        i;
1004:   Mat             Pr=NULL;

1007:   if (sz==2 && theta[1]==0.0) sz = 1;
1008:   MatShellGetContext(pjd->Pshell,(void**)&matctx);
1009:   PCShellGetContext(pjd->pcshell,(void**)&pcctx);
1010:   if (matctx->Pr && matctx->theta[0]==theta[0] && ((!matctx->Pi && sz==1) || (sz==2 && matctx->theta[1]==theta[1]))) {
1011:     if (pcctx->n == pjd->nlock) return(0);
1012:     skipmat = PETSC_TRUE;
1013:   }
1014:   if (!skipmat) {
1015:     PetscMalloc2(pep->nmat,&vals,pep->nmat,&valsi);
1016:     STGetMatStructure(pep->st,&str);
1017:     PEPEvaluateBasis(pep,theta[0],theta[1],vals,valsi);
1018:     if (!matctx->Pr) {
1019:       MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pr);
1020:     } else {
1021:       MatCopy(pep->A[0],matctx->Pr,str);
1022:     }
1023:     for (i=1;i<pep->nmat;i++) {
1024:       MatAXPY(matctx->Pr,vals[i],pep->A[i],str);
1025:     }
1026:     if (!pjd->reusepc) {
1027:       if (pcctx->PPr && sz==2) {
1028:         MatCopy(matctx->Pr,pcctx->PPr,str);
1029:         Pr = pcctx->PPr;
1030:       } else Pr = matctx->Pr;
1031:     }
1032:     matctx->theta[0] = theta[0];
1033: #if !defined(PETSC_USE_COMPLEX)
1034:     if (sz==2) {
1035:       if (!matctx->Pi) {
1036:         MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pi);
1037:       } else {
1038:         MatCopy(pep->A[1],matctx->Pi,str);
1039:       }
1040:       MatScale(matctx->Pi,valsi[1]);
1041:       for (i=2;i<pep->nmat;i++) {
1042:         MatAXPY(matctx->Pi,valsi[i],pep->A[i],str);
1043:       }
1044:       matctx->theta[1] = theta[1];
1045:     }
1046: #endif
1047:     PetscFree2(vals,valsi);
1048:   }
1049:   if (!pjd->reusepc) {
1050:     if (!skipmat) {
1051:       PCSetOperators(pcctx->pc,Pr,Pr);
1052:       PCSetUp(pcctx->pc);
1053:     }
1054:     PEPJDUpdateExtendedPC(pep,theta[0]);
1055:   }
1056:   return(0);
1057: }

1059: static PetscErrorCode PEPJDCreateShellPC(PEP pep,Vec *ww)
1060: {
1061:   PEP_JD          *pjd = (PEP_JD*)pep->data;
1062:   PEP_JD_PCSHELL  *pcctx;
1063:   PEP_JD_MATSHELL *matctx;
1064:   KSP             ksp;
1065:   PetscInt        nloc,mloc,kspsf=1;
1066:   Vec             v[2];
1067:   PetscScalar     target[2];
1068:   PetscErrorCode  ierr;
1069:   Mat             Pr;

1072:   /* Create the reference vector */
1073:   BVGetColumn(pjd->V,0,&v[0]);
1074:   v[1] = v[0];
1075: #if !defined (PETSC_USE_COMPLEX)
1076:   kspsf = 2;
1077: #endif
1078:   VecCreateCompWithVecs(v,kspsf,NULL,&pjd->vtempl);
1079:   BVRestoreColumn(pjd->V,0,&v[0]);
1080:   PetscLogObjectParent((PetscObject)pep,(PetscObject)pjd->vtempl);

1082:   /* Replace preconditioner with one containing projectors */
1083:   PCCreate(PetscObjectComm((PetscObject)pep),&pjd->pcshell);
1084:   PCSetType(pjd->pcshell,PCSHELL);
1085:   PCShellSetName(pjd->pcshell,"PCPEPJD");
1086:   PCShellSetApply(pjd->pcshell,PCShellApply_PEPJD);
1087:   PetscNew(&pcctx);
1088:   PCShellSetContext(pjd->pcshell,pcctx);
1089:   STGetKSP(pep->st,&ksp);
1090:   BVCreateVec(pjd->V,&pcctx->Bp[0]);
1091:   VecDuplicate(pcctx->Bp[0],&pcctx->Bp[1]);
1092:   KSPGetPC(ksp,&pcctx->pc);
1093:   PetscObjectReference((PetscObject)pcctx->pc);
1094:   MatGetLocalSize(pep->A[0],&mloc,&nloc);
1095:   if (pjd->ld>1) {
1096:     nloc += pjd->ld-1; mloc += pjd->ld-1;
1097:   }
1098:   PetscNew(&matctx);
1099:   MatCreateShell(PetscObjectComm((PetscObject)pep),kspsf*nloc,kspsf*mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&pjd->Pshell);
1100:   MatShellSetOperation(pjd->Pshell,MATOP_MULT,(void(*)(void))PEPJDShellMatMult);
1101:   MatShellSetOperation(pjd->Pshell,MATOP_CREATE_VECS,(void(*)(void))PEPJDSellMatCreateVecs);
1102:   matctx->pep = pep;
1103:   target[0] = pep->target; target[1] = 0.0;
1104:   PEPJDMatSetUp(pep,1,target);
1105:   Pr = matctx->Pr;
1106:   pcctx->PPr = NULL;
1107: #if !defined(PETSC_USE_COMPLEX)
1108:   if (!pjd->reusepc) {
1109:     MatDuplicate(matctx->Pr,MAT_COPY_VALUES,&pcctx->PPr);
1110:     Pr = pcctx->PPr;
1111:   }
1112: #endif
1113:   PCSetOperators(pcctx->pc,Pr,Pr);
1114:   PCSetErrorIfFailure(pcctx->pc,PETSC_TRUE);
1115:   KSPSetPC(ksp,pjd->pcshell);
1116:   if (pjd->reusepc) {
1117:     PCSetReusePreconditioner(pcctx->pc,PETSC_TRUE);
1118:     KSPSetReusePreconditioner(ksp,PETSC_TRUE);
1119:   }
1120:   KSPSetOperators(ksp,pjd->Pshell,pjd->Pshell);
1121:   KSPSetUp(ksp);
1122:   if (pjd->ld>1) {
1123:     PetscMalloc2(pjd->ld*pjd->ld,&pcctx->M,pjd->ld*pjd->ld,&pcctx->ps);
1124:     pcctx->pep = pep;
1125:   }
1126:   matctx->work = ww;
1127:   pcctx->work  = ww;
1128:   return(0);
1129: }

1131: static PetscErrorCode PEPJDEigenvectors(PEP pep)
1132: {
1133: #if defined(SLEPC_MISSING_LAPACK_TREVC)
1135:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable");
1136: #else
1138:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1139:   PetscBLASInt   ld,nconv,info,nc;
1140:   PetscScalar    *Z,*w;
1141:   PetscReal      *wr,norm;
1142:   PetscInt       i;
1143:   Mat            U;
1144: #if !defined(PETSC_USE_COMPLEX)
1145:   Vec            v,v1;
1146: #endif

1149:   PetscMalloc3(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr,2*pep->ncv,&w);
1150:   PetscBLASIntCast(pep->ncv,&ld);
1151:   PetscBLASIntCast(pep->nconv,&nconv);
1152: #if !defined(PETSC_USE_COMPLEX)
1153:   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,wr,&info));
1154: #else
1155:   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,w,wr,&info));
1156: #endif
1157:   SlepcCheckLapackInfo("trevc",info);
1158:   MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,Z,&U);
1159:   BVSetActiveColumns(pjd->X,0,pep->nconv);
1160:   BVMultInPlace(pjd->X,U,0,pep->nconv);
1161:   for (i=0;i<pep->nconv;i++) {
1162: #if !defined(PETSC_USE_COMPLEX)
1163:     if (pep->eigi[i]!=0.0) {   /* first eigenvalue of a complex conjugate pair */
1164:       BVGetColumn(pjd->X,i,&v);
1165:       BVGetColumn(pjd->X,i+1,&v1);
1166:       VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
1167:       BVRestoreColumn(pjd->X,i,&v);
1168:       BVRestoreColumn(pjd->X,i+1,&v1);
1169:       i++;
1170:     } else   /* real eigenvalue */
1171: #endif
1172:     {
1173:       BVNormColumn(pjd->X,i,NORM_2,&norm);
1174:       BVScaleColumn(pjd->X,i,1.0/norm);
1175:     }
1176:   }
1177:   MatDestroy(&U);
1178:   PetscFree3(Z,wr,w);
1179:   return(0);
1180: #endif
1181: }

1183: static PetscErrorCode PEPJDLockConverged(PEP pep,PetscInt *nv,PetscInt sz)
1184: {
1186:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1187:   PetscInt       j,i,*P,ldds,rk=0,nvv=*nv;
1188:   Vec            v,x,w;
1189:   PetscScalar    *R,*r,*pX,target[2];
1190:   Mat            X;
1191:   PetscBLASInt   sz_,rk_,nv_,info;
1192:   PetscMPIInt    np;

1195:   /* update AX and XpX */
1196:   for (i=sz;i>0;i--) {
1197:     BVGetColumn(pjd->X,pjd->nlock-i,&x);
1198:     for (j=0;j<pep->nmat;j++) {
1199:       BVGetColumn(pjd->AX[j],pjd->nlock-i,&v);
1200:       MatMult(pep->A[j],x,v);
1201:       BVRestoreColumn(pjd->AX[j],pjd->nlock-i,&v);
1202:       BVSetActiveColumns(pjd->AX[j],0,pjd->nlock-i+1);
1203:     }
1204:     BVRestoreColumn(pjd->X,pjd->nlock-i,&x);
1205:     BVDotColumn(pjd->X,(pjd->nlock-i),pjd->XpX+(pjd->nlock-i)*(pjd->ld));
1206:     pjd->XpX[(pjd->nlock-i)*(1+pjd->ld)] = 1.0;
1207:     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]);
1208:   }

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

1213:   /* evaluate the polynomial basis in T */
1214:   PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat);
1215:   for (j=0;j<pep->nmat;j++) {
1216:     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);
1217:   }

1219:   /* Extend search space */
1220:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
1221:   PetscCalloc3(nvv,&P,nvv*nvv,&R,nvv*sz,&r);
1222:   DSGetLeadingDimension(pep->ds,&ldds);
1223:   DSGetArray(pep->ds,DS_MAT_X,&pX);
1224:   PEPJDOrthogonalize(nvv,nvv,pX,ldds,&rk,P,R,nvv);
1225:   for (j=0;j<sz;j++) {
1226:     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 */
1227:   }
1228:   PetscBLASIntCast(rk,&rk_);
1229:   PetscBLASIntCast(sz,&sz_);
1230:   PetscBLASIntCast(nvv,&nv_);
1231:   PetscStackCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&rk_,R,&nv_,&info));
1232:   SlepcCheckLapackInfo("trtri",info);
1233:   for (i=0;i<sz;i++) PetscStackCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&rk_,R,&nv_,r+i,&sz_));
1234:   for (i=0;i<sz*rk;i++) r[i] = PetscConj(r[i])/PetscSqrtReal(np); /* revert */
1235:   BVSetActiveColumns(pjd->V,0,nvv);
1236:   rk -= sz;
1237:   for (j=0;j<rk;j++) {
1238:     PetscArraycpy(R+j*nvv,pX+(j+sz)*ldds,nvv);
1239:   }
1240:   DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1241:   MatCreateSeqDense(PETSC_COMM_SELF,nvv,rk,R,&X);
1242:   BVMultInPlace(pjd->V,X,0,rk);
1243:   MatDestroy(&X);
1244:   BVSetActiveColumns(pjd->V,0,rk);
1245:   for (j=0;j<rk;j++) {
1246:     BVGetColumn(pjd->V,j,&v);
1247:     PEPJDCopyToExtendedVec(pep,NULL,r+sz*(j+sz),sz,pjd->nlock-sz,v,PETSC_FALSE);
1248:     BVRestoreColumn(pjd->V,j,&v);
1249:   }
1250:   BVOrthogonalize(pjd->V,NULL);

1252:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1253:     for (j=0;j<rk;j++) {
1254:       /* W = P(target)*V */
1255:       BVGetColumn(pjd->W,j,&w);
1256:       BVGetColumn(pjd->V,j,&v);
1257:       target[0] = pep->target; target[1] = 0.0;
1258:       PEPJDComputeResidual(pep,PETSC_FALSE,1,&v,target,&w,pep->work);
1259:       BVRestoreColumn(pjd->V,j,&v);
1260:       BVRestoreColumn(pjd->W,j,&w);
1261:     }
1262:     BVSetActiveColumns(pjd->W,0,rk);
1263:     BVOrthogonalize(pjd->W,NULL);
1264:   }
1265:   *nv = rk;
1266:   PetscFree3(P,R,r);
1267:   return(0);
1268: }

1270: PetscErrorCode PEPJDSystemSetUp(PEP pep,PetscInt sz,PetscScalar *theta,Vec *u,Vec *p,Vec *ww)
1271: {
1273:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1274:   PEP_JD_PCSHELL *pcctx;
1275: #if !defined(PETSC_USE_COMPLEX)
1276:   PetscScalar    s[2];
1277: #endif

1280:   PCShellGetContext(pjd->pcshell,(void**)&pcctx);
1281:   PEPJDMatSetUp(pep,sz,theta);
1282:   pcctx->u[0] = u[0]; pcctx->u[1] = u[1];
1283:   /* Compute r'. p is a work space vector */
1284:   PEPJDComputeResidual(pep,PETSC_TRUE,sz,u,theta,p,ww);
1285:   PEPJDExtendedPCApply(pjd->pcshell,p[0],pcctx->Bp[0]);
1286:   VecDot(pcctx->Bp[0],u[0],pcctx->gamma);
1287: #if !defined(PETSC_USE_COMPLEX)
1288:   if (sz==2) {
1289:     PEPJDExtendedPCApply(pjd->pcshell,p[1],pcctx->Bp[1]);
1290:     VecDot(pcctx->Bp[0],u[1],pcctx->gamma+1);
1291:     VecMDot(pcctx->Bp[1],2,u,s);
1292:     pcctx->gamma[0] += s[1];
1293:     pcctx->gamma[1] = -pcctx->gamma[1]+s[0];
1294:   }
1295: #endif
1296:   if (sz==1) {
1297:     VecZeroEntries(pcctx->Bp[1]);
1298:     pcctx->gamma[1] = 0.0;
1299:   }
1300:   return(0);
1301: }

1303: PetscErrorCode PEPSolve_JD(PEP pep)
1304: {
1305:   PetscErrorCode  ierr;
1306:   PEP_JD          *pjd = (PEP_JD*)pep->data;
1307:   PetscInt        k,nv,nvc,ld,minv,dim,bupdated=0,sz=1,kspsf=1,idx,off,maxits,nloc;
1308:   PetscMPIInt     np,count;
1309:   PetscScalar     theta[2]={0.0,0.0},ritz[2]={0.0,0.0},*pX,*eig,*eigi,*array;
1310:   PetscReal       norm,*res,tol=0.0,rtol,abstol, dtol;
1311:   PetscBool       lindep,ini=PETSC_TRUE;
1312:   Vec             tc,t[2]={NULL,NULL},u[2]={NULL,NULL},p[2]={NULL,NULL};
1313:   Vec             rc,rr[2],r[2]={NULL,NULL},*ww=pep->work,v[2];
1314:   Mat             G,X,Y;
1315:   KSP             ksp;
1316:   PEP_JD_PCSHELL  *pcctx;
1317:   PEP_JD_MATSHELL *matctx;
1318: #if !defined(PETSC_USE_COMPLEX)
1319:   PetscReal       norm1;
1320: #endif

1323:   PetscCitationsRegister(citation,&cited);
1324:   MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np);
1325:   BVGetSizes(pep->V,&nloc,NULL,NULL);
1326:   DSGetLeadingDimension(pep->ds,&ld);
1327:   PetscCalloc3(pep->ncv+pep->nev,&eig,pep->ncv+pep->nev,&eigi,pep->ncv+pep->nev,&res);
1328:   pjd->nlock = 0;
1329:   STGetKSP(pep->st,&ksp);
1330:   KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits);
1331: #if !defined (PETSC_USE_COMPLEX)
1332:   kspsf = 2;
1333: #endif
1334:   PEPJDProcessInitialSpace(pep,ww);
1335:   nv = (pep->nini)?pep->nini:1;

1337:   /* Replace preconditioner with one containing projectors */
1338:   PEPJDCreateShellPC(pep,ww);
1339:   PCShellGetContext(pjd->pcshell,(void**)&pcctx);

1341:   /* Create auxiliar vectors */
1342:   BVCreateVec(pjd->V,&u[0]);
1343:   VecDuplicate(u[0],&p[0]);
1344:   VecDuplicate(u[0],&r[0]);
1345: #if !defined (PETSC_USE_COMPLEX)
1346:   VecDuplicate(u[0],&u[1]);
1347:   VecDuplicate(u[0],&p[1]);
1348:   VecDuplicate(u[0],&r[1]);
1349: #endif

1351:   /* Restart loop */
1352:   while (pep->reason == PEP_CONVERGED_ITERATING) {
1353:     pep->its++;
1354:     DSSetDimensions(pep->ds,nv,0,0,0);
1355:     BVSetActiveColumns(pjd->V,bupdated,nv);
1356:     PEPJDUpdateTV(pep,bupdated,nv,ww);
1357:     if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) { BVSetActiveColumns(pjd->W,bupdated,nv); }
1358:     for (k=0;k<pep->nmat;k++) {
1359:       BVSetActiveColumns(pjd->TV[k],bupdated,nv);
1360:       DSGetMat(pep->ds,DSMatExtra[k],&G);
1361:       BVMatProject(pjd->TV[k],NULL,pjd->W,G);
1362:       DSRestoreMat(pep->ds,DSMatExtra[k],&G);
1363:     }
1364:     BVSetActiveColumns(pjd->V,0,nv);
1365:     BVSetActiveColumns(pjd->W,0,nv);

1367:     /* Solve projected problem */
1368:     DSSetState(pep->ds,DS_STATE_RAW);
1369:     DSSolve(pep->ds,pep->eigr,pep->eigi);
1370:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
1371:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);
1372:     idx = 0;
1373:     do {
1374:       ritz[0] = pep->eigr[idx];
1375: #if !defined(PETSC_USE_COMPLEX)
1376:       ritz[1] = pep->eigi[idx];
1377:       sz = (ritz[1]==0.0)?1:2;
1378: #endif
1379:       /* Compute Ritz vector u=V*X(:,1) */
1380:       DSGetArray(pep->ds,DS_MAT_X,&pX);
1381:       BVSetActiveColumns(pjd->V,0,nv);
1382:       BVMultVec(pjd->V,1.0,0.0,u[0],pX+idx*ld);
1383: #if !defined(PETSC_USE_COMPLEX)
1384:       if (sz==2) {
1385:         BVMultVec(pjd->V,1.0,0.0,u[1],pX+(idx+1)*ld);
1386:       }
1387: #endif
1388:       DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1389:       PEPJDComputeResidual(pep,PETSC_FALSE,sz,u,ritz,r,ww);
1390:       /* Check convergence */
1391:       VecNorm(r[0],NORM_2,&norm);
1392: #if !defined(PETSC_USE_COMPLEX)
1393:       if (sz==2) {
1394:         VecNorm(r[1],NORM_2,&norm1);
1395:         norm = SlepcAbs(norm,norm1);
1396:       }
1397: #endif
1398:       (*pep->converged)(pep,ritz[0],ritz[1],norm,&pep->errest[pep->nconv],pep->convergedctx);
1399:       if (sz==2) pep->errest[pep->nconv+1] = pep->errest[pep->nconv];
1400:       if (ini) {
1401:         tol = PetscMin(.1,pep->errest[pep->nconv]); ini = PETSC_FALSE;
1402:       } else tol = PetscMin(pep->errest[pep->nconv],tol/2);
1403:       (*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);
1404:       if (pep->errest[pep->nconv]<pep->tol) {
1405:         /* Ritz pair converged */
1406:         ini = PETSC_TRUE;
1407:         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1408:         if (pjd->ld>1) {
1409:           BVGetColumn(pjd->X,pep->nconv,&v[0]);
1410:           PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*pep->nconv,pjd->ld-1,0,u[0],PETSC_TRUE);
1411:           BVRestoreColumn(pjd->X,pep->nconv,&v[0]);
1412:           BVSetActiveColumns(pjd->X,0,pep->nconv+1);
1413:           BVNormColumn(pjd->X,pep->nconv,NORM_2,&norm);
1414:           BVScaleColumn(pjd->X,pep->nconv,1.0/norm);
1415:           for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*pep->nconv+k] *= PetscSqrtReal(np)/norm;
1416:           pjd->T[(pep->ncv+1)*pep->nconv] = ritz[0];
1417:           eig[pep->nconv] = ritz[0];
1418:           idx++;
1419: #if !defined(PETSC_USE_COMPLEX)
1420:           if (sz==2) {
1421:             BVGetColumn(pjd->X,pep->nconv+1,&v[0]);
1422:             PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*(pep->nconv+1),pjd->ld-1,0,u[1],PETSC_TRUE);
1423:             BVRestoreColumn(pjd->X,pep->nconv+1,&v[0]);
1424:             BVSetActiveColumns(pjd->X,0,pep->nconv+2);
1425:             BVNormColumn(pjd->X,pep->nconv+1,NORM_2,&norm1);
1426:             BVScaleColumn(pjd->X,pep->nconv+1,1.0/norm1);
1427:             for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*(pep->nconv+1)+k] *= PetscSqrtReal(np)/norm1;
1428:             pjd->T[(pep->ncv+1)*(pep->nconv+1)] = ritz[0];
1429:             pjd->T[(pep->ncv+1)*pep->nconv+1] = -ritz[1]*norm1/norm;
1430:             pjd->T[(pep->ncv+1)*(pep->nconv+1)-1] = ritz[1]*norm/norm1;
1431:             eig[pep->nconv+1] = ritz[0];
1432:             eigi[pep->nconv] = ritz[1]; eigi[pep->nconv+1] = -ritz[1];
1433:             idx++;
1434:           }
1435: #endif
1436:         } else {
1437:           BVInsertVec(pep->V,pep->nconv,u[0]);
1438:         }
1439:         pep->nconv += sz;
1440:       }
1441:     } while (pep->errest[pep->nconv]<pep->tol && pep->nconv<nv);

1443:     if (pep->reason==PEP_CONVERGED_ITERATING) {
1444:       nvc = nv;
1445:       if (idx) {
1446:         pjd->nlock +=idx;
1447:         PEPJDLockConverged(pep,&nv,idx);
1448:       }
1449:       if (nv+sz>=pep->ncv-1) {
1450:         /* Basis full, force restart */
1451:         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1452:         DSGetDimensions(pep->ds,&dim,NULL,NULL,NULL,NULL);
1453:         DSGetArray(pep->ds,DS_MAT_X,&pX);
1454:         PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld);
1455:         DSRestoreArray(pep->ds,DS_MAT_X,&pX);
1456:         DSGetArray(pep->ds,DS_MAT_Y,&pX);
1457:         PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld);
1458:         DSRestoreArray(pep->ds,DS_MAT_Y,&pX);
1459:         DSGetMat(pep->ds,DS_MAT_X,&X);
1460:         BVMultInPlace(pjd->V,X,0,minv);
1461:         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1462:          DSGetMat(pep->ds,DS_MAT_Y,&Y);
1463:          BVMultInPlace(pjd->W,Y,pep->nconv,minv);
1464:          DSRestoreMat(pep->ds,DS_MAT_Y,&Y);
1465:         }
1466:         MatDestroy(&X);
1467:         nv = minv;
1468:         bupdated = 0;
1469:       } else {
1470:         if (!idx && pep->errest[pep->nconv]<pjd->fix) {theta[0] = ritz[0]; theta[1] = ritz[1];}
1471:         else {theta[0] = pep->target; theta[1] = 0.0;}
1472:         /* Update system mat */
1473:         PEPJDSystemSetUp(pep,sz,theta,u,p,ww);
1474:         /* Solve correction equation to expand basis */
1475:         BVGetColumn(pjd->V,nv,&t[0]);
1476:         rr[0] = r[0];
1477:         if (sz==2) {
1478:           BVGetColumn(pjd->V,nv+1,&t[1]);
1479:           rr[1] = r[1];
1480:         } else {
1481:           t[1] = NULL;
1482:           rr[1] = NULL;
1483:         }
1484:         VecCreateCompWithVecs(t,kspsf,pjd->vtempl,&tc);
1485:         VecCreateCompWithVecs(rr,kspsf,pjd->vtempl,&rc);
1486:         VecCompSetSubVecs(pjd->vtempl,sz,NULL);
1487:         tol  = PetscMax(rtol,tol/2);
1488:         KSPSetTolerances(ksp,tol,abstol,dtol,maxits);
1489:         KSPSolve(ksp,rc,tc);
1490:         VecDestroy(&tc);
1491:         VecDestroy(&rc);
1492:         VecGetArray(t[0],&array);
1493:         PetscMPIIntCast(pep->nconv,&count);
1494:         MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
1495:         VecRestoreArray(t[0],&array);
1496:         BVRestoreColumn(pjd->V,nv,&t[0]);
1497:         BVOrthogonalizeColumn(pjd->V,nv,NULL,&norm,&lindep);
1498:         if (lindep || norm==0.0) {
1499:           if (sz==1) SETERRQ(PETSC_COMM_SELF,1,"Linearly dependent continuation vector");
1500:           else off = 1;
1501:         } else {
1502:           off = 0;
1503:           BVScaleColumn(pjd->V,nv,1.0/norm);
1504:         }
1505: #if !defined(PETSC_USE_COMPLEX)
1506:         if (sz==2) {
1507:           VecGetArray(t[1],&array);
1508:           MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep));
1509:           VecRestoreArray(t[1],&array);
1510:           BVRestoreColumn(pjd->V,nv+1,&t[1]);
1511:           if (off) {
1512:             BVCopyColumn(pjd->V,nv+1,nv);
1513:           }
1514:           BVOrthogonalizeColumn(pjd->V,nv+1-off,NULL,&norm,&lindep);
1515:           if (lindep || norm==0.0) {
1516:             if (off) SETERRQ(PETSC_COMM_SELF,1,"Linearly dependent continuation vector");
1517:             else off = 1;
1518:           } else {
1519:             BVScaleColumn(pjd->V,nv+1-off,1.0/norm);
1520:           }
1521:         }
1522: #endif
1523:         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1524:           BVInsertVec(pjd->W,nv,r[0]);
1525:           if (sz==2 && !off) {
1526:             BVInsertVec(pjd->W,nv+1,r[1]);
1527:           }
1528:           BVOrthogonalizeColumn(pjd->W,nv,NULL,&norm,&lindep);
1529:           if (lindep || norm==0.0) SETERRQ(PETSC_COMM_SELF,1,"Linearly dependent continuation vector");
1530:           BVScaleColumn(pjd->W,nv,1.0/norm);
1531:           if (sz==2 && !off) {
1532:             BVOrthogonalizeColumn(pjd->W,nv+1,NULL,&norm,&lindep);
1533:             if (lindep || norm==0.0) SETERRQ(PETSC_COMM_SELF,1,"Linearly dependent continuation vector");
1534:             BVScaleColumn(pjd->W,nv+1,1.0/norm);
1535:           }
1536:         }
1537:         bupdated = idx?0:nv;
1538:         nv += sz-off;
1539:       }
1540:       for (k=0;k<nvc;k++) {
1541:         eig[pep->nconv-idx+k] = pep->eigr[k];
1542: #if !defined(PETSC_USE_COMPLEX)
1543:         eigi[pep->nconv-idx+k] = pep->eigi[k];
1544: #endif
1545:       }
1546:       PEPMonitor(pep,pep->its,pep->nconv,eig,eigi,pep->errest,pep->nconv+1);
1547:     }
1548:   }
1549:   if (pjd->ld>1) {
1550:     for (k=0;k<pep->nconv;k++) {
1551:       pep->eigr[k] = eig[k];
1552:       pep->eigi[k] = eigi[k];
1553:     }
1554:     if (pep->nconv>0) { PEPJDEigenvectors(pep); }
1555:     PetscFree2(pcctx->M,pcctx->ps);
1556:   }
1557:   VecDestroy(&u[0]);
1558:   VecDestroy(&r[0]);
1559:   VecDestroy(&p[0]);
1560: #if !defined (PETSC_USE_COMPLEX)
1561:   VecDestroy(&u[1]);
1562:   VecDestroy(&r[1]);
1563:   VecDestroy(&p[1]);
1564: #endif
1565:   KSPSetTolerances(ksp,rtol,abstol,dtol,maxits);
1566:   KSPSetPC(ksp,pcctx->pc);
1567:   VecDestroy(&pcctx->Bp[0]);
1568:   VecDestroy(&pcctx->Bp[1]);
1569:   MatShellGetContext(pjd->Pshell,(void**)&matctx);
1570:   MatDestroy(&matctx->Pr);
1571:   MatDestroy(&matctx->Pi);
1572:   MatDestroy(&pjd->Pshell);
1573:   MatDestroy(&pcctx->PPr);
1574:   PCDestroy(&pcctx->pc);
1575:   PetscFree(pcctx);
1576:   PetscFree(matctx);
1577:   PCDestroy(&pjd->pcshell);
1578:   PetscFree3(eig,eigi,res);
1579:   VecDestroy(&pjd->vtempl);
1580:   return(0);
1581: }

1583: PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
1584: {
1585:   PEP_JD *pjd = (PEP_JD*)pep->data;

1588:   if (keep==PETSC_DEFAULT) pjd->keep = 0.5;
1589:   else {
1590:     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]");
1591:     pjd->keep = keep;
1592:   }
1593:   return(0);
1594: }

1596: /*@
1597:    PEPJDSetRestart - Sets the restart parameter for the Jacobi-Davidson
1598:    method, in particular the proportion of basis vectors that must be kept
1599:    after restart.

1601:    Logically Collective on pep

1603:    Input Parameters:
1604: +  pep  - the eigenproblem solver context
1605: -  keep - the number of vectors to be kept at restart

1607:    Options Database Key:
1608: .  -pep_jd_restart - Sets the restart parameter

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

1613:    Level: advanced

1615: .seealso: PEPJDGetRestart()
1616: @*/
1617: PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
1618: {

1624:   PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
1625:   return(0);
1626: }

1628: PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
1629: {
1630:   PEP_JD *pjd = (PEP_JD*)pep->data;

1633:   *keep = pjd->keep;
1634:   return(0);
1635: }

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

1640:    Not Collective

1642:    Input Parameter:
1643: .  pep - the eigenproblem solver context

1645:    Output Parameter:
1646: .  keep - the restart parameter

1648:    Level: advanced

1650: .seealso: PEPJDSetRestart()
1651: @*/
1652: PetscErrorCode PEPJDGetRestart(PEP pep,PetscReal *keep)
1653: {

1659:   PetscUseMethod(pep,"PEPJDGetRestart_C",(PEP,PetscReal*),(pep,keep));
1660:   return(0);
1661: }

1663: PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
1664: {
1665:   PEP_JD *pjd = (PEP_JD*)pep->data;

1668:   if (fix == PETSC_DEFAULT || fix == PETSC_DECIDE) pjd->fix = 0.01;
1669:   else {
1670:     if (fix < 0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid fix value");
1671:     pjd->fix = fix;
1672:   }
1673:   return(0);
1674: }

1676: /*@
1677:    PEPJDSetFix - Sets the threshold for changing the target in the correction
1678:    equation.

1680:    Logically Collective on pep

1682:    Input Parameters:
1683: +  pep - the eigenproblem solver context
1684: -  fix - threshold for changing the target

1686:    Options Database Key:
1687: .  -pep_jd_fix - the fix value

1689:    Note:
1690:    The target in the correction equation is fixed at the first iterations.
1691:    When the norm of the residual vector is lower than the fix value,
1692:    the target is set to the corresponding eigenvalue.

1694:    Level: advanced

1696: .seealso: PEPJDGetFix()
1697: @*/
1698: PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
1699: {

1705:   PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
1706:   return(0);
1707: }

1709: PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
1710: {
1711:   PEP_JD *pjd = (PEP_JD*)pep->data;

1714:   *fix = pjd->fix;
1715:   return(0);
1716: }

1718: /*@
1719:    PEPJDGetFix - Returns the threshold for changing the target in the correction
1720:    equation.

1722:    Not Collective

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

1727:    Output Parameter:
1728: .  fix - threshold for changing the target

1730:    Note:
1731:    The target in the correction equation is fixed at the first iterations.
1732:    When the norm of the residual vector is lower than the fix value,
1733:    the target is set to the corresponding eigenvalue.

1735:    Level: advanced

1737: .seealso: PEPJDSetFix()
1738: @*/
1739: PetscErrorCode PEPJDGetFix(PEP pep,PetscReal *fix)
1740: {

1746:   PetscUseMethod(pep,"PEPJDGetFix_C",(PEP,PetscReal*),(pep,fix));
1747:   return(0);
1748: }

1750: PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
1751: {
1752:   PEP_JD *pjd = (PEP_JD*)pep->data;

1755:   pjd->reusepc = reusepc;
1756:   return(0);
1757: }

1759: /*@
1760:    PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
1761:    must be reused or not.

1763:    Logically Collective on pep

1765:    Input Parameters:
1766: +  pep     - the eigenproblem solver context
1767: -  reusepc - the reuse flag

1769:    Options Database Key:
1770: .  -pep_jd_reuse_preconditioner - the reuse flag

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

1777:    Level: advanced

1779: .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
1780: @*/
1781: PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
1782: {

1788:   PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
1789:   return(0);
1790: }

1792: PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
1793: {
1794:   PEP_JD *pjd = (PEP_JD*)pep->data;

1797:   *reusepc = pjd->reusepc;
1798:   return(0);
1799: }

1801: /*@
1802:    PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.

1804:    Not Collective

1806:    Input Parameter:
1807: .  pep - the eigenproblem solver context

1809:    Output Parameter:
1810: .  reusepc - the reuse flag

1812:    Level: advanced

1814: .seealso: PEPJDSetReusePreconditioner()
1815: @*/
1816: PetscErrorCode PEPJDGetReusePreconditioner(PEP pep,PetscBool *reusepc)
1817: {

1823:   PetscUseMethod(pep,"PEPJDGetReusePreconditioner_C",(PEP,PetscBool*),(pep,reusepc));
1824:   return(0);
1825: }

1827: PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
1828: {
1829:   PEP_JD *pjd = (PEP_JD*)pep->data;

1832:   if (mmidx == PETSC_DEFAULT || mmidx == PETSC_DECIDE) pjd->mmidx = 1;
1833:   else {
1834:     if (mmidx < 1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mmidx value");
1835:     pjd->mmidx = mmidx;
1836:     pep->state = PEP_STATE_INITIAL;
1837:   }
1838:   return(0);
1839: }

1841: /*@
1842:    PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.

1844:    Logically Collective on pep

1846:    Input Parameters:
1847: +  pep   - the eigenproblem solver context
1848: -  mmidx - maximum minimality index

1850:    Options Database Key:
1851: .  -pep_jd_minimality_index - the minimality index value

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

1857:    Level: advanced

1859: .seealso: PEPJDGetMinimalityIndex()
1860: @*/
1861: PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
1862: {

1868:   PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
1869:   return(0);
1870: }

1872: PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
1873: {
1874:   PEP_JD *pjd = (PEP_JD*)pep->data;

1877:   *mmidx = pjd->mmidx;
1878:   return(0);
1879: }

1881: /*@
1882:    PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
1883:    index.

1885:    Not Collective

1887:    Input Parameter:
1888: .  pep - the eigenproblem solver context

1890:    Output Parameter:
1891: .  mmidx - minimality index

1893:    Level: advanced

1895: .seealso: PEPJDSetMinimalityIndex()
1896: @*/
1897: PetscErrorCode PEPJDGetMinimalityIndex(PEP pep,PetscInt *mmidx)
1898: {

1904:   PetscUseMethod(pep,"PEPJDGetMinimalityIndex_C",(PEP,PetscInt*),(pep,mmidx));
1905:   return(0);
1906: }

1908: PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
1909: {
1910:   PEP_JD *pjd = (PEP_JD*)pep->data;

1913:   switch (proj) {
1914:     case PEP_JD_PROJECTION_HARMONIC:
1915:     case PEP_JD_PROJECTION_ORTHOGONAL:
1916:       if (pjd->proj != proj) {
1917:         pep->state = PEP_STATE_INITIAL;
1918:         pjd->proj = proj;
1919:       }
1920:       break;
1921:     default:
1922:       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'proj' value");
1923:   }
1924:   return(0);
1925: }

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

1930:    Logically Collective on pep

1932:    Input Parameters:
1933: +  pep  - the eigenproblem solver context
1934: -  proj - the type of projection

1936:    Options Database Key:
1937: .  -pep_jd_projection - the projection type, either orthogonal or harmonic

1939:    Level: advanced

1941: .seealso: PEPJDGetProjection()
1942: @*/
1943: PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
1944: {

1950:   PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
1951:   return(0);
1952: }

1954: PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
1955: {
1956:   PEP_JD *pjd = (PEP_JD*)pep->data;

1959:   *proj = pjd->proj;
1960:   return(0);
1961: }

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

1966:    Not Collective

1968:    Input Parameter:
1969: .  pep - the eigenproblem solver context

1971:    Output Parameter:
1972: .  proj - the type of projection

1974:    Level: advanced

1976: .seealso: PEPJDSetProjection()
1977: @*/
1978: PetscErrorCode PEPJDGetProjection(PEP pep,PEPJDProjection *proj)
1979: {

1985:   PetscUseMethod(pep,"PEPJDGetProjection_C",(PEP,PEPJDProjection*),(pep,proj));
1986:   return(0);
1987: }

1989: PetscErrorCode PEPSetFromOptions_JD(PetscOptionItems *PetscOptionsObject,PEP pep)
1990: {
1991:   PetscErrorCode  ierr;
1992:   PetscBool       flg,b1;
1993:   PetscReal       r1;
1994:   PetscInt        i1;
1995:   PEPJDProjection proj;

1998:   PetscOptionsHead(PetscOptionsObject,"PEP JD Options");

2000:     PetscOptionsReal("-pep_jd_restart","Proportion of vectors kept after restart","PEPJDSetRestart",0.5,&r1,&flg);
2001:     if (flg) { PEPJDSetRestart(pep,r1); }

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

2006:     PetscOptionsBool("-pep_jd_reuse_preconditioner","Whether to reuse the preconditioner","PEPJDSetReusePreconditoiner",PETSC_FALSE,&b1,&flg);
2007:     if (flg) { PEPJDSetReusePreconditioner(pep,b1); }

2009:     PetscOptionsInt("-pep_jd_minimality_index","Maximum allowed minimality index","PEPJDSetMinimalityIndex",1,&i1,&flg);
2010:     if (flg) { PEPJDSetMinimalityIndex(pep,i1); }

2012:     PetscOptionsEnum("-pep_jd_projection","Type of projection","PEPJDSetProjection",PEPJDProjectionTypes,(PetscEnum)PEP_JD_PROJECTION_HARMONIC,(PetscEnum*)&proj,&flg);
2013:     if (flg) { PEPJDSetProjection(pep,proj); }

2015:   PetscOptionsTail();
2016:   return(0);
2017: }

2019: PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
2020: {
2022:   PEP_JD         *pjd = (PEP_JD*)pep->data;
2023:   PetscBool      isascii;

2026:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2027:   if (isascii) {
2028:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*pjd->keep));
2029:     PetscViewerASCIIPrintf(viewer,"  threshold for changing the target in the correction equation (fix): %g\n",(double)pjd->fix);
2030:     PetscViewerASCIIPrintf(viewer,"  projection type: %s\n",PEPJDProjectionTypes[pjd->proj]);
2031:     PetscViewerASCIIPrintf(viewer,"  maximum allowed minimality index: %d\n",pjd->mmidx);
2032:     if (pjd->reusepc) { PetscViewerASCIIPrintf(viewer,"  reusing the preconditioner\n"); }
2033:   }
2034:   return(0);
2035: }

2037: PetscErrorCode PEPSetDefaultST_JD(PEP pep)
2038: {
2040:   KSP            ksp;

2043:   if (!((PetscObject)pep->st)->type_name) {
2044:     STSetType(pep->st,STPRECOND);
2045:     STPrecondSetKSPHasMat(pep->st,PETSC_TRUE);
2046:   }
2047:   STSetTransform(pep->st,PETSC_FALSE);
2048:   STGetKSP(pep->st,&ksp);
2049:   if (!((PetscObject)ksp)->type_name) {
2050:     KSPSetType(ksp,KSPBCGSL);
2051:     KSPSetTolerances(ksp,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,100);
2052:   }
2053:   return(0);
2054: }

2056: PetscErrorCode PEPReset_JD(PEP pep)
2057: {
2059:   PEP_JD         *pjd = (PEP_JD*)pep->data;
2060:   PetscInt       i;

2063:   for (i=0;i<pep->nmat;i++) {
2064:     BVDestroy(pjd->TV+i);
2065:   }
2066:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) { BVDestroy(&pjd->W); }
2067:   if (pjd->ld>1) {
2068:     BVDestroy(&pjd->V);
2069:     for (i=0;i<pep->nmat;i++) {
2070:       BVDestroy(pjd->AX+i);
2071:     }
2072:     BVDestroy(&pjd->N[0]);
2073:     BVDestroy(&pjd->N[1]);
2074:     PetscFree3(pjd->XpX,pjd->T,pjd->Tj);
2075:   }
2076:   PetscFree2(pjd->TV,pjd->AX);
2077:   return(0);
2078: }

2080: PetscErrorCode PEPDestroy_JD(PEP pep)
2081: {

2085:   PetscFree(pep->data);
2086:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",NULL);
2087:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",NULL);
2088:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",NULL);
2089:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",NULL);
2090:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",NULL);
2091:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",NULL);
2092:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",NULL);
2093:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",NULL);
2094:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",NULL);
2095:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",NULL);
2096:   return(0);
2097: }

2099: SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
2100: {
2101:   PEP_JD         *pjd;

2105:   PetscNewLog(pep,&pjd);
2106:   pep->data = (void*)pjd;

2108:   pjd->fix   = 0.01;
2109:   pjd->mmidx = 0;

2111:   pep->ops->solve          = PEPSolve_JD;
2112:   pep->ops->setup          = PEPSetUp_JD;
2113:   pep->ops->setfromoptions = PEPSetFromOptions_JD;
2114:   pep->ops->destroy        = PEPDestroy_JD;
2115:   pep->ops->reset          = PEPReset_JD;
2116:   pep->ops->view           = PEPView_JD;
2117:   pep->ops->setdefaultst   = PEPSetDefaultST_JD;

2119:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",PEPJDSetRestart_JD);
2120:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",PEPJDGetRestart_JD);
2121:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",PEPJDSetFix_JD);
2122:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",PEPJDGetFix_JD);
2123:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",PEPJDSetReusePreconditioner_JD);
2124:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",PEPJDGetReusePreconditioner_JD);
2125:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",PEPJDSetMinimalityIndex_JD);
2126:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",PEPJDGetMinimalityIndex_JD);
2127:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",PEPJDSetProjection_JD);
2128:   PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",PEPJDGetProjection_JD);
2129:   return(0);
2130: }

Coincidencia en el fichero binario (entrada estándar)