Actual source code: pjd.c

slepc-main 2024-11-15
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;
108:   VecType            vtype;

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

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

134:   PetscFunctionBegin;
135:   PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
136:   if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
137:   if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;
138:   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()");

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1112:   /* evaluate the polynomial basis in T */
1113:   PetscCall(PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat));
1114:   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));

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

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

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

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

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

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

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

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

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

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

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

1468: static PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
1469: {
1470:   PEP_JD *pjd = (PEP_JD*)pep->data;

1472:   PetscFunctionBegin;
1473:   if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) pjd->keep = 0.5;
1474:   else {
1475:     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]");
1476:     pjd->keep = keep;
1477:   }
1478:   PetscFunctionReturn(PETSC_SUCCESS);
1479: }

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

1486:    Logically Collective

1488:    Input Parameters:
1489: +  pep  - the eigenproblem solver context
1490: -  keep - the number of vectors to be kept at restart

1492:    Options Database Key:
1493: .  -pep_jd_restart - Sets the restart parameter

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

1498:    Level: advanced

1500: .seealso: PEPJDGetRestart()
1501: @*/
1502: PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
1503: {
1504:   PetscFunctionBegin;
1507:   PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
1508:   PetscFunctionReturn(PETSC_SUCCESS);
1509: }

1511: static PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
1512: {
1513:   PEP_JD *pjd = (PEP_JD*)pep->data;

1515:   PetscFunctionBegin;
1516:   *keep = pjd->keep;
1517:   PetscFunctionReturn(PETSC_SUCCESS);
1518: }

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

1523:    Not Collective

1525:    Input Parameter:
1526: .  pep - the eigenproblem solver context

1528:    Output Parameter:
1529: .  keep - the restart parameter

1531:    Level: advanced

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

1544: static PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
1545: {
1546:   PEP_JD *pjd = (PEP_JD*)pep->data;

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

1557: /*@
1558:    PEPJDSetFix - Sets the threshold for changing the target in the correction
1559:    equation.

1561:    Logically Collective

1563:    Input Parameters:
1564: +  pep - the eigenproblem solver context
1565: -  fix - threshold for changing the target

1567:    Options Database Key:
1568: .  -pep_jd_fix - the fix value

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

1575:    Level: advanced

1577: .seealso: PEPJDGetFix()
1578: @*/
1579: PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
1580: {
1581:   PetscFunctionBegin;
1584:   PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
1585:   PetscFunctionReturn(PETSC_SUCCESS);
1586: }

1588: static PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
1589: {
1590:   PEP_JD *pjd = (PEP_JD*)pep->data;

1592:   PetscFunctionBegin;
1593:   *fix = pjd->fix;
1594:   PetscFunctionReturn(PETSC_SUCCESS);
1595: }

1597: /*@
1598:    PEPJDGetFix - Returns the threshold for changing the target in the correction
1599:    equation.

1601:    Not Collective

1603:    Input Parameter:
1604: .  pep - the eigenproblem solver context

1606:    Output Parameter:
1607: .  fix - threshold for changing the target

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

1614:    Level: advanced

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

1627: static PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
1628: {
1629:   PEP_JD *pjd = (PEP_JD*)pep->data;

1631:   PetscFunctionBegin;
1632:   pjd->reusepc = reusepc;
1633:   PetscFunctionReturn(PETSC_SUCCESS);
1634: }

1636: /*@
1637:    PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
1638:    must be reused or not.

1640:    Logically Collective

1642:    Input Parameters:
1643: +  pep     - the eigenproblem solver context
1644: -  reusepc - the reuse flag

1646:    Options Database Key:
1647: .  -pep_jd_reuse_preconditioner - the reuse flag

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

1654:    Level: advanced

1656: .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
1657: @*/
1658: PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
1659: {
1660:   PetscFunctionBegin;
1663:   PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
1664:   PetscFunctionReturn(PETSC_SUCCESS);
1665: }

1667: static PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
1668: {
1669:   PEP_JD *pjd = (PEP_JD*)pep->data;

1671:   PetscFunctionBegin;
1672:   *reusepc = pjd->reusepc;
1673:   PetscFunctionReturn(PETSC_SUCCESS);
1674: }

1676: /*@
1677:    PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.

1679:    Not Collective

1681:    Input Parameter:
1682: .  pep - the eigenproblem solver context

1684:    Output Parameter:
1685: .  reusepc - the reuse flag

1687:    Level: advanced

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

1700: static PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
1701: {
1702:   PEP_JD *pjd = (PEP_JD*)pep->data;

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

1716: /*@
1717:    PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.

1719:    Logically Collective

1721:    Input Parameters:
1722: +  pep   - the eigenproblem solver context
1723: -  mmidx - maximum minimality index

1725:    Options Database Key:
1726: .  -pep_jd_minimality_index - the minimality index value

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

1732:    Level: advanced

1734: .seealso: PEPJDGetMinimalityIndex()
1735: @*/
1736: PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
1737: {
1738:   PetscFunctionBegin;
1741:   PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
1742:   PetscFunctionReturn(PETSC_SUCCESS);
1743: }

1745: static PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
1746: {
1747:   PEP_JD *pjd = (PEP_JD*)pep->data;

1749:   PetscFunctionBegin;
1750:   *mmidx = pjd->mmidx;
1751:   PetscFunctionReturn(PETSC_SUCCESS);
1752: }

1754: /*@
1755:    PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
1756:    index.

1758:    Not Collective

1760:    Input Parameter:
1761: .  pep - the eigenproblem solver context

1763:    Output Parameter:
1764: .  mmidx - minimality index

1766:    Level: advanced

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

1779: static PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
1780: {
1781:   PEP_JD *pjd = (PEP_JD*)pep->data;

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

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

1801:    Logically Collective

1803:    Input Parameters:
1804: +  pep  - the eigenproblem solver context
1805: -  proj - the type of projection

1807:    Options Database Key:
1808: .  -pep_jd_projection - the projection type, either orthogonal or harmonic

1810:    Level: advanced

1812: .seealso: PEPJDGetProjection()
1813: @*/
1814: PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
1815: {
1816:   PetscFunctionBegin;
1819:   PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
1820:   PetscFunctionReturn(PETSC_SUCCESS);
1821: }

1823: static PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
1824: {
1825:   PEP_JD *pjd = (PEP_JD*)pep->data;

1827:   PetscFunctionBegin;
1828:   *proj = pjd->proj;
1829:   PetscFunctionReturn(PETSC_SUCCESS);
1830: }

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

1835:    Not Collective

1837:    Input Parameter:
1838: .  pep - the eigenproblem solver context

1840:    Output Parameter:
1841: .  proj - the type of projection

1843:    Level: advanced

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

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

1863:   PetscFunctionBegin;
1864:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP JD Options");

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

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

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

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

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

1881:   PetscOptionsHeadEnd();
1882:   PetscFunctionReturn(PETSC_SUCCESS);
1883: }

1885: static PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
1886: {
1887:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1888:   PetscBool      isascii;

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

1902: static PetscErrorCode PEPSetDefaultST_JD(PEP pep)
1903: {
1904:   KSP            ksp;

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

1920: static PetscErrorCode PEPReset_JD(PEP pep)
1921: {
1922:   PEP_JD         *pjd = (PEP_JD*)pep->data;
1923:   PetscInt       i;

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

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

1956: SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
1957: {
1958:   PEP_JD         *pjd;

1960:   PetscFunctionBegin;
1961:   PetscCall(PetscNew(&pjd));
1962:   pep->data = (void*)pjd;

1964:   pep->lineariz = PETSC_FALSE;
1965:   pjd->fix      = 0.01;
1966:   pjd->mmidx    = 0;

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

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