Actual source code: pjd.c

slepc-3.20.1 2023-11-27
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:   PetscFunctionBegin;
110:   if (pjd->ld>1) {
111:     PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),basis));
112:     PetscCall(BVGetSizes(pep->V,&nloc,NULL,&m));
113:     nloc += pjd->ld-1;
114:     PetscCall(BVSetSizes(*basis,nloc,PETSC_DECIDE,m));
115:     PetscCall(BVGetType(pep->V,&type));
116:     PetscCall(BVSetType(*basis,type));
117:     PetscCall(BVGetOrthogonalization(pep->V,&otype,&oref,&oeta,&oblock));
118:     PetscCall(BVSetOrthogonalization(*basis,otype,oref,oeta,oblock));
119:     PetscCall(PetscObjectStateIncrease((PetscObject)*basis));
120:   } else PetscCall(BVDuplicate(pep->V,basis));
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

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

131:   PetscFunctionBegin;
132:   PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
133:   if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
134:   if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;
135:   PetscCheck(pep->which==PEP_TARGET_MAGNITUDE || pep->which==PEP_TARGET_REAL || pep->which==PEP_TARGET_IMAGINARY,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver supports only target which, see PEPSetWhichEigenpairs()");

137:   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STPRECOND,&isprecond));
138:   PetscCheck(isprecond,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver only works with PRECOND spectral transformation");

140:   PetscCall(STGetTransform(pep->st,&flg));
141:   PetscCheck(!flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver requires the ST transform flag unset, see STSetTransform()");
142:   PEPCheckIgnored(pep,PEP_FEATURE_EXTRACT);

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

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

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

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

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

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

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

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

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

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

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

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

407:   PetscFunctionBegin;
408:   PetscCall(PCShellGetContext(pc,&ctx));
409:   PetscCall(VecCompGetSubVecs(x,&sz,&xs));
410:   PetscCall(VecCompGetSubVecs(y,NULL,&ys));
411:   /* y = B\x apply extended PC */
412:   PetscCall(PEPJDExtendedPCApply(pc,xs[0],ys[0]));
413: #if !defined(PETSC_USE_COMPLEX)
414:   if (sz==2) PetscCall(PEPJDExtendedPCApply(pc,xs[1],ys[1]));
415: #endif

417:   /* Compute eta = u'*y / u'*Bp */
418:   PetscCall(VecDot(ys[0],ctx->u[0],&rr));
419:   eta  = -rr*ctx->gamma[0];
420: #if !defined(PETSC_USE_COMPLEX)
421:   if (sz==2) {
422:     PetscCall(VecDot(ys[0],ctx->u[1],&xr));
423:     PetscCall(VecDot(ys[1],ctx->u[0],&rx));
424:     PetscCall(VecDot(ys[1],ctx->u[1],&xx));
425:     eta += -ctx->gamma[0]*xx-ctx->gamma[1]*(-xr+rx);
426:   }
427: #endif
428:   eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];

430:   /* y = y - eta*Bp */
431:   PetscCall(VecAXPY(ys[0],eta,ctx->Bp[0]));
432: #if !defined(PETSC_USE_COMPLEX)
433:   if (sz==2) {
434:     PetscCall(VecAXPY(ys[1],eta,ctx->Bp[1]));
435:     eta = -ctx->gamma[1]*(rr+xx)+ctx->gamma[0]*(-xr+rx);
436:     eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
437:     PetscCall(VecAXPY(ys[0],eta,ctx->Bp[1]));
438:     PetscCall(VecAXPY(ys[1],-eta,ctx->Bp[0]));
439:   }
440: #endif
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

444: static PetscErrorCode PEPJDCopyToExtendedVec(PEP pep,Vec v,PetscScalar *a,PetscInt na,PetscInt off,Vec vex,PetscBool back)
445: {
446:   PetscMPIInt    np,rk,count;
447:   PetscScalar    *array1,*array2;
448:   PetscInt       nloc;

450:   PetscFunctionBegin;
451:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
452:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
453:   PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
454:   if (v) {
455:     PetscCall(VecGetArray(v,&array1));
456:     PetscCall(VecGetArray(vex,&array2));
457:     if (back) PetscCall(PetscArraycpy(array1,array2,nloc));
458:     else PetscCall(PetscArraycpy(array2,array1,nloc));
459:     PetscCall(VecRestoreArray(v,&array1));
460:     PetscCall(VecRestoreArray(vex,&array2));
461:   }
462:   if (a) {
463:     PetscCall(VecGetArray(vex,&array2));
464:     if (back) {
465:       PetscCall(PetscArraycpy(a,array2+nloc+off,na));
466:       PetscCall(PetscMPIIntCast(na,&count));
467:       PetscCallMPI(MPI_Bcast(a,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
468:     } else {
469:       PetscCall(PetscArraycpy(array2+nloc+off,a,na));
470:       PetscCall(PetscMPIIntCast(na,&count));
471:       PetscCallMPI(MPI_Bcast(array2+nloc+off,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
472:     }
473:     PetscCall(VecRestoreArray(vex,&array2));
474:   }
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

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

488:   PetscFunctionBegin;
489:   a = pep->pbc; b=a+pep->nmat; g=b+pep->nmat;
490:   PetscCall(PetscBLASIntCast(n,&n_));
491:   PetscCall(PetscBLASIntCast(ldh,&ldh_));
492:   if (idx<1) PetscCall(PetscArrayzero(q,t?n:n*n));
493:   else if (idx==1) {
494:     if (t) {for (j=0;j<n;j++) q[j] = t[j]*beval[idx-1]/a[0];}
495:     else {
496:       PetscCall(PetscArrayzero(q,n*n));
497:       for (j=0;j<n;j++) q[(j+1)*n] = beval[idx-1]/a[0];
498:     }
499:   } else {
500:     if (t) {
501:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,H,&ldh_,qp,&one,&zero,q,&one));
502:       for (j=0;j<n;j++) {
503:         q[j] += beval[idx-1]*t[j]-b[idx-1]*qp[j]-g[idx-1]*qpp[j];
504:         q[j] /= a[idx-1];
505:       }
506:     } else {
507:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,H,&ldh_,qp,&n_,&zero,q,&n_));
508:       for (j=0;j<n;j++) {
509:         q[j+n*j] += beval[idx-1];
510:         for (i=0;i<n;i++) {
511:           q[i+n*j] += -b[idx-1]*qp[j*n+i]-g[idx-1]*qpp[j*n+i];
512:           q[i+n*j] /= a[idx-1];
513:         }
514:       }
515:     }
516:   }
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

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

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

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

660: static PetscErrorCode PEPJDProcessInitialSpace(PEP pep,Vec *w)
661: {
662:   PEP_JD         *pjd = (PEP_JD*)pep->data;
663:   PetscScalar    *tt,target[2];
664:   Vec            vg,wg;
665:   PetscInt       i;
666:   PetscReal      norm;

668:   PetscFunctionBegin;
669:   PetscCall(PetscMalloc1(pjd->ld-1,&tt));
670:   PetscCheck(pep->nini==0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for initial vectors not implemented yet");
671:   PetscCall(BVSetRandomColumn(pjd->V,0));
672:   for (i=0;i<pjd->ld-1;i++) tt[i] = 0.0;
673:   PetscCall(BVGetColumn(pjd->V,0,&vg));
674:   PetscCall(PEPJDCopyToExtendedVec(pep,NULL,tt,pjd->ld-1,0,vg,PETSC_FALSE));
675:   PetscCall(BVRestoreColumn(pjd->V,0,&vg));
676:   PetscCall(BVNormColumn(pjd->V,0,NORM_2,&norm));
677:   PetscCall(BVScaleColumn(pjd->V,0,1.0/norm));
678:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
679:     PetscCall(BVGetColumn(pjd->V,0,&vg));
680:     PetscCall(BVGetColumn(pjd->W,0,&wg));
681:     PetscCall(VecSet(wg,0.0));
682:     target[0] = pep->target; target[1] = 0.0;
683:     PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,1,&vg,target,&wg,w));
684:     PetscCall(BVRestoreColumn(pjd->W,0,&wg));
685:     PetscCall(BVRestoreColumn(pjd->V,0,&vg));
686:     PetscCall(BVNormColumn(pjd->W,0,NORM_2,&norm));
687:     PetscCall(BVScaleColumn(pjd->W,0,1.0/norm));
688:   }
689:   PetscCall(PetscFree(tt));
690:   PetscFunctionReturn(PETSC_SUCCESS);
691: }

693: static PetscErrorCode MatMult_PEPJD(Mat P,Vec x,Vec y)
694: {
695:   PEP_JD_MATSHELL *matctx;
696:   PEP_JD          *pjd;
697:   PetscInt        i,j,nconv,nloc,nmat,ldt,ncv,sz;
698:   Vec             tx,ty;
699:   const Vec       *xs,*ys;
700:   PetscScalar     *array1,*array2,*x2=NULL,*y2,*tt=NULL,*xx=NULL,*xxi,theta[2],sone=1.0,*qj,*val,*vali=NULL;
701:   PetscBLASInt    n,ld,one=1;
702:   PetscMPIInt     np;
703: #if !defined(PETSC_USE_COMPLEX)
704:   Vec             txi=NULL,tyi=NULL;
705:   PetscScalar     *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
706: #endif

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

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

826: static PetscErrorCode MatCreateVecs_PEPJD(Mat A,Vec *right,Vec *left)
827: {
828:   PEP_JD_MATSHELL *matctx;
829:   PEP_JD          *pjd;
830:   PetscInt        kspsf=1,i;
831:   Vec             v[2];

833:   PetscFunctionBegin;
834:   PetscCall(MatShellGetContext(A,&matctx));
835:   pjd   = (PEP_JD*)(matctx->pep->data);
836: #if !defined (PETSC_USE_COMPLEX)
837:   kspsf = 2;
838: #endif
839:   for (i=0;i<kspsf;i++) PetscCall(BVCreateVec(pjd->V,v+i));
840:   if (right) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,right));
841:   if (left) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,left));
842:   for (i=0;i<kspsf;i++) PetscCall(VecDestroy(&v[i]));
843:   PetscFunctionReturn(PETSC_SUCCESS);
844: }

846: static PetscErrorCode PEPJDUpdateExtendedPC(PEP pep,PetscScalar theta)
847: {
848:   PEP_JD         *pjd = (PEP_JD*)pep->data;
849:   PEP_JD_PCSHELL *pcctx;
850:   PetscInt       i,j,k,n=pjd->nlock,ld=pjd->ld,deg=pep->nmat-1;
851:   PetscScalar    *M,*ps,*work,*U,*V,*S,*Sp,*Spp,snone=-1.0,sone=1.0,zero=0.0,*val;
852:   PetscReal      tol,maxeig=0.0,*sg,*rwork;
853:   PetscBLASInt   n_,info,ld_,*p,lw_,rk=0;

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

893:       /* compute M */
894:       PetscCall(PEPEvaluateBasis(pep,theta,0.0,val,NULL));
895:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&snone,pjd->XpX,&ld_,ps,&ld_,&zero,M,&ld_));
896:       PetscCall(PetscArrayzero(S,2*n*n));
897:       Sp = S+n*n;
898:       for (j=0;j<n;j++) S[j*(n+1)] = 1.0;
899:       for (k=1;k<pjd->midx;k++) {
900:         for (j=0;j<n;j++) for (i=0;i<n;i++) V[j*n+i] = S[j*n+i] - ps[j*ld+i]*val[k];
901:         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pjd->XpX,&ld_,V,&n_,&zero,U,&n_));
902:         PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+k*ld*ld,&ld_,U,&n_,&sone,M,&ld_));
903:         Spp = Sp; Sp = S;
904:         PetscCall(PEPJDEvaluateHatBasis(pep,n,pjd->T,pep->ncv,val,NULL,k+1,Spp,Sp,S));
905:       }
906:     }
907:     /* inverse */
908:     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,M,&ld_,p,&info));
909:     SlepcCheckLapackInfo("getrf",info);
910:     PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&n_,M,&ld_,p,work,&n_,&info));
911:     SlepcCheckLapackInfo("getri",info);
912:     PetscCall(PetscFPTrapPop());
913:     if (pjd->midx==1) PetscCall(PetscFree2(work,p));
914:     else PetscCall(PetscFree7(U,S,sg,work,rwork,p,val));
915:   }
916:   PetscFunctionReturn(PETSC_SUCCESS);
917: }

919: static PetscErrorCode PEPJDMatSetUp(PEP pep,PetscInt sz,PetscScalar *theta)
920: {
921:   PEP_JD          *pjd = (PEP_JD*)pep->data;
922:   PEP_JD_MATSHELL *matctx;
923:   PEP_JD_PCSHELL  *pcctx;
924:   MatStructure    str;
925:   PetscScalar     *vals,*valsi;
926:   PetscBool       skipmat=PETSC_FALSE;
927:   PetscInt        i;
928:   Mat             Pr=NULL;

930:   PetscFunctionBegin;
931:   if (sz==2 && theta[1]==0.0) sz = 1;
932:   PetscCall(MatShellGetContext(pjd->Pshell,&matctx));
933:   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
934:   if (matctx->Pr && matctx->theta[0]==theta[0] && ((!matctx->Pi && sz==1) || (sz==2 && matctx->theta[1]==theta[1]))) {
935:     if (pcctx->n == pjd->nlock) PetscFunctionReturn(PETSC_SUCCESS);
936:     skipmat = PETSC_TRUE;
937:   }
938:   if (!skipmat) {
939:     PetscCall(PetscMalloc2(pep->nmat,&vals,pep->nmat,&valsi));
940:     PetscCall(STGetMatStructure(pep->st,&str));
941:     PetscCall(PEPEvaluateBasis(pep,theta[0],theta[1],vals,valsi));
942:     if (!matctx->Pr) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pr));
943:     else PetscCall(MatCopy(pep->A[0],matctx->Pr,str));
944:     for (i=1;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pr,vals[i],pep->A[i],str));
945:     if (!pjd->reusepc) {
946:       if (pcctx->PPr && sz==2) {
947:         PetscCall(MatCopy(matctx->Pr,pcctx->PPr,str));
948:         Pr = pcctx->PPr;
949:       } else Pr = matctx->Pr;
950:     }
951:     matctx->theta[0] = theta[0];
952: #if !defined(PETSC_USE_COMPLEX)
953:     if (sz==2) {
954:       if (!matctx->Pi) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pi));
955:       else PetscCall(MatCopy(pep->A[1],matctx->Pi,str));
956:       PetscCall(MatScale(matctx->Pi,valsi[1]));
957:       for (i=2;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pi,valsi[i],pep->A[i],str));
958:       matctx->theta[1] = theta[1];
959:     }
960: #endif
961:     PetscCall(PetscFree2(vals,valsi));
962:   }
963:   if (!pjd->reusepc) {
964:     if (!skipmat) {
965:       PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
966:       PetscCall(PCSetUp(pcctx->pc));
967:     }
968:     PetscCall(PEPJDUpdateExtendedPC(pep,theta[0]));
969:   }
970:   PetscFunctionReturn(PETSC_SUCCESS);
971: }

973: static PetscErrorCode PEPJDCreateShellPC(PEP pep,Vec *ww)
974: {
975:   PEP_JD          *pjd = (PEP_JD*)pep->data;
976:   PEP_JD_PCSHELL  *pcctx;
977:   PEP_JD_MATSHELL *matctx;
978:   KSP             ksp;
979:   PetscInt        nloc,mloc,kspsf=1;
980:   Vec             v[2];
981:   PetscScalar     target[2];
982:   Mat             Pr;

984:   PetscFunctionBegin;
985:   /* Create the reference vector */
986:   PetscCall(BVGetColumn(pjd->V,0,&v[0]));
987:   v[1] = v[0];
988: #if !defined (PETSC_USE_COMPLEX)
989:   kspsf = 2;
990: #endif
991:   PetscCall(VecCreateCompWithVecs(v,kspsf,NULL,&pjd->vtempl));
992:   PetscCall(BVRestoreColumn(pjd->V,0,&v[0]));

994:   /* Replace preconditioner with one containing projectors */
995:   PetscCall(PCCreate(PetscObjectComm((PetscObject)pep),&pjd->pcshell));
996:   PetscCall(PCSetType(pjd->pcshell,PCSHELL));
997:   PetscCall(PCShellSetName(pjd->pcshell,"PCPEPJD"));
998:   PetscCall(PCShellSetApply(pjd->pcshell,PCShellApply_PEPJD));
999:   PetscCall(PetscNew(&pcctx));
1000:   PetscCall(PCShellSetContext(pjd->pcshell,pcctx));
1001:   PetscCall(STGetKSP(pep->st,&ksp));
1002:   PetscCall(BVCreateVec(pjd->V,&pcctx->Bp[0]));
1003:   PetscCall(VecDuplicate(pcctx->Bp[0],&pcctx->Bp[1]));
1004:   PetscCall(KSPGetPC(ksp,&pcctx->pc));
1005:   PetscCall(PetscObjectReference((PetscObject)pcctx->pc));
1006:   PetscCall(MatGetLocalSize(pep->A[0],&mloc,&nloc));
1007:   if (pjd->ld>1) {
1008:     nloc += pjd->ld-1; mloc += pjd->ld-1;
1009:   }
1010:   PetscCall(PetscNew(&matctx));
1011:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),kspsf*nloc,kspsf*mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&pjd->Pshell));
1012:   PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_MULT,(void(*)(void))MatMult_PEPJD));
1013:   PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_PEPJD));
1014:   matctx->pep = pep;
1015:   target[0] = pep->target; target[1] = 0.0;
1016:   PetscCall(PEPJDMatSetUp(pep,1,target));
1017:   Pr = matctx->Pr;
1018:   pcctx->PPr = NULL;
1019: #if !defined(PETSC_USE_COMPLEX)
1020:   if (!pjd->reusepc) {
1021:     PetscCall(MatDuplicate(matctx->Pr,MAT_COPY_VALUES,&pcctx->PPr));
1022:     Pr = pcctx->PPr;
1023:   }
1024: #endif
1025:   PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
1026:   PetscCall(PCSetErrorIfFailure(pcctx->pc,PETSC_TRUE));
1027:   PetscCall(KSPSetPC(ksp,pjd->pcshell));
1028:   if (pjd->reusepc) {
1029:     PetscCall(PCSetReusePreconditioner(pcctx->pc,PETSC_TRUE));
1030:     PetscCall(KSPSetReusePreconditioner(ksp,PETSC_TRUE));
1031:   }
1032:   PetscCall(PEP_KSPSetOperators(ksp,pjd->Pshell,pjd->Pshell));
1033:   PetscCall(KSPSetUp(ksp));
1034:   if (pjd->ld>1) {
1035:     PetscCall(PetscMalloc2(pjd->ld*pjd->ld,&pcctx->M,pjd->ld*pjd->ld,&pcctx->ps));
1036:     pcctx->pep = pep;
1037:   }
1038:   matctx->work = ww;
1039:   pcctx->work  = ww;
1040:   PetscFunctionReturn(PETSC_SUCCESS);
1041: }

1043: static PetscErrorCode PEPJDEigenvectors(PEP pep)
1044: {
1045:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1046:   PetscBLASInt   ld,nconv,info,nc;
1047:   PetscScalar    *Z;
1048:   PetscReal      *wr;
1049:   Mat            U;
1050: #if defined(PETSC_USE_COMPLEX)
1051:   PetscScalar    *w;
1052: #endif

1054:   PetscFunctionBegin;
1055:   PetscCall(PetscBLASIntCast(pep->ncv,&ld));
1056:   PetscCall(PetscBLASIntCast(pep->nconv,&nconv));
1057:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1058: #if !defined(PETSC_USE_COMPLEX)
1059:   PetscCall(PetscMalloc2(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr));
1060:   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,wr,&info));
1061: #else
1062:   PetscCall(PetscMalloc3(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr,2*pep->ncv,&w));
1063:   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,w,wr,&info));
1064: #endif
1065:   PetscCall(PetscFPTrapPop());
1066:   SlepcCheckLapackInfo("trevc",info);
1067:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,Z,&U));
1068:   PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv));
1069:   PetscCall(BVMultInPlace(pjd->X,U,0,pep->nconv));
1070:   PetscCall(BVNormalize(pjd->X,pep->eigi));
1071:   PetscCall(MatDestroy(&U));
1072: #if !defined(PETSC_USE_COMPLEX)
1073:   PetscCall(PetscFree2(Z,wr));
1074: #else
1075:   PetscCall(PetscFree3(Z,wr,w));
1076: #endif
1077:   PetscFunctionReturn(PETSC_SUCCESS);
1078: }

1080: static PetscErrorCode PEPJDLockConverged(PEP pep,PetscInt *nv,PetscInt sz)
1081: {
1082:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1083:   PetscInt       j,i,*P,ldds,rk=0,nvv=*nv;
1084:   Vec            v,x,w;
1085:   PetscScalar    *R,*r,*pX,target[2];
1086:   Mat            X;
1087:   PetscBLASInt   sz_,rk_,nv_,info;
1088:   PetscMPIInt    np;

1090:   PetscFunctionBegin;
1091:   /* update AX and XpX */
1092:   for (i=sz;i>0;i--) {
1093:     PetscCall(BVGetColumn(pjd->X,pjd->nlock-i,&x));
1094:     for (j=0;j<pep->nmat;j++) {
1095:       PetscCall(BVGetColumn(pjd->AX[j],pjd->nlock-i,&v));
1096:       PetscCall(MatMult(pep->A[j],x,v));
1097:       PetscCall(BVRestoreColumn(pjd->AX[j],pjd->nlock-i,&v));
1098:       PetscCall(BVSetActiveColumns(pjd->AX[j],0,pjd->nlock-i+1));
1099:     }
1100:     PetscCall(BVRestoreColumn(pjd->X,pjd->nlock-i,&x));
1101:     PetscCall(BVDotColumn(pjd->X,(pjd->nlock-i),pjd->XpX+(pjd->nlock-i)*(pjd->ld)));
1102:     pjd->XpX[(pjd->nlock-i)*(1+pjd->ld)] = 1.0;
1103:     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]);
1104:   }

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

1109:   /* evaluate the polynomial basis in T */
1110:   PetscCall(PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat));
1111:   for (j=0;j<pep->nmat;j++) PetscCall(PEPEvaluateBasisMat(pep,pjd->nlock,pjd->T,pep->ncv,j,(j>1)?pjd->Tj+(j-2)*pjd->ld*pjd->ld:NULL,pjd->ld,j?pjd->Tj+(j-1)*pjd->ld*pjd->ld:NULL,pjd->ld,pjd->Tj+j*pjd->ld*pjd->ld,pjd->ld));

1113:   /* Extend search space */
1114:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
1115:   PetscCall(PetscCalloc3(nvv,&P,nvv*nvv,&R,nvv*sz,&r));
1116:   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
1117:   PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
1118:   PetscCall(PEPJDOrthogonalize(nvv,nvv,pX,ldds,&rk,P,R,nvv));
1119:   for (j=0;j<sz;j++) {
1120:     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 */
1121:   }
1122:   PetscCall(PetscBLASIntCast(rk,&rk_));
1123:   PetscCall(PetscBLASIntCast(sz,&sz_));
1124:   PetscCall(PetscBLASIntCast(nvv,&nv_));
1125:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1126:   PetscCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&rk_,R,&nv_,&info));
1127:   PetscCall(PetscFPTrapPop());
1128:   SlepcCheckLapackInfo("trtri",info);
1129:   for (i=0;i<sz;i++) PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&rk_,R,&nv_,r+i,&sz_));
1130:   for (i=0;i<sz*rk;i++) r[i] = PetscConj(r[i])/PetscSqrtReal(np); /* revert */
1131:   PetscCall(BVSetActiveColumns(pjd->V,0,nvv));
1132:   rk -= sz;
1133:   for (j=0;j<rk;j++) PetscCall(PetscArraycpy(R+j*nvv,pX+(j+sz)*ldds,nvv));
1134:   PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
1135:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nvv,rk,R,&X));
1136:   PetscCall(BVMultInPlace(pjd->V,X,0,rk));
1137:   PetscCall(MatDestroy(&X));
1138:   PetscCall(BVSetActiveColumns(pjd->V,0,rk));
1139:   for (j=0;j<rk;j++) {
1140:     PetscCall(BVGetColumn(pjd->V,j,&v));
1141:     PetscCall(PEPJDCopyToExtendedVec(pep,NULL,r+sz*(j+sz),sz,pjd->nlock-sz,v,PETSC_FALSE));
1142:     PetscCall(BVRestoreColumn(pjd->V,j,&v));
1143:   }
1144:   PetscCall(BVOrthogonalize(pjd->V,NULL));

1146:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
1147:     for (j=0;j<rk;j++) {
1148:       /* W = P(target)*V */
1149:       PetscCall(BVGetColumn(pjd->W,j,&w));
1150:       PetscCall(BVGetColumn(pjd->V,j,&v));
1151:       target[0] = pep->target; target[1] = 0.0;
1152:       PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,1,&v,target,&w,pep->work));
1153:       PetscCall(BVRestoreColumn(pjd->V,j,&v));
1154:       PetscCall(BVRestoreColumn(pjd->W,j,&w));
1155:     }
1156:     PetscCall(BVSetActiveColumns(pjd->W,0,rk));
1157:     PetscCall(BVOrthogonalize(pjd->W,NULL));
1158:   }
1159:   *nv = rk;
1160:   PetscCall(PetscFree3(P,R,r));
1161:   PetscFunctionReturn(PETSC_SUCCESS);
1162: }

1164: static PetscErrorCode PEPJDSystemSetUp(PEP pep,PetscInt sz,PetscScalar *theta,Vec *u,Vec *p,Vec *ww)
1165: {
1166:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1167:   PEP_JD_PCSHELL *pcctx;
1168: #if !defined(PETSC_USE_COMPLEX)
1169:   PetscScalar    s[2];
1170: #endif

1172:   PetscFunctionBegin;
1173:   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
1174:   PetscCall(PEPJDMatSetUp(pep,sz,theta));
1175:   pcctx->u[0] = u[0]; pcctx->u[1] = u[1];
1176:   /* Compute r'. p is a work space vector */
1177:   PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,sz,u,theta,p,ww));
1178:   PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[0],pcctx->Bp[0]));
1179:   PetscCall(VecDot(pcctx->Bp[0],u[0],pcctx->gamma));
1180: #if !defined(PETSC_USE_COMPLEX)
1181:   if (sz==2) {
1182:     PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[1],pcctx->Bp[1]));
1183:     PetscCall(VecDot(pcctx->Bp[0],u[1],pcctx->gamma+1));
1184:     PetscCall(VecMDot(pcctx->Bp[1],2,u,s));
1185:     pcctx->gamma[0] += s[1];
1186:     pcctx->gamma[1] = -pcctx->gamma[1]+s[0];
1187:   }
1188: #endif
1189:   if (sz==1) {
1190:     PetscCall(VecZeroEntries(pcctx->Bp[1]));
1191:     pcctx->gamma[1] = 0.0;
1192:   }
1193:   PetscFunctionReturn(PETSC_SUCCESS);
1194: }

1196: static PetscErrorCode PEPSolve_JD(PEP pep)
1197: {
1198:   PEP_JD          *pjd = (PEP_JD*)pep->data;
1199:   PetscInt        k,nv,nvc,ld,minv,dim,bupdated=0,sz=1,kspsf=1,idx,off,maxits,nloc;
1200:   PetscMPIInt     np,count;
1201:   PetscScalar     theta[2]={0.0,0.0},ritz[2]={0.0,0.0},*pX,*eig,*eigi,*array;
1202:   PetscReal       norm,*res,tol=0.0,rtol,abstol, dtol;
1203:   PetscBool       lindep,ini=PETSC_TRUE;
1204:   Vec             tc,t[2]={NULL,NULL},u[2]={NULL,NULL},p[2]={NULL,NULL};
1205:   Vec             rc,rr[2],r[2]={NULL,NULL},*ww=pep->work,v[2];
1206:   Mat             G,X,Y;
1207:   KSP             ksp;
1208:   PEP_JD_PCSHELL  *pcctx;
1209:   PEP_JD_MATSHELL *matctx;
1210: #if !defined(PETSC_USE_COMPLEX)
1211:   PetscReal       norm1;
1212: #endif

1214:   PetscFunctionBegin;
1215:   PetscCall(PetscCitationsRegister(citation,&cited));
1216:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
1217:   PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
1218:   PetscCall(DSGetLeadingDimension(pep->ds,&ld));
1219:   PetscCall(PetscCalloc3(pep->ncv+pep->nev,&eig,pep->ncv+pep->nev,&eigi,pep->ncv+pep->nev,&res));
1220:   pjd->nlock = 0;
1221:   PetscCall(STGetKSP(pep->st,&ksp));
1222:   PetscCall(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
1223: #if !defined (PETSC_USE_COMPLEX)
1224:   kspsf = 2;
1225: #endif
1226:   PetscCall(PEPJDProcessInitialSpace(pep,ww));
1227:   nv = (pep->nini)?pep->nini:1;

1229:   /* Replace preconditioner with one containing projectors */
1230:   PetscCall(PEPJDCreateShellPC(pep,ww));
1231:   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));

1233:   /* Create auxiliary vectors */
1234:   PetscCall(BVCreateVec(pjd->V,&u[0]));
1235:   PetscCall(VecDuplicate(u[0],&p[0]));
1236:   PetscCall(VecDuplicate(u[0],&r[0]));
1237: #if !defined (PETSC_USE_COMPLEX)
1238:   PetscCall(VecDuplicate(u[0],&u[1]));
1239:   PetscCall(VecDuplicate(u[0],&p[1]));
1240:   PetscCall(VecDuplicate(u[0],&r[1]));
1241: #endif

1243:   /* Restart loop */
1244:   while (pep->reason == PEP_CONVERGED_ITERATING) {
1245:     pep->its++;
1246:     PetscCall(DSSetDimensions(pep->ds,nv,0,0));
1247:     PetscCall(BVSetActiveColumns(pjd->V,bupdated,nv));
1248:     PetscCall(PEPJDUpdateTV(pep,bupdated,nv,ww));
1249:     if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVSetActiveColumns(pjd->W,bupdated,nv));
1250:     for (k=0;k<pep->nmat;k++) {
1251:       PetscCall(BVSetActiveColumns(pjd->TV[k],bupdated,nv));
1252:       PetscCall(DSGetMat(pep->ds,DSMatExtra[k],&G));
1253:       PetscCall(BVMatProject(pjd->TV[k],NULL,pjd->W,G));
1254:       PetscCall(DSRestoreMat(pep->ds,DSMatExtra[k],&G));
1255:     }
1256:     PetscCall(BVSetActiveColumns(pjd->V,0,nv));
1257:     PetscCall(BVSetActiveColumns(pjd->W,0,nv));

1259:     /* Solve projected problem */
1260:     PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
1261:     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
1262:     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
1263:     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
1264:     idx = 0;
1265:     do {
1266:       ritz[0] = pep->eigr[idx];
1267: #if !defined(PETSC_USE_COMPLEX)
1268:       ritz[1] = pep->eigi[idx];
1269:       sz = (ritz[1]==0.0)?1:2;
1270: #endif
1271:       /* Compute Ritz vector u=V*X(:,1) */
1272:       PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
1273:       PetscCall(BVSetActiveColumns(pjd->V,0,nv));
1274:       PetscCall(BVMultVec(pjd->V,1.0,0.0,u[0],pX+idx*ld));
1275: #if !defined(PETSC_USE_COMPLEX)
1276:       if (sz==2) PetscCall(BVMultVec(pjd->V,1.0,0.0,u[1],pX+(idx+1)*ld));
1277: #endif
1278:       PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
1279:       PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,sz,u,ritz,r,ww));
1280:       /* Check convergence */
1281:       PetscCall(VecNorm(r[0],NORM_2,&norm));
1282: #if !defined(PETSC_USE_COMPLEX)
1283:       if (sz==2) {
1284:         PetscCall(VecNorm(r[1],NORM_2,&norm1));
1285:         norm = SlepcAbs(norm,norm1);
1286:       }
1287: #endif
1288:       PetscCall((*pep->converged)(pep,ritz[0],ritz[1],norm,&pep->errest[pep->nconv],pep->convergedctx));
1289:       if (sz==2) pep->errest[pep->nconv+1] = pep->errest[pep->nconv];
1290:       if (ini) {
1291:         tol = PetscMin(.1,pep->errest[pep->nconv]); ini = PETSC_FALSE;
1292:       } else tol = PetscMin(pep->errest[pep->nconv],tol/2);
1293:       PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,(pep->errest[pep->nconv]<pep->tol)?pep->nconv+sz:pep->nconv,pep->nev,&pep->reason,pep->stoppingctx));
1294:       if (pep->errest[pep->nconv]<pep->tol) {
1295:         /* Ritz pair converged */
1296:         ini = PETSC_TRUE;
1297:         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
1298:         if (pjd->ld>1) {
1299:           PetscCall(BVGetColumn(pjd->X,pep->nconv,&v[0]));
1300:           PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*pep->nconv,pjd->ld-1,0,u[0],PETSC_TRUE));
1301:           PetscCall(BVRestoreColumn(pjd->X,pep->nconv,&v[0]));
1302:           PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+1));
1303:           PetscCall(BVNormColumn(pjd->X,pep->nconv,NORM_2,&norm));
1304:           PetscCall(BVScaleColumn(pjd->X,pep->nconv,1.0/norm));
1305:           for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*pep->nconv+k] *= PetscSqrtReal(np)/norm;
1306:           pjd->T[(pep->ncv+1)*pep->nconv] = ritz[0];
1307:           eig[pep->nconv] = ritz[0];
1308:           idx++;
1309: #if !defined(PETSC_USE_COMPLEX)
1310:           if (sz==2) {
1311:             PetscCall(BVGetColumn(pjd->X,pep->nconv+1,&v[0]));
1312:             PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*(pep->nconv+1),pjd->ld-1,0,u[1],PETSC_TRUE));
1313:             PetscCall(BVRestoreColumn(pjd->X,pep->nconv+1,&v[0]));
1314:             PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+2));
1315:             PetscCall(BVNormColumn(pjd->X,pep->nconv+1,NORM_2,&norm1));
1316:             PetscCall(BVScaleColumn(pjd->X,pep->nconv+1,1.0/norm1));
1317:             for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*(pep->nconv+1)+k] *= PetscSqrtReal(np)/norm1;
1318:             pjd->T[(pep->ncv+1)*(pep->nconv+1)] = ritz[0];
1319:             pjd->T[(pep->ncv+1)*pep->nconv+1] = -ritz[1]*norm1/norm;
1320:             pjd->T[(pep->ncv+1)*(pep->nconv+1)-1] = ritz[1]*norm/norm1;
1321:             eig[pep->nconv+1] = ritz[0];
1322:             eigi[pep->nconv] = ritz[1]; eigi[pep->nconv+1] = -ritz[1];
1323:             idx++;
1324:           }
1325: #endif
1326:         } else PetscCall(BVInsertVec(pep->V,pep->nconv,u[0]));
1327:         pep->nconv += sz;
1328:       }
1329:     } while (pep->errest[pep->nconv]<pep->tol && pep->nconv<nv);

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

1465: static PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
1466: {
1467:   PEP_JD *pjd = (PEP_JD*)pep->data;

1469:   PetscFunctionBegin;
1470:   if (keep==(PetscReal)PETSC_DEFAULT) pjd->keep = 0.5;
1471:   else {
1472:     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
1473:     pjd->keep = keep;
1474:   }
1475:   PetscFunctionReturn(PETSC_SUCCESS);
1476: }

1478: /*@
1479:    PEPJDSetRestart - Sets the restart parameter for the Jacobi-Davidson
1480:    method, in particular the proportion of basis vectors that must be kept
1481:    after restart.

1483:    Logically Collective

1485:    Input Parameters:
1486: +  pep  - the eigenproblem solver context
1487: -  keep - the number of vectors to be kept at restart

1489:    Options Database Key:
1490: .  -pep_jd_restart - Sets the restart parameter

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

1495:    Level: advanced

1497: .seealso: PEPJDGetRestart()
1498: @*/
1499: PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
1500: {
1501:   PetscFunctionBegin;
1504:   PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
1505:   PetscFunctionReturn(PETSC_SUCCESS);
1506: }

1508: static PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
1509: {
1510:   PEP_JD *pjd = (PEP_JD*)pep->data;

1512:   PetscFunctionBegin;
1513:   *keep = pjd->keep;
1514:   PetscFunctionReturn(PETSC_SUCCESS);
1515: }

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

1520:    Not Collective

1522:    Input Parameter:
1523: .  pep - the eigenproblem solver context

1525:    Output Parameter:
1526: .  keep - the restart parameter

1528:    Level: advanced

1530: .seealso: PEPJDSetRestart()
1531: @*/
1532: PetscErrorCode PEPJDGetRestart(PEP pep,PetscReal *keep)
1533: {
1534:   PetscFunctionBegin;
1536:   PetscAssertPointer(keep,2);
1537:   PetscUseMethod(pep,"PEPJDGetRestart_C",(PEP,PetscReal*),(pep,keep));
1538:   PetscFunctionReturn(PETSC_SUCCESS);
1539: }

1541: static PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
1542: {
1543:   PEP_JD *pjd = (PEP_JD*)pep->data;

1545:   PetscFunctionBegin;
1546:   if (fix == (PetscReal)PETSC_DEFAULT || fix == (PetscReal)PETSC_DECIDE) pjd->fix = 0.01;
1547:   else {
1548:     PetscCheck(fix>=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid fix value, must be >0");
1549:     pjd->fix = fix;
1550:   }
1551:   PetscFunctionReturn(PETSC_SUCCESS);
1552: }

1554: /*@
1555:    PEPJDSetFix - Sets the threshold for changing the target in the correction
1556:    equation.

1558:    Logically Collective

1560:    Input Parameters:
1561: +  pep - the eigenproblem solver context
1562: -  fix - threshold for changing the target

1564:    Options Database Key:
1565: .  -pep_jd_fix - the fix value

1567:    Note:
1568:    The target in the correction equation is fixed at the first iterations.
1569:    When the norm of the residual vector is lower than the fix value,
1570:    the target is set to the corresponding eigenvalue.

1572:    Level: advanced

1574: .seealso: PEPJDGetFix()
1575: @*/
1576: PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
1577: {
1578:   PetscFunctionBegin;
1581:   PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
1582:   PetscFunctionReturn(PETSC_SUCCESS);
1583: }

1585: static PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
1586: {
1587:   PEP_JD *pjd = (PEP_JD*)pep->data;

1589:   PetscFunctionBegin;
1590:   *fix = pjd->fix;
1591:   PetscFunctionReturn(PETSC_SUCCESS);
1592: }

1594: /*@
1595:    PEPJDGetFix - Returns the threshold for changing the target in the correction
1596:    equation.

1598:    Not Collective

1600:    Input Parameter:
1601: .  pep - the eigenproblem solver context

1603:    Output Parameter:
1604: .  fix - threshold for changing the target

1606:    Note:
1607:    The target in the correction equation is fixed at the first iterations.
1608:    When the norm of the residual vector is lower than the fix value,
1609:    the target is set to the corresponding eigenvalue.

1611:    Level: advanced

1613: .seealso: PEPJDSetFix()
1614: @*/
1615: PetscErrorCode PEPJDGetFix(PEP pep,PetscReal *fix)
1616: {
1617:   PetscFunctionBegin;
1619:   PetscAssertPointer(fix,2);
1620:   PetscUseMethod(pep,"PEPJDGetFix_C",(PEP,PetscReal*),(pep,fix));
1621:   PetscFunctionReturn(PETSC_SUCCESS);
1622: }

1624: static PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
1625: {
1626:   PEP_JD *pjd = (PEP_JD*)pep->data;

1628:   PetscFunctionBegin;
1629:   pjd->reusepc = reusepc;
1630:   PetscFunctionReturn(PETSC_SUCCESS);
1631: }

1633: /*@
1634:    PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
1635:    must be reused or not.

1637:    Logically Collective

1639:    Input Parameters:
1640: +  pep     - the eigenproblem solver context
1641: -  reusepc - the reuse flag

1643:    Options Database Key:
1644: .  -pep_jd_reuse_preconditioner - the reuse flag

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

1651:    Level: advanced

1653: .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
1654: @*/
1655: PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
1656: {
1657:   PetscFunctionBegin;
1660:   PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
1661:   PetscFunctionReturn(PETSC_SUCCESS);
1662: }

1664: static PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
1665: {
1666:   PEP_JD *pjd = (PEP_JD*)pep->data;

1668:   PetscFunctionBegin;
1669:   *reusepc = pjd->reusepc;
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

1673: /*@
1674:    PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.

1676:    Not Collective

1678:    Input Parameter:
1679: .  pep - the eigenproblem solver context

1681:    Output Parameter:
1682: .  reusepc - the reuse flag

1684:    Level: advanced

1686: .seealso: PEPJDSetReusePreconditioner()
1687: @*/
1688: PetscErrorCode PEPJDGetReusePreconditioner(PEP pep,PetscBool *reusepc)
1689: {
1690:   PetscFunctionBegin;
1692:   PetscAssertPointer(reusepc,2);
1693:   PetscUseMethod(pep,"PEPJDGetReusePreconditioner_C",(PEP,PetscBool*),(pep,reusepc));
1694:   PetscFunctionReturn(PETSC_SUCCESS);
1695: }

1697: static PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
1698: {
1699:   PEP_JD *pjd = (PEP_JD*)pep->data;

1701:   PetscFunctionBegin;
1702:   if (mmidx == PETSC_DEFAULT || mmidx == PETSC_DECIDE) {
1703:     if (pjd->mmidx != 1) pep->state = PEP_STATE_INITIAL;
1704:     pjd->mmidx = 1;
1705:   } else {
1706:     PetscCheck(mmidx>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mmidx value, should be >0");
1707:     if (pjd->mmidx != mmidx) pep->state = PEP_STATE_INITIAL;
1708:     pjd->mmidx = mmidx;
1709:   }
1710:   PetscFunctionReturn(PETSC_SUCCESS);
1711: }

1713: /*@
1714:    PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.

1716:    Logically Collective

1718:    Input Parameters:
1719: +  pep   - the eigenproblem solver context
1720: -  mmidx - maximum minimality index

1722:    Options Database Key:
1723: .  -pep_jd_minimality_index - the minimality index value

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

1729:    Level: advanced

1731: .seealso: PEPJDGetMinimalityIndex()
1732: @*/
1733: PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
1734: {
1735:   PetscFunctionBegin;
1738:   PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
1739:   PetscFunctionReturn(PETSC_SUCCESS);
1740: }

1742: static PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
1743: {
1744:   PEP_JD *pjd = (PEP_JD*)pep->data;

1746:   PetscFunctionBegin;
1747:   *mmidx = pjd->mmidx;
1748:   PetscFunctionReturn(PETSC_SUCCESS);
1749: }

1751: /*@
1752:    PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
1753:    index.

1755:    Not Collective

1757:    Input Parameter:
1758: .  pep - the eigenproblem solver context

1760:    Output Parameter:
1761: .  mmidx - minimality index

1763:    Level: advanced

1765: .seealso: PEPJDSetMinimalityIndex()
1766: @*/
1767: PetscErrorCode PEPJDGetMinimalityIndex(PEP pep,PetscInt *mmidx)
1768: {
1769:   PetscFunctionBegin;
1771:   PetscAssertPointer(mmidx,2);
1772:   PetscUseMethod(pep,"PEPJDGetMinimalityIndex_C",(PEP,PetscInt*),(pep,mmidx));
1773:   PetscFunctionReturn(PETSC_SUCCESS);
1774: }

1776: static PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
1777: {
1778:   PEP_JD *pjd = (PEP_JD*)pep->data;

1780:   PetscFunctionBegin;
1781:   switch (proj) {
1782:     case PEP_JD_PROJECTION_HARMONIC:
1783:     case PEP_JD_PROJECTION_ORTHOGONAL:
1784:       if (pjd->proj != proj) {
1785:         pep->state = PEP_STATE_INITIAL;
1786:         pjd->proj = proj;
1787:       }
1788:       break;
1789:     default:
1790:       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'proj' value");
1791:   }
1792:   PetscFunctionReturn(PETSC_SUCCESS);
1793: }

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

1798:    Logically Collective

1800:    Input Parameters:
1801: +  pep  - the eigenproblem solver context
1802: -  proj - the type of projection

1804:    Options Database Key:
1805: .  -pep_jd_projection - the projection type, either orthogonal or harmonic

1807:    Level: advanced

1809: .seealso: PEPJDGetProjection()
1810: @*/
1811: PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
1812: {
1813:   PetscFunctionBegin;
1816:   PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
1817:   PetscFunctionReturn(PETSC_SUCCESS);
1818: }

1820: static PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
1821: {
1822:   PEP_JD *pjd = (PEP_JD*)pep->data;

1824:   PetscFunctionBegin;
1825:   *proj = pjd->proj;
1826:   PetscFunctionReturn(PETSC_SUCCESS);
1827: }

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

1832:    Not Collective

1834:    Input Parameter:
1835: .  pep - the eigenproblem solver context

1837:    Output Parameter:
1838: .  proj - the type of projection

1840:    Level: advanced

1842: .seealso: PEPJDSetProjection()
1843: @*/
1844: PetscErrorCode PEPJDGetProjection(PEP pep,PEPJDProjection *proj)
1845: {
1846:   PetscFunctionBegin;
1848:   PetscAssertPointer(proj,2);
1849:   PetscUseMethod(pep,"PEPJDGetProjection_C",(PEP,PEPJDProjection*),(pep,proj));
1850:   PetscFunctionReturn(PETSC_SUCCESS);
1851: }

1853: static PetscErrorCode PEPSetFromOptions_JD(PEP pep,PetscOptionItems *PetscOptionsObject)
1854: {
1855:   PetscBool       flg,b1;
1856:   PetscReal       r1;
1857:   PetscInt        i1;
1858:   PEPJDProjection proj;

1860:   PetscFunctionBegin;
1861:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP JD Options");

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

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

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

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

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

1878:   PetscOptionsHeadEnd();
1879:   PetscFunctionReturn(PETSC_SUCCESS);
1880: }

1882: static PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
1883: {
1884:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1885:   PetscBool      isascii;

1887:   PetscFunctionBegin;
1888:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1889:   if (isascii) {
1890:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*pjd->keep)));
1891:     PetscCall(PetscViewerASCIIPrintf(viewer,"  threshold for changing the target in the correction equation (fix): %g\n",(double)pjd->fix));
1892:     PetscCall(PetscViewerASCIIPrintf(viewer,"  projection type: %s\n",PEPJDProjectionTypes[pjd->proj]));
1893:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum allowed minimality index: %" PetscInt_FMT "\n",pjd->mmidx));
1894:     if (pjd->reusepc) PetscCall(PetscViewerASCIIPrintf(viewer,"  reusing the preconditioner\n"));
1895:   }
1896:   PetscFunctionReturn(PETSC_SUCCESS);
1897: }

1899: static PetscErrorCode PEPSetDefaultST_JD(PEP pep)
1900: {
1901:   KSP            ksp;

1903:   PetscFunctionBegin;
1904:   if (!((PetscObject)pep->st)->type_name) {
1905:     PetscCall(STSetType(pep->st,STPRECOND));
1906:     PetscCall(STPrecondSetKSPHasMat(pep->st,PETSC_TRUE));
1907:   }
1908:   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
1909:   PetscCall(STGetKSP(pep->st,&ksp));
1910:   if (!((PetscObject)ksp)->type_name) {
1911:     PetscCall(KSPSetType(ksp,KSPBCGSL));
1912:     PetscCall(KSPSetTolerances(ksp,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,100));
1913:   }
1914:   PetscFunctionReturn(PETSC_SUCCESS);
1915: }

1917: static PetscErrorCode PEPReset_JD(PEP pep)
1918: {
1919:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1920:   PetscInt       i;

1922:   PetscFunctionBegin;
1923:   for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->TV+i));
1924:   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVDestroy(&pjd->W));
1925:   if (pjd->ld>1) {
1926:     PetscCall(BVDestroy(&pjd->V));
1927:     for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->AX+i));
1928:     PetscCall(BVDestroy(&pjd->N[0]));
1929:     PetscCall(BVDestroy(&pjd->N[1]));
1930:     PetscCall(PetscFree3(pjd->XpX,pjd->T,pjd->Tj));
1931:   }
1932:   PetscCall(PetscFree2(pjd->TV,pjd->AX));
1933:   PetscFunctionReturn(PETSC_SUCCESS);
1934: }

1936: static PetscErrorCode PEPDestroy_JD(PEP pep)
1937: {
1938:   PetscFunctionBegin;
1939:   PetscCall(PetscFree(pep->data));
1940:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",NULL));
1941:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",NULL));
1942:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",NULL));
1943:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",NULL));
1944:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",NULL));
1945:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",NULL));
1946:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",NULL));
1947:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",NULL));
1948:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",NULL));
1949:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",NULL));
1950:   PetscFunctionReturn(PETSC_SUCCESS);
1951: }

1953: SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
1954: {
1955:   PEP_JD         *pjd;

1957:   PetscFunctionBegin;
1958:   PetscCall(PetscNew(&pjd));
1959:   pep->data = (void*)pjd;

1961:   pep->lineariz = PETSC_FALSE;
1962:   pjd->fix      = 0.01;
1963:   pjd->mmidx    = 0;

1965:   pep->ops->solve          = PEPSolve_JD;
1966:   pep->ops->setup          = PEPSetUp_JD;
1967:   pep->ops->setfromoptions = PEPSetFromOptions_JD;
1968:   pep->ops->destroy        = PEPDestroy_JD;
1969:   pep->ops->reset          = PEPReset_JD;
1970:   pep->ops->view           = PEPView_JD;
1971:   pep->ops->setdefaultst   = PEPSetDefaultST_JD;

1973:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",PEPJDSetRestart_JD));
1974:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",PEPJDGetRestart_JD));
1975:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",PEPJDSetFix_JD));
1976:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",PEPJDGetFix_JD));
1977:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",PEPJDSetReusePreconditioner_JD));
1978:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",PEPJDGetReusePreconditioner_JD));
1979:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",PEPJDSetMinimalityIndex_JD));
1980:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",PEPJDGetMinimalityIndex_JD));
1981:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",PEPJDSetProjection_JD));
1982:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",PEPJDGetProjection_JD));
1983:   PetscFunctionReturn(PETSC_SUCCESS);
1984: }