Actual source code: ciss.c

slepc-3.20.2 2024-03-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 eigensolver: "ciss"

 13:    Method: Contour Integral Spectral Slicing

 15:    Algorithm:

 17:        Contour integral based on Sakurai-Sugiura method to construct a
 18:        subspace, with various eigenpair extractions (Rayleigh-Ritz,
 19:        explicit moment).

 21:    Based on code contributed by Y. Maeda, T. Sakurai.

 23:    References:

 25:        [1] T. Sakurai and H. Sugiura, "A projection method for generalized
 26:            eigenvalue problems", J. Comput. Appl. Math. 159:119-128, 2003.

 28:        [2] T. Sakurai and H. Tadano, "CIRR: a Rayleigh-Ritz type method with
 29:            contour integral for generalized eigenvalue problems", Hokkaido
 30:            Math. J. 36:745-757, 2007.
 31: */

 33: #include <slepc/private/epsimpl.h>
 34: #include <slepc/private/slepccontour.h>
 35: #include <slepcblaslapack.h>

 37: typedef struct {
 38:   /* user parameters */
 39:   PetscInt          N;          /* number of integration points (32) */
 40:   PetscInt          L;          /* block size (16) */
 41:   PetscInt          M;          /* moment degree (N/4 = 4) */
 42:   PetscReal         delta;      /* threshold of singular value (1e-12) */
 43:   PetscInt          L_max;      /* maximum number of columns of the source matrix V */
 44:   PetscReal         spurious_threshold; /* discard spurious eigenpairs */
 45:   PetscBool         isreal;     /* A and B are real */
 46:   PetscInt          npart;      /* number of partitions */
 47:   PetscInt          refine_inner;
 48:   PetscInt          refine_blocksize;
 49:   EPSCISSQuadRule   quad;
 50:   EPSCISSExtraction extraction;
 51:   PetscBool         usest;
 52:   /* private data */
 53:   SlepcContourData  contour;
 54:   PetscReal         *sigma;     /* threshold for numerical rank */
 55:   PetscScalar       *weight;
 56:   PetscScalar       *omega;
 57:   PetscScalar       *pp;
 58:   BV                V;
 59:   BV                S;
 60:   BV                pV;
 61:   BV                Y;
 62:   PetscBool         useconj;
 63:   PetscBool         usest_set;  /* whether the user set the usest flag or not */
 64:   PetscObjectId     rgid;
 65:   PetscObjectState  rgstate;
 66: } EPS_CISS;

 68: /*
 69:   Set up KSP solvers for every integration point, only called if !ctx->usest
 70: */
 71: static PetscErrorCode EPSCISSSetUp(EPS eps,Mat A,Mat B,Mat Pa,Mat Pb)
 72: {
 73:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
 74:   SlepcContourData contour;
 75:   PetscInt         i,p_id,nsplit;
 76:   Mat              Amat,Pmat;
 77:   MatStructure     str,strp;

 79:   PetscFunctionBegin;
 80:   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
 81:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
 82:   contour = ctx->contour;
 83:   PetscCall(STGetMatStructure(eps->st,&str));
 84:   PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,&strp));
 85:   for (i=0;i<contour->npoints;i++) {
 86:     p_id = i*contour->subcomm->n + contour->subcomm->color;
 87:     PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Amat));
 88:     if (B) PetscCall(MatAXPY(Amat,-ctx->omega[p_id],B,str));
 89:     else PetscCall(MatShift(Amat,-ctx->omega[p_id]));
 90:     if (nsplit) {
 91:       PetscCall(MatDuplicate(Pa,MAT_COPY_VALUES,&Pmat));
 92:       if (Pb) PetscCall(MatAXPY(Pmat,-ctx->omega[p_id],Pb,strp));
 93:       else PetscCall(MatShift(Pmat,-ctx->omega[p_id]));
 94:     } else Pmat = Amat;
 95:     PetscCall(EPS_KSPSetOperators(contour->ksp[i],Amat,Amat));
 96:     PetscCall(MatDestroy(&Amat));
 97:     if (nsplit) PetscCall(MatDestroy(&Pmat));
 98:   }
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }

102: /*
103:   Y_i = (A-z_i B)^{-1}BV for every integration point, Y=[Y_i] is in the context
104: */
105: static PetscErrorCode EPSCISSSolve(EPS eps,Mat B,BV V,PetscInt L_start,PetscInt L_end)
106: {
107:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
108:   SlepcContourData contour;
109:   PetscInt         i,p_id;
110:   Mat              MV,BMV=NULL,MC;
111:   KSP              ksp;

113:   PetscFunctionBegin;
114:   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
115:   contour = ctx->contour;
116:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
117:   PetscCall(BVSetActiveColumns(V,L_start,L_end));
118:   PetscCall(BVGetMat(V,&MV));
119:   for (i=0;i<contour->npoints;i++) {
120:     p_id = i*contour->subcomm->n + contour->subcomm->color;
121:     if (ctx->usest)  {
122:       PetscCall(STSetShift(eps->st,ctx->omega[p_id]));
123:       PetscCall(STGetKSP(eps->st,&ksp));
124:     } else ksp = contour->ksp[i];
125:     PetscCall(BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end));
126:     PetscCall(BVGetMat(ctx->Y,&MC));
127:     if (B) {
128:       if (!i) {
129:         PetscCall(MatProductCreate(B,MV,NULL,&BMV));
130:         PetscCall(MatProductSetType(BMV,MATPRODUCT_AB));
131:         PetscCall(MatProductSetFromOptions(BMV));
132:         PetscCall(MatProductSymbolic(BMV));
133:       }
134:       PetscCall(MatProductNumeric(BMV));
135:       PetscCall(KSPMatSolve(ksp,BMV,MC));
136:     } else PetscCall(KSPMatSolve(ksp,MV,MC));
137:     PetscCall(BVRestoreMat(ctx->Y,&MC));
138:     if (ctx->usest && i<contour->npoints-1) PetscCall(KSPReset(ksp));
139:   }
140:   PetscCall(MatDestroy(&BMV));
141:   PetscCall(BVRestoreMat(V,&MV));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
146: {
147:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
148:   PetscInt       i;
149:   PetscScalar    center;
150:   PetscReal      radius,a,b,c,d,rgscale;
151: #if defined(PETSC_USE_COMPLEX)
152:   PetscReal      start_ang,end_ang,vscale,theta;
153: #endif
154:   PetscBool      isring,isellipse,isinterval;

156:   PetscFunctionBegin;
157:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
158:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
159:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
160:   PetscCall(RGGetScale(eps->rg,&rgscale));
161:   if (isinterval) {
162:     PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
163:     if (c==d) {
164:       for (i=0;i<nv;i++) {
165: #if defined(PETSC_USE_COMPLEX)
166:         eps->eigr[i] = PetscRealPart(eps->eigr[i]);
167: #else
168:         eps->eigi[i] = 0;
169: #endif
170:       }
171:     }
172:   }
173:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
174:     if (isellipse) {
175:       PetscCall(RGEllipseGetParameters(eps->rg,&center,&radius,NULL));
176:       for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
177:     } else if (isinterval) {
178:       PetscCall(RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d));
179:       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
180:         for (i=0;i<nv;i++) {
181:           if (c==d) eps->eigr[i] = ((b-a)*(eps->eigr[i]+1.0)/2.0+a)*rgscale;
182:           if (a==b) {
183: #if defined(PETSC_USE_COMPLEX)
184:             eps->eigr[i] = ((d-c)*(eps->eigr[i]+1.0)/2.0+c)*rgscale*PETSC_i;
185: #else
186:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
187: #endif
188:           }
189:         }
190:       } else {
191:         center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
192:         radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
193:         for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
194:       }
195:     } else if (isring) {  /* only supported in complex scalars */
196: #if defined(PETSC_USE_COMPLEX)
197:       PetscCall(RGRingGetParameters(eps->rg,&center,&radius,&vscale,&start_ang,&end_ang,NULL));
198:       if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
199:         for (i=0;i<nv;i++) {
200:           theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
201:           eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
202:         }
203:       } else {
204:         for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
205:       }
206: #endif
207:     }
208:   }
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode EPSSetUp_CISS(EPS eps)
213: {
214:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
215:   SlepcContourData contour;
216:   PetscBool        istrivial,isring,isellipse,isinterval,flg;
217:   PetscReal        c,d;
218:   PetscInt         nsplit;
219:   PetscRandom      rand;
220:   PetscObjectId    id;
221:   PetscObjectState state;
222:   Mat              A[2],Psplit[2];
223:   Vec              v0;

225:   PetscFunctionBegin;
226:   if (eps->ncv==PETSC_DEFAULT) {
227:     eps->ncv = ctx->L_max*ctx->M;
228:     if (eps->ncv>eps->n) {
229:       eps->ncv = eps->n;
230:       ctx->L_max = eps->ncv/ctx->M;
231:       PetscCheck(ctx->L_max,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
232:     }
233:   } else {
234:     PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
235:     ctx->L_max = eps->ncv/ctx->M;
236:     if (!ctx->L_max) {
237:       ctx->L_max = 1;
238:       eps->ncv = ctx->L_max*ctx->M;
239:     }
240:   }
241:   ctx->L = PetscMin(ctx->L,ctx->L_max);
242:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 5;
243:   if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
244:   if (!eps->which) eps->which = EPS_ALL;
245:   PetscCheck(eps->which==EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
246:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);

248:   /* check region */
249:   PetscCall(RGIsTrivial(eps->rg,&istrivial));
250:   PetscCheck(!istrivial,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
251:   PetscCall(RGGetComplement(eps->rg,&flg));
252:   PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
253:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
254:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
255:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
256:   PetscCheck(isellipse || isring || isinterval,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");

258:   /* if the region has changed, then reset contour data */
259:   PetscCall(PetscObjectGetId((PetscObject)eps->rg,&id));
260:   PetscCall(PetscObjectStateGet((PetscObject)eps->rg,&state));
261:   if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
262:     PetscCall(SlepcContourDataDestroy(&ctx->contour));
263:     PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of region\n"));
264:     ctx->rgid = id; ctx->rgstate = state;
265:   }

267: #if !defined(PETSC_USE_COMPLEX)
268:   PetscCheck(!isring,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
269: #endif
270:   if (isinterval) {
271:     PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
272: #if !defined(PETSC_USE_COMPLEX)
273:     PetscCheck(c==d && c==0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"In real scalars, endpoints of the imaginary axis must be both zero");
274: #endif
275:     if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
276:   }
277:   if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;

279:   /* create contour data structure */
280:   if (!ctx->contour) {
281:     PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
282:     PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
283:   }

285:   PetscCall(EPSAllocateSolution(eps,0));
286:   PetscCall(BVGetRandomContext(eps->V,&rand));  /* make sure the random context is available when duplicating */
287:   if (ctx->weight) PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
288:   PetscCall(PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma));

290:   /* allocate basis vectors */
291:   PetscCall(BVDestroy(&ctx->S));
292:   PetscCall(BVDuplicateResize(eps->V,ctx->L*ctx->M,&ctx->S));
293:   PetscCall(BVDestroy(&ctx->V));
294:   PetscCall(BVDuplicateResize(eps->V,ctx->L,&ctx->V));

296:   PetscCall(STGetMatrix(eps->st,0,&A[0]));
297:   PetscCall(MatIsShell(A[0],&flg));
298:   PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
299:   if (eps->isgeneralized) PetscCall(STGetMatrix(eps->st,1,&A[1]));

301:   if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
302:   PetscCheck(!ctx->usest || ctx->npart==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");

304:   /* check if a user-defined split preconditioner has been set */
305:   PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
306:   if (nsplit) {
307:     PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Psplit[0]));
308:     if (eps->isgeneralized) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Psplit[1]));
309:   }

311:   contour = ctx->contour;
312:   PetscCall(SlepcContourRedundantMat(contour,eps->isgeneralized?2:1,A,nsplit?Psplit:NULL));
313:   if (contour->pA) {
314:     PetscCall(BVGetColumn(ctx->V,0,&v0));
315:     PetscCall(SlepcContourScatterCreate(contour,v0));
316:     PetscCall(BVRestoreColumn(ctx->V,0,&v0));
317:     PetscCall(BVDestroy(&ctx->pV));
318:     PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV));
319:     PetscCall(BVSetSizesFromVec(ctx->pV,contour->xsub,eps->n));
320:     PetscCall(BVSetFromOptions(ctx->pV));
321:     PetscCall(BVResize(ctx->pV,ctx->L,PETSC_FALSE));
322:   }

324:   EPSCheckDefinite(eps);
325:   EPSCheckSinvertCondition(eps,ctx->usest," (with the usest flag set)");

327:   PetscCall(BVDestroy(&ctx->Y));
328:   if (contour->pA) {
329:     PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y));
330:     PetscCall(BVSetSizesFromVec(ctx->Y,contour->xsub,eps->n));
331:     PetscCall(BVSetFromOptions(ctx->Y));
332:     PetscCall(BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE));
333:   } else PetscCall(BVDuplicateResize(eps->V,contour->npoints*ctx->L,&ctx->Y));

335:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(DSSetType(eps->ds,DSGNHEP));
336:   else if (eps->isgeneralized) {
337:     if (eps->ishermitian && eps->ispositive) PetscCall(DSSetType(eps->ds,DSGHEP));
338:     else PetscCall(DSSetType(eps->ds,DSGNHEP));
339:   } else {
340:     if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
341:     else PetscCall(DSSetType(eps->ds,DSNHEP));
342:   }
343:   PetscCall(DSAllocate(eps->ds,eps->ncv));

345: #if !defined(PETSC_USE_COMPLEX)
346:   PetscCall(EPSSetWorkVecs(eps,3));
347:   if (!eps->ishermitian) PetscCall(PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"));
348: #else
349:   PetscCall(EPSSetWorkVecs(eps,2));
350: #endif
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

354: static PetscErrorCode EPSSetUpSort_CISS(EPS eps)
355: {
356:   SlepcSC        sc;

358:   PetscFunctionBegin;
359:   /* fill sorting criterion context */
360:   eps->sc->comparison    = SlepcCompareSmallestReal;
361:   eps->sc->comparisonctx = NULL;
362:   eps->sc->map           = NULL;
363:   eps->sc->mapobj        = NULL;

365:   /* fill sorting criterion for DS */
366:   PetscCall(DSGetSlepcSC(eps->ds,&sc));
367:   sc->comparison    = SlepcCompareLargestMagnitude;
368:   sc->comparisonctx = NULL;
369:   sc->map           = NULL;
370:   sc->mapobj        = NULL;
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: static PetscErrorCode EPSSolve_CISS(EPS eps)
375: {
376:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
377:   SlepcContourData contour = ctx->contour;
378:   Mat              A,B,X,M,pA,pB,T,J,Pa=NULL,Pb=NULL;
379:   BV               V;
380:   PetscInt         i,j,ld,nmat,L_add=0,nv=0,L_base=ctx->L,inner,nlocal,*inside,nsplit;
381:   PetscScalar      *Mu,*H0,*H1=NULL,*rr,*temp;
382:   PetscReal        error,max_error,norm;
383:   PetscBool        *fl1;
384:   Vec              si,si1=NULL,w[3];
385:   PetscRandom      rand;
386: #if defined(PETSC_USE_COMPLEX)
387:   PetscBool        isellipse;
388:   PetscReal        est_eig,eta;
389: #else
390:   PetscReal        normi;
391: #endif

393:   PetscFunctionBegin;
394:   w[0] = eps->work[0];
395: #if defined(PETSC_USE_COMPLEX)
396:   w[1] = NULL;
397: #else
398:   w[1] = eps->work[2];
399: #endif
400:   w[2] = eps->work[1];
401:   PetscCall(VecGetLocalSize(w[0],&nlocal));
402:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
403:   PetscCall(RGComputeQuadrature(eps->rg,ctx->quad==EPS_CISS_QUADRULE_CHEBYSHEV?RG_QUADRULE_CHEBYSHEV:RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight));
404:   PetscCall(STGetNumMatrices(eps->st,&nmat));
405:   PetscCall(STGetMatrix(eps->st,0,&A));
406:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
407:   else B = NULL;
408:   J = (contour->pA && nmat>1)? contour->pA[1]: B;
409:   V = contour->pA? ctx->pV: ctx->V;
410:   if (!ctx->usest) {
411:     T = contour->pA? contour->pA[0]: A;
412:     PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
413:     if (nsplit) {
414:       if (contour->pA) {
415:         Pa = contour->pP[0];
416:         if (nsplit>1) Pb = contour->pP[1];
417:       } else {
418:         PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Pa));
419:         if (nsplit>1) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Pb));
420:       }
421:     }
422:     PetscCall(EPSCISSSetUp(eps,T,J,Pa,Pb));
423:   }
424:   PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
425:   PetscCall(BVSetRandomSign(ctx->V));
426:   PetscCall(BVGetRandomContext(ctx->V,&rand));

428:   if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
429:   PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
430: #if defined(PETSC_USE_COMPLEX)
431:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
432:   if (isellipse) {
433:     PetscCall(BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig));
434:     PetscCall(PetscInfo(eps,"Estimated eigenvalue count: %f\n",(double)est_eig));
435:     eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
436:     L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
437:     if (L_add>ctx->L_max-ctx->L) {
438:       PetscCall(PetscInfo(eps,"Number of eigenvalues inside the contour path may be too large\n"));
439:       L_add = ctx->L_max-ctx->L;
440:     }
441:   }
442: #endif
443:   if (L_add>0) {
444:     PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add));
445:     PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
446:     PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
447:     PetscCall(BVSetRandomSign(ctx->V));
448:     if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
449:     ctx->L += L_add;
450:     PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
451:   }
452:   PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
453:   for (i=0;i<ctx->refine_blocksize;i++) {
454:     PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
455:     PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
456:     PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
457:     PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
458:     PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
459:     if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
460:     L_add = L_base;
461:     if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
462:     PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add));
463:     PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
464:     PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
465:     PetscCall(BVSetRandomSign(ctx->V));
466:     if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
467:     ctx->L += L_add;
468:     PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
469:     if (L_add) {
470:       PetscCall(PetscFree2(Mu,H0));
471:       PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
472:     }
473:   }
474:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1));

476:   while (eps->reason == EPS_CONVERGED_ITERATING) {
477:     eps->its++;
478:     for (inner=0;inner<=ctx->refine_inner;inner++) {
479:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
480:         PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
481:         PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
482:         PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
483:         PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
484:         PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
485:         break;
486:       } else {
487:         PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
488:         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
489:         PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
490:         PetscCall(BVCopy(ctx->S,ctx->V));
491:         PetscCall(BVSVDAndRank(ctx->S,ctx->M,ctx->L,ctx->delta,BV_SVD_METHOD_REFINE,H0,ctx->sigma,&nv));
492:         if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
493:           if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
494:           PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
495:         } else break;
496:       }
497:     }
498:     eps->nconv = 0;
499:     if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
500:     else {
501:       PetscCall(DSSetDimensions(eps->ds,nv,0,0));
502:       PetscCall(DSSetState(eps->ds,DS_STATE_RAW));

504:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
505:         PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
506:         PetscCall(CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1));
507:         PetscCall(DSGetArray(eps->ds,DS_MAT_A,&temp));
508:         for (j=0;j<nv;j++) {
509:           for (i=0;i<nv;i++) {
510:             temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
511:           }
512:         }
513:         PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&temp));
514:         PetscCall(DSGetArray(eps->ds,DS_MAT_B,&temp));
515:         for (j=0;j<nv;j++) {
516:           for (i=0;i<nv;i++) {
517:             temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
518:           }
519:         }
520:         PetscCall(DSRestoreArray(eps->ds,DS_MAT_B,&temp));
521:       } else {
522:         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
523:         PetscCall(DSGetMat(eps->ds,DS_MAT_A,&pA));
524:         PetscCall(MatZeroEntries(pA));
525:         PetscCall(BVMatProject(ctx->S,A,ctx->S,pA));
526:         PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&pA));
527:         if (B) {
528:           PetscCall(DSGetMat(eps->ds,DS_MAT_B,&pB));
529:           PetscCall(MatZeroEntries(pB));
530:           PetscCall(BVMatProject(ctx->S,B,ctx->S,pB));
531:           PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&pB));
532:         }
533:       }

535:       PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
536:       PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

538:       PetscCall(PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr));
539:       PetscCall(rescale_eig(eps,nv));
540:       PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
541:       PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
542:       PetscCall(SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1));
543:       PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
544:       PetscCall(RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside));
545:       for (i=0;i<nv;i++) {
546:         if (fl1[i] && inside[i]>=0) {
547:           rr[i] = 1.0;
548:           eps->nconv++;
549:         } else rr[i] = 0.0;
550:       }
551:       PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv));
552:       PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
553:       PetscCall(rescale_eig(eps,nv));
554:       PetscCall(PetscFree3(fl1,inside,rr));
555:       PetscCall(BVSetActiveColumns(eps->V,0,nv));
556:       if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
557:         PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
558:         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
559:         PetscCall(BVCopy(ctx->S,ctx->V));
560:         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
561:       }
562:       PetscCall(BVCopy(ctx->S,eps->V));

564:       PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
565:       PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
566:       PetscCall(BVMultInPlace(ctx->S,X,0,eps->nconv));
567:       if (eps->ishermitian) PetscCall(BVMultInPlace(eps->V,X,0,eps->nconv));
568:       PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
569:       max_error = 0.0;
570:       for (i=0;i<eps->nconv;i++) {
571:         PetscCall(BVGetColumn(ctx->S,i,&si));
572: #if !defined(PETSC_USE_COMPLEX)
573:         if (eps->eigi[i]!=0.0) PetscCall(BVGetColumn(ctx->S,i+1,&si1));
574: #endif
575:         PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,si1,w,&error));
576:         if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {  /* vector is not normalized */
577:           PetscCall(VecNorm(si,NORM_2,&norm));
578: #if !defined(PETSC_USE_COMPLEX)
579:           if (eps->eigi[i]!=0.0) {
580:             PetscCall(VecNorm(si1,NORM_2,&normi));
581:             norm = SlepcAbsEigenvalue(norm,normi);
582:           }
583: #endif
584:           error /= norm;
585:         }
586:         PetscCall((*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx));
587:         PetscCall(BVRestoreColumn(ctx->S,i,&si));
588: #if !defined(PETSC_USE_COMPLEX)
589:         if (eps->eigi[i]!=0.0) {
590:           PetscCall(BVRestoreColumn(ctx->S,i+1,&si1));
591:           i++;
592:         }
593: #endif
594:         max_error = PetscMax(max_error,error);
595:       }

597:       if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
598:       else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
599:       else {
600:         if (eps->nconv > ctx->L) nv = eps->nconv;
601:         else if (ctx->L > nv) nv = ctx->L;
602:         nv = PetscMin(nv,ctx->L*ctx->M);
603:         PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M));
604:         PetscCall(MatSetRandom(M,rand));
605:         PetscCall(BVSetActiveColumns(ctx->S,0,nv));
606:         PetscCall(BVMultInPlace(ctx->S,M,0,ctx->L));
607:         PetscCall(MatDestroy(&M));
608:         PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
609:         PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
610:         PetscCall(BVCopy(ctx->S,ctx->V));
611:         if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
612:         PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
613:       }
614:     }
615:   }
616:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscFree(H1));
617:   PetscCall(PetscFree2(Mu,H0));
618:   PetscFunctionReturn(PETSC_SUCCESS);
619: }

621: static PetscErrorCode EPSComputeVectors_CISS(EPS eps)
622: {
623:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
624:   PetscInt       n;
625:   Mat            Z,B=NULL;

627:   PetscFunctionBegin;
628:   if (eps->ishermitian) {
629:     if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
630:     else PetscCall(EPSComputeVectors_Hermitian(eps));
631:     if (eps->isgeneralized && eps->ispositive && ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
632:       /* normalize to have unit B-norm */
633:       PetscCall(STGetMatrix(eps->st,1,&B));
634:       PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
635:       PetscCall(BVNormalize(eps->V,NULL));
636:       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
637:     }
638:     PetscFunctionReturn(PETSC_SUCCESS);
639:   }
640:   PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
641:   PetscCall(BVSetActiveColumns(eps->V,0,n));

643:   /* right eigenvectors */
644:   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));

646:   /* V = V * Z */
647:   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
648:   PetscCall(BVMultInPlace(eps->V,Z,0,n));
649:   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
650:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));

652:   /* normalize */
653:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(BVNormalize(eps->V,NULL));
654:   PetscFunctionReturn(PETSC_SUCCESS);
655: }

657: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
658: {
659:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
660:   PetscInt       oN,oL,oM,oLmax,onpart;
661:   PetscMPIInt    size;

663:   PetscFunctionBegin;
664:   oN = ctx->N;
665:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
666:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
667:   } else {
668:     PetscCheck(ip>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
669:     PetscCheck(ip%2==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
670:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
671:   }
672:   oL = ctx->L;
673:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
674:     ctx->L = 16;
675:   } else {
676:     PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
677:     ctx->L = bs;
678:   }
679:   oM = ctx->M;
680:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
681:     ctx->M = ctx->N/4;
682:   } else {
683:     PetscCheck(ms>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
684:     PetscCheck(ms<=ctx->N,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
685:     ctx->M = ms;
686:   }
687:   onpart = ctx->npart;
688:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
689:     ctx->npart = 1;
690:   } else {
691:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
692:     PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
693:     ctx->npart = npart;
694:   }
695:   oLmax = ctx->L_max;
696:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
697:     ctx->L_max = 64;
698:   } else {
699:     PetscCheck(bsmax>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
700:     ctx->L_max = PetscMax(bsmax,ctx->L);
701:   }
702:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
703:     PetscCall(SlepcContourDataDestroy(&ctx->contour));
704:     PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of parameters\n"));
705:     eps->state = EPS_STATE_INITIAL;
706:   }
707:   ctx->isreal = realmats;
708:   if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) eps->state = EPS_STATE_INITIAL;
709:   PetscFunctionReturn(PETSC_SUCCESS);
710: }

712: /*@
713:    EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.

715:    Logically Collective

717:    Input Parameters:
718: +  eps   - the eigenproblem solver context
719: .  ip    - number of integration points
720: .  bs    - block size
721: .  ms    - moment size
722: .  npart - number of partitions when splitting the communicator
723: .  bsmax - max block size
724: -  realmats - A and B are real

726:    Options Database Keys:
727: +  -eps_ciss_integration_points - Sets the number of integration points
728: .  -eps_ciss_blocksize - Sets the block size
729: .  -eps_ciss_moments - Sets the moment size
730: .  -eps_ciss_partitions - Sets the number of partitions
731: .  -eps_ciss_maxblocksize - Sets the maximum block size
732: -  -eps_ciss_realmats - A and B are real

734:    Note:
735:    The default number of partitions is 1. This means the internal KSP object is shared
736:    among all processes of the EPS communicator. Otherwise, the communicator is split
737:    into npart communicators, so that npart KSP solves proceed simultaneously.

739:    Level: advanced

741: .seealso: EPSCISSGetSizes()
742: @*/
743: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
744: {
745:   PetscFunctionBegin;
753:   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
754:   PetscFunctionReturn(PETSC_SUCCESS);
755: }

757: static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
758: {
759:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

761:   PetscFunctionBegin;
762:   if (ip) *ip = ctx->N;
763:   if (bs) *bs = ctx->L;
764:   if (ms) *ms = ctx->M;
765:   if (npart) *npart = ctx->npart;
766:   if (bsmax) *bsmax = ctx->L_max;
767:   if (realmats) *realmats = ctx->isreal;
768:   PetscFunctionReturn(PETSC_SUCCESS);
769: }

771: /*@
772:    EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.

774:    Not Collective

776:    Input Parameter:
777: .  eps - the eigenproblem solver context

779:    Output Parameters:
780: +  ip    - number of integration points
781: .  bs    - block size
782: .  ms    - moment size
783: .  npart - number of partitions when splitting the communicator
784: .  bsmax - max block size
785: -  realmats - A and B are real

787:    Level: advanced

789: .seealso: EPSCISSSetSizes()
790: @*/
791: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
792: {
793:   PetscFunctionBegin;
795:   PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
796:   PetscFunctionReturn(PETSC_SUCCESS);
797: }

799: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
800: {
801:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

803:   PetscFunctionBegin;
804:   if (delta == (PetscReal)PETSC_DEFAULT) {
805:     ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
806:   } else {
807:     PetscCheck(delta>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
808:     ctx->delta = delta;
809:   }
810:   if (spur == (PetscReal)PETSC_DEFAULT) {
811:     ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
812:   } else {
813:     PetscCheck(spur>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
814:     ctx->spurious_threshold = spur;
815:   }
816:   PetscFunctionReturn(PETSC_SUCCESS);
817: }

819: /*@
820:    EPSCISSSetThreshold - Sets the values of various threshold parameters in
821:    the CISS solver.

823:    Logically Collective

825:    Input Parameters:
826: +  eps   - the eigenproblem solver context
827: .  delta - threshold for numerical rank
828: -  spur  - spurious threshold (to discard spurious eigenpairs)

830:    Options Database Keys:
831: +  -eps_ciss_delta - Sets the delta
832: -  -eps_ciss_spurious_threshold - Sets the spurious threshold

834:    Level: advanced

836: .seealso: EPSCISSGetThreshold()
837: @*/
838: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
839: {
840:   PetscFunctionBegin;
844:   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
845:   PetscFunctionReturn(PETSC_SUCCESS);
846: }

848: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
849: {
850:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

852:   PetscFunctionBegin;
853:   if (delta) *delta = ctx->delta;
854:   if (spur)  *spur = ctx->spurious_threshold;
855:   PetscFunctionReturn(PETSC_SUCCESS);
856: }

858: /*@
859:    EPSCISSGetThreshold - Gets the values of various threshold parameters
860:    in the CISS solver.

862:    Not Collective

864:    Input Parameter:
865: .  eps - the eigenproblem solver context

867:    Output Parameters:
868: +  delta - threshold for numerical rank
869: -  spur  - spurious threshold (to discard spurious eigenpairs)

871:    Level: advanced

873: .seealso: EPSCISSSetThreshold()
874: @*/
875: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
876: {
877:   PetscFunctionBegin;
879:   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
880:   PetscFunctionReturn(PETSC_SUCCESS);
881: }

883: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
884: {
885:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

887:   PetscFunctionBegin;
888:   if (inner == PETSC_DEFAULT) {
889:     ctx->refine_inner = 0;
890:   } else {
891:     PetscCheck(inner>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
892:     ctx->refine_inner = inner;
893:   }
894:   if (blsize == PETSC_DEFAULT) {
895:     ctx->refine_blocksize = 0;
896:   } else {
897:     PetscCheck(blsize>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
898:     ctx->refine_blocksize = blsize;
899:   }
900:   PetscFunctionReturn(PETSC_SUCCESS);
901: }

903: /*@
904:    EPSCISSSetRefinement - Sets the values of various refinement parameters
905:    in the CISS solver.

907:    Logically Collective

909:    Input Parameters:
910: +  eps    - the eigenproblem solver context
911: .  inner  - number of iterative refinement iterations (inner loop)
912: -  blsize - number of iterative refinement iterations (blocksize loop)

914:    Options Database Keys:
915: +  -eps_ciss_refine_inner - Sets number of inner iterations
916: -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations

918:    Level: advanced

920: .seealso: EPSCISSGetRefinement()
921: @*/
922: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
923: {
924:   PetscFunctionBegin;
928:   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
929:   PetscFunctionReturn(PETSC_SUCCESS);
930: }

932: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
933: {
934:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

936:   PetscFunctionBegin;
937:   if (inner)  *inner = ctx->refine_inner;
938:   if (blsize) *blsize = ctx->refine_blocksize;
939:   PetscFunctionReturn(PETSC_SUCCESS);
940: }

942: /*@
943:    EPSCISSGetRefinement - Gets the values of various refinement parameters
944:    in the CISS solver.

946:    Not Collective

948:    Input Parameter:
949: .  eps - the eigenproblem solver context

951:    Output Parameters:
952: +  inner  - number of iterative refinement iterations (inner loop)
953: -  blsize - number of iterative refinement iterations (blocksize loop)

955:    Level: advanced

957: .seealso: EPSCISSSetRefinement()
958: @*/
959: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
960: {
961:   PetscFunctionBegin;
963:   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
964:   PetscFunctionReturn(PETSC_SUCCESS);
965: }

967: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
968: {
969:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

971:   PetscFunctionBegin;
972:   ctx->usest     = usest;
973:   ctx->usest_set = PETSC_TRUE;
974:   eps->state     = EPS_STATE_INITIAL;
975:   PetscFunctionReturn(PETSC_SUCCESS);
976: }

978: /*@
979:    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
980:    use the ST object for the linear solves.

982:    Logically Collective

984:    Input Parameters:
985: +  eps    - the eigenproblem solver context
986: -  usest  - boolean flag to use the ST object or not

988:    Options Database Keys:
989: .  -eps_ciss_usest <bool> - whether the ST object will be used or not

991:    Level: advanced

993: .seealso: EPSCISSGetUseST()
994: @*/
995: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
996: {
997:   PetscFunctionBegin;
1000:   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1001:   PetscFunctionReturn(PETSC_SUCCESS);
1002: }

1004: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1005: {
1006:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1008:   PetscFunctionBegin;
1009:   *usest = ctx->usest;
1010:   PetscFunctionReturn(PETSC_SUCCESS);
1011: }

1013: /*@
1014:    EPSCISSGetUseST - Gets the flag for using the ST object
1015:    in the CISS solver.

1017:    Not Collective

1019:    Input Parameter:
1020: .  eps - the eigenproblem solver context

1022:    Output Parameters:
1023: .  usest - boolean flag indicating if the ST object is being used

1025:    Level: advanced

1027: .seealso: EPSCISSSetUseST()
1028: @*/
1029: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1030: {
1031:   PetscFunctionBegin;
1033:   PetscAssertPointer(usest,2);
1034:   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1035:   PetscFunctionReturn(PETSC_SUCCESS);
1036: }

1038: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1039: {
1040:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1042:   PetscFunctionBegin;
1043:   if (ctx->quad != quad) {
1044:     ctx->quad  = quad;
1045:     eps->state = EPS_STATE_INITIAL;
1046:   }
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: /*@
1051:    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.

1053:    Logically Collective

1055:    Input Parameters:
1056: +  eps  - the eigenproblem solver context
1057: -  quad - the quadrature rule

1059:    Options Database Key:
1060: .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1061:                            'chebyshev')

1063:    Notes:
1064:    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).

1066:    If the 'chebyshev' option is specified (EPS_CISS_QUADRULE_CHEBYSHEV), then
1067:    Chebyshev points are used as quadrature points.

1069:    Level: advanced

1071: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1072: @*/
1073: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1074: {
1075:   PetscFunctionBegin;
1078:   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1079:   PetscFunctionReturn(PETSC_SUCCESS);
1080: }

1082: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1083: {
1084:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1086:   PetscFunctionBegin;
1087:   *quad = ctx->quad;
1088:   PetscFunctionReturn(PETSC_SUCCESS);
1089: }

1091: /*@
1092:    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.

1094:    Not Collective

1096:    Input Parameter:
1097: .  eps - the eigenproblem solver context

1099:    Output Parameters:
1100: .  quad - quadrature rule

1102:    Level: advanced

1104: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1105: @*/
1106: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1107: {
1108:   PetscFunctionBegin;
1110:   PetscAssertPointer(quad,2);
1111:   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1112:   PetscFunctionReturn(PETSC_SUCCESS);
1113: }

1115: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1116: {
1117:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1119:   PetscFunctionBegin;
1120:   if (ctx->extraction != extraction) {
1121:     ctx->extraction = extraction;
1122:     eps->state      = EPS_STATE_INITIAL;
1123:   }
1124:   PetscFunctionReturn(PETSC_SUCCESS);
1125: }

1127: /*@
1128:    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.

1130:    Logically Collective

1132:    Input Parameters:
1133: +  eps        - the eigenproblem solver context
1134: -  extraction - the extraction technique

1136:    Options Database Key:
1137: .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1138:                            'hankel')

1140:    Notes:
1141:    By default, the Rayleigh-Ritz extraction is used (EPS_CISS_EXTRACTION_RITZ).

1143:    If the 'hankel' option is specified (EPS_CISS_EXTRACTION_HANKEL), then
1144:    the Block Hankel method is used for extracting eigenpairs.

1146:    Level: advanced

1148: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1149: @*/
1150: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1151: {
1152:   PetscFunctionBegin;
1155:   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1156:   PetscFunctionReturn(PETSC_SUCCESS);
1157: }

1159: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1160: {
1161:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1163:   PetscFunctionBegin;
1164:   *extraction = ctx->extraction;
1165:   PetscFunctionReturn(PETSC_SUCCESS);
1166: }

1168: /*@
1169:    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.

1171:    Not Collective

1173:    Input Parameter:
1174: .  eps - the eigenproblem solver context

1176:    Output Parameters:
1177: .  extraction - extraction technique

1179:    Level: advanced

1181: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1182: @*/
1183: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1184: {
1185:   PetscFunctionBegin;
1187:   PetscAssertPointer(extraction,2);
1188:   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1189:   PetscFunctionReturn(PETSC_SUCCESS);
1190: }

1192: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1193: {
1194:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
1195:   SlepcContourData contour;
1196:   PetscInt         i,nsplit;
1197:   PC               pc;
1198:   MPI_Comm         child;

1200:   PetscFunctionBegin;
1201:   if (!ctx->contour) {  /* initialize contour data structure first */
1202:     PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
1203:     PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
1204:   }
1205:   contour = ctx->contour;
1206:   if (!contour->ksp) {
1207:     PetscCall(PetscMalloc1(contour->npoints,&contour->ksp));
1208:     PetscCall(EPSGetST(eps,&eps->st));
1209:     PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
1210:     PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
1211:     for (i=0;i<contour->npoints;i++) {
1212:       PetscCall(KSPCreate(child,&contour->ksp[i]));
1213:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1));
1214:       PetscCall(KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix));
1215:       PetscCall(KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_"));
1216:       PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options));
1217:       PetscCall(KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE));
1218:       PetscCall(KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
1219:       PetscCall(KSPGetPC(contour->ksp[i],&pc));
1220:       if (nsplit) {
1221:         PetscCall(KSPSetType(contour->ksp[i],KSPBCGS));
1222:         PetscCall(PCSetType(pc,PCBJACOBI));
1223:       } else {
1224:         PetscCall(KSPSetType(contour->ksp[i],KSPPREONLY));
1225:         PetscCall(PCSetType(pc,PCLU));
1226:       }
1227:     }
1228:   }
1229:   if (nsolve) *nsolve = contour->npoints;
1230:   if (ksp)    *ksp    = contour->ksp;
1231:   PetscFunctionReturn(PETSC_SUCCESS);
1232: }

1234: /*@C
1235:    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1236:    the CISS solver.

1238:    Not Collective

1240:    Input Parameter:
1241: .  eps - the eigenproblem solver solver

1243:    Output Parameters:
1244: +  nsolve - number of solver objects
1245: -  ksp - array of linear solver object

1247:    Notes:
1248:    The number of KSP solvers is equal to the number of integration points divided by
1249:    the number of partitions. This value is halved in the case of real matrices with
1250:    a region centered at the real axis.

1252:    Level: advanced

1254: .seealso: EPSCISSSetSizes()
1255: @*/
1256: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1257: {
1258:   PetscFunctionBegin;
1260:   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1261:   PetscFunctionReturn(PETSC_SUCCESS);
1262: }

1264: static PetscErrorCode EPSReset_CISS(EPS eps)
1265: {
1266:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1268:   PetscFunctionBegin;
1269:   PetscCall(BVDestroy(&ctx->S));
1270:   PetscCall(BVDestroy(&ctx->V));
1271:   PetscCall(BVDestroy(&ctx->Y));
1272:   if (!ctx->usest) PetscCall(SlepcContourDataReset(ctx->contour));
1273:   PetscCall(BVDestroy(&ctx->pV));
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }

1277: static PetscErrorCode EPSSetFromOptions_CISS(EPS eps,PetscOptionItems *PetscOptionsObject)
1278: {
1279:   PetscReal         r3,r4;
1280:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
1281:   PetscBool         b1,b2,flg,flg2,flg3,flg4,flg5,flg6;
1282:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
1283:   EPSCISSQuadRule   quad;
1284:   EPSCISSExtraction extraction;

1286:   PetscFunctionBegin;
1287:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS CISS Options");

1289:     PetscCall(EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1));
1290:     PetscCall(PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg));
1291:     PetscCall(PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2));
1292:     PetscCall(PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3));
1293:     PetscCall(PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4));
1294:     PetscCall(PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5));
1295:     PetscCall(PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6));
1296:     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1));

1298:     PetscCall(EPSCISSGetThreshold(eps,&r3,&r4));
1299:     PetscCall(PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg));
1300:     PetscCall(PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2));
1301:     if (flg || flg2) PetscCall(EPSCISSSetThreshold(eps,r3,r4));

1303:     PetscCall(EPSCISSGetRefinement(eps,&i6,&i7));
1304:     PetscCall(PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg));
1305:     PetscCall(PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2));
1306:     if (flg || flg2) PetscCall(EPSCISSSetRefinement(eps,i6,i7));

1308:     PetscCall(EPSCISSGetUseST(eps,&b2));
1309:     PetscCall(PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg));
1310:     if (flg) PetscCall(EPSCISSSetUseST(eps,b2));

1312:     PetscCall(PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg));
1313:     if (flg) PetscCall(EPSCISSSetQuadRule(eps,quad));

1315:     PetscCall(PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg));
1316:     if (flg) PetscCall(EPSCISSSetExtraction(eps,extraction));

1318:   PetscOptionsHeadEnd();

1320:   if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
1321:   PetscCall(RGSetFromOptions(eps->rg)); /* this is necessary here to set useconj */
1322:   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
1323:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1324:   for (i=0;i<ctx->contour->npoints;i++) PetscCall(KSPSetFromOptions(ctx->contour->ksp[i]));
1325:   PetscCall(PetscSubcommSetFromOptions(ctx->contour->subcomm));
1326:   PetscFunctionReturn(PETSC_SUCCESS);
1327: }

1329: static PetscErrorCode EPSDestroy_CISS(EPS eps)
1330: {
1331:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1333:   PetscFunctionBegin;
1334:   PetscCall(SlepcContourDataDestroy(&ctx->contour));
1335:   PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
1336:   PetscCall(PetscFree(eps->data));
1337:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL));
1338:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL));
1339:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL));
1340:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL));
1341:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL));
1342:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL));
1343:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL));
1344:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL));
1345:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL));
1346:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL));
1347:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL));
1348:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL));
1349:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL));
1350:   PetscFunctionReturn(PETSC_SUCCESS);
1351: }

1353: static PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1354: {
1355:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1356:   PetscBool      isascii;
1357:   PetscViewer    sviewer;

1359:   PetscFunctionBegin;
1360:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1361:   if (isascii) {
1362:     PetscCall(PetscViewerASCIIPrintf(viewer,"  sizes { integration points: %" PetscInt_FMT ", block size: %" PetscInt_FMT ", moment size: %" PetscInt_FMT ", partitions: %" PetscInt_FMT ", maximum block size: %" PetscInt_FMT " }\n",ctx->N,ctx->L,ctx->M,ctx->npart,ctx->L_max));
1363:     if (ctx->isreal) PetscCall(PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n"));
1364:     PetscCall(PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold));
1365:     PetscCall(PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize));
1366:     PetscCall(PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]));
1367:     PetscCall(PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]));
1368:     if (ctx->usest) PetscCall(PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n"));
1369:     else {
1370:       if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
1371:       PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1372:       PetscCall(PetscViewerASCIIPushTab(viewer));
1373:       if (ctx->npart>1 && ctx->contour->subcomm) {
1374:         PetscCall(PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1375:         if (!ctx->contour->subcomm->color) PetscCall(KSPView(ctx->contour->ksp[0],sviewer));
1376:         PetscCall(PetscViewerFlush(sviewer));
1377:         PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1378:         PetscCall(PetscViewerFlush(viewer));
1379:         /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1380:         PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1381:       } else PetscCall(KSPView(ctx->contour->ksp[0],viewer));
1382:       PetscCall(PetscViewerASCIIPopTab(viewer));
1383:     }
1384:   }
1385:   PetscFunctionReturn(PETSC_SUCCESS);
1386: }

1388: static PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1389: {
1390:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1391:   PetscBool      usest = ctx->usest;
1392:   KSP            ksp;
1393:   PC             pc;

1395:   PetscFunctionBegin;
1396:   if (!((PetscObject)eps->st)->type_name) {
1397:     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1398:     if (usest) PetscCall(STSetType(eps->st,STSINVERT));
1399:     else {
1400:       /* we are not going to use ST, so avoid factorizing the matrix */
1401:       PetscCall(STSetType(eps->st,STSHIFT));
1402:       if (eps->isgeneralized) {
1403:         PetscCall(STGetKSP(eps->st,&ksp));
1404:         PetscCall(KSPGetPC(ksp,&pc));
1405:         PetscCall(PCSetType(pc,PCNONE));
1406:       }
1407:     }
1408:   }
1409:   PetscFunctionReturn(PETSC_SUCCESS);
1410: }

1412: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1413: {
1414:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1416:   PetscFunctionBegin;
1417:   PetscCall(PetscNew(&ctx));
1418:   eps->data = ctx;

1420:   eps->useds = PETSC_TRUE;
1421:   eps->categ = EPS_CATEGORY_CONTOUR;

1423:   eps->ops->solve          = EPSSolve_CISS;
1424:   eps->ops->setup          = EPSSetUp_CISS;
1425:   eps->ops->setupsort      = EPSSetUpSort_CISS;
1426:   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1427:   eps->ops->destroy        = EPSDestroy_CISS;
1428:   eps->ops->reset          = EPSReset_CISS;
1429:   eps->ops->view           = EPSView_CISS;
1430:   eps->ops->computevectors = EPSComputeVectors_CISS;
1431:   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;

1433:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS));
1434:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS));
1435:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS));
1436:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS));
1437:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS));
1438:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS));
1439:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS));
1440:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS));
1441:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS));
1442:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS));
1443:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS));
1444:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS));
1445:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS));

1447:   /* set default values of parameters */
1448:   ctx->N                  = 32;
1449:   ctx->L                  = 16;
1450:   ctx->M                  = ctx->N/4;
1451:   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
1452:   ctx->L_max              = 64;
1453:   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1454:   ctx->usest              = PETSC_TRUE;
1455:   ctx->usest_set          = PETSC_FALSE;
1456:   ctx->isreal             = PETSC_FALSE;
1457:   ctx->refine_inner       = 0;
1458:   ctx->refine_blocksize   = 0;
1459:   ctx->npart              = 1;
1460:   ctx->quad               = (EPSCISSQuadRule)0;
1461:   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
1462:   PetscFunctionReturn(PETSC_SUCCESS);
1463: }