Actual source code: ciss.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

716:    Logically Collective

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

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

735:    Note:
736:    For all integer arguments, you can use PETSC_CURRENT to keep the current value, and
737:    PETSC_DETERMINE to set them to a default value.

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

743:    Level: advanced

745: .seealso: EPSCISSGetSizes()
746: @*/
747: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
748: {
749:   PetscFunctionBegin;
757:   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
758:   PetscFunctionReturn(PETSC_SUCCESS);
759: }

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

765:   PetscFunctionBegin;
766:   if (ip) *ip = ctx->N;
767:   if (bs) *bs = ctx->L;
768:   if (ms) *ms = ctx->M;
769:   if (npart) *npart = ctx->npart;
770:   if (bsmax) *bsmax = ctx->L_max;
771:   if (realmats) *realmats = ctx->isreal;
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

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

778:    Not Collective

780:    Input Parameter:
781: .  eps - the eigenproblem solver context

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

791:    Level: advanced

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

803: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
804: {
805:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

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

823: /*@
824:    EPSCISSSetThreshold - Sets the values of various threshold parameters in
825:    the CISS solver.

827:    Logically Collective

829:    Input Parameters:
830: +  eps   - the eigenproblem solver context
831: .  delta - threshold for numerical rank
832: -  spur  - spurious threshold (to discard spurious eigenpairs)

834:    Options Database Keys:
835: +  -eps_ciss_delta - Sets the delta
836: -  -eps_ciss_spurious_threshold - Sets the spurious threshold

838:    Note:
839:    PETSC_CURRENT can be used to preserve the current value of any of the
840:    arguments, and PETSC_DETERMINE to set them to a default value.

842:    Level: advanced

844: .seealso: EPSCISSGetThreshold()
845: @*/
846: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
847: {
848:   PetscFunctionBegin;
852:   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
853:   PetscFunctionReturn(PETSC_SUCCESS);
854: }

856: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
857: {
858:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

860:   PetscFunctionBegin;
861:   if (delta) *delta = ctx->delta;
862:   if (spur)  *spur = ctx->spurious_threshold;
863:   PetscFunctionReturn(PETSC_SUCCESS);
864: }

866: /*@
867:    EPSCISSGetThreshold - Gets the values of various threshold parameters
868:    in the CISS solver.

870:    Not Collective

872:    Input Parameter:
873: .  eps - the eigenproblem solver context

875:    Output Parameters:
876: +  delta - threshold for numerical rank
877: -  spur  - spurious threshold (to discard spurious eigenpairs)

879:    Level: advanced

881: .seealso: EPSCISSSetThreshold()
882: @*/
883: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
884: {
885:   PetscFunctionBegin;
887:   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
888:   PetscFunctionReturn(PETSC_SUCCESS);
889: }

891: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
892: {
893:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

895:   PetscFunctionBegin;
896:   if (inner == PETSC_DETERMINE) {
897:     ctx->refine_inner = 0;
898:   } else if (inner != PETSC_CURRENT) {
899:     PetscCheck(inner>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
900:     ctx->refine_inner = inner;
901:   }
902:   if (blsize == PETSC_DETERMINE) {
903:     ctx->refine_blocksize = 0;
904:   } else if (blsize != PETSC_CURRENT) {
905:     PetscCheck(blsize>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
906:     ctx->refine_blocksize = blsize;
907:   }
908:   PetscFunctionReturn(PETSC_SUCCESS);
909: }

911: /*@
912:    EPSCISSSetRefinement - Sets the values of various refinement parameters
913:    in the CISS solver.

915:    Logically Collective

917:    Input Parameters:
918: +  eps    - the eigenproblem solver context
919: .  inner  - number of iterative refinement iterations (inner loop)
920: -  blsize - number of iterative refinement iterations (blocksize loop)

922:    Options Database Keys:
923: +  -eps_ciss_refine_inner - Sets number of inner iterations
924: -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations

926:    Note:
927:    PETSC_CURRENT can be used to preserve the current value of any of the
928:    arguments, and PETSC_DETERMINE to set them to a default of 0.

930:    Level: advanced

932: .seealso: EPSCISSGetRefinement()
933: @*/
934: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
935: {
936:   PetscFunctionBegin;
940:   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
941:   PetscFunctionReturn(PETSC_SUCCESS);
942: }

944: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
945: {
946:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

948:   PetscFunctionBegin;
949:   if (inner)  *inner = ctx->refine_inner;
950:   if (blsize) *blsize = ctx->refine_blocksize;
951:   PetscFunctionReturn(PETSC_SUCCESS);
952: }

954: /*@
955:    EPSCISSGetRefinement - Gets the values of various refinement parameters
956:    in the CISS solver.

958:    Not Collective

960:    Input Parameter:
961: .  eps - the eigenproblem solver context

963:    Output Parameters:
964: +  inner  - number of iterative refinement iterations (inner loop)
965: -  blsize - number of iterative refinement iterations (blocksize loop)

967:    Level: advanced

969: .seealso: EPSCISSSetRefinement()
970: @*/
971: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
972: {
973:   PetscFunctionBegin;
975:   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
976:   PetscFunctionReturn(PETSC_SUCCESS);
977: }

979: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
980: {
981:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

983:   PetscFunctionBegin;
984:   ctx->usest     = usest;
985:   ctx->usest_set = PETSC_TRUE;
986:   eps->state     = EPS_STATE_INITIAL;
987:   PetscFunctionReturn(PETSC_SUCCESS);
988: }

990: /*@
991:    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
992:    use the ST object for the linear solves.

994:    Logically Collective

996:    Input Parameters:
997: +  eps    - the eigenproblem solver context
998: -  usest  - boolean flag to use the ST object or not

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

1003:    Level: advanced

1005: .seealso: EPSCISSGetUseST()
1006: @*/
1007: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
1008: {
1009:   PetscFunctionBegin;
1012:   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1013:   PetscFunctionReturn(PETSC_SUCCESS);
1014: }

1016: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1017: {
1018:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1020:   PetscFunctionBegin;
1021:   *usest = ctx->usest;
1022:   PetscFunctionReturn(PETSC_SUCCESS);
1023: }

1025: /*@
1026:    EPSCISSGetUseST - Gets the flag for using the ST object
1027:    in the CISS solver.

1029:    Not Collective

1031:    Input Parameter:
1032: .  eps - the eigenproblem solver context

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

1037:    Level: advanced

1039: .seealso: EPSCISSSetUseST()
1040: @*/
1041: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1042: {
1043:   PetscFunctionBegin;
1045:   PetscAssertPointer(usest,2);
1046:   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1051: {
1052:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1054:   PetscFunctionBegin;
1055:   if (ctx->quad != quad) {
1056:     ctx->quad  = quad;
1057:     eps->state = EPS_STATE_INITIAL;
1058:   }
1059:   PetscFunctionReturn(PETSC_SUCCESS);
1060: }

1062: /*@
1063:    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.

1065:    Logically Collective

1067:    Input Parameters:
1068: +  eps  - the eigenproblem solver context
1069: -  quad - the quadrature rule

1071:    Options Database Key:
1072: .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1073:                            'chebyshev')

1075:    Notes:
1076:    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).

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

1081:    Level: advanced

1083: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1084: @*/
1085: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1086: {
1087:   PetscFunctionBegin;
1090:   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

1094: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1095: {
1096:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1098:   PetscFunctionBegin;
1099:   *quad = ctx->quad;
1100:   PetscFunctionReturn(PETSC_SUCCESS);
1101: }

1103: /*@
1104:    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.

1106:    Not Collective

1108:    Input Parameter:
1109: .  eps - the eigenproblem solver context

1111:    Output Parameters:
1112: .  quad - quadrature rule

1114:    Level: advanced

1116: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1117: @*/
1118: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1119: {
1120:   PetscFunctionBegin;
1122:   PetscAssertPointer(quad,2);
1123:   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1124:   PetscFunctionReturn(PETSC_SUCCESS);
1125: }

1127: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1128: {
1129:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1131:   PetscFunctionBegin;
1132:   if (ctx->extraction != extraction) {
1133:     ctx->extraction = extraction;
1134:     eps->state      = EPS_STATE_INITIAL;
1135:   }
1136:   PetscFunctionReturn(PETSC_SUCCESS);
1137: }

1139: /*@
1140:    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.

1142:    Logically Collective

1144:    Input Parameters:
1145: +  eps        - the eigenproblem solver context
1146: -  extraction - the extraction technique

1148:    Options Database Key:
1149: .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1150:                            'hankel')

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

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

1158:    Level: advanced

1160: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1161: @*/
1162: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1163: {
1164:   PetscFunctionBegin;
1167:   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1168:   PetscFunctionReturn(PETSC_SUCCESS);
1169: }

1171: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1172: {
1173:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1175:   PetscFunctionBegin;
1176:   *extraction = ctx->extraction;
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: /*@
1181:    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.

1183:    Not Collective

1185:    Input Parameter:
1186: .  eps - the eigenproblem solver context

1188:    Output Parameters:
1189: .  extraction - extraction technique

1191:    Level: advanced

1193: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1194: @*/
1195: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1196: {
1197:   PetscFunctionBegin;
1199:   PetscAssertPointer(extraction,2);
1200:   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1201:   PetscFunctionReturn(PETSC_SUCCESS);
1202: }

1204: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1205: {
1206:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
1207:   SlepcContourData contour;
1208:   PetscInt         i,nsplit;
1209:   PC               pc;
1210:   MPI_Comm         child;

1212:   PetscFunctionBegin;
1213:   if (!ctx->contour) {  /* initialize contour data structure first */
1214:     PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
1215:     PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
1216:   }
1217:   contour = ctx->contour;
1218:   if (!contour->ksp) {
1219:     PetscCall(PetscMalloc1(contour->npoints,&contour->ksp));
1220:     PetscCall(EPSGetST(eps,&eps->st));
1221:     PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
1222:     PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
1223:     for (i=0;i<contour->npoints;i++) {
1224:       PetscCall(KSPCreate(child,&contour->ksp[i]));
1225:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1));
1226:       PetscCall(KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix));
1227:       PetscCall(KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_"));
1228:       PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options));
1229:       PetscCall(KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE));
1230:       PetscCall(KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1231:       PetscCall(KSPGetPC(contour->ksp[i],&pc));
1232:       if (nsplit) {
1233:         PetscCall(KSPSetType(contour->ksp[i],KSPBCGS));
1234:         PetscCall(PCSetType(pc,PCBJACOBI));
1235:       } else {
1236:         PetscCall(KSPSetType(contour->ksp[i],KSPPREONLY));
1237:         PetscCall(PCSetType(pc,PCLU));
1238:       }
1239:     }
1240:   }
1241:   if (nsolve) *nsolve = contour->npoints;
1242:   if (ksp)    *ksp    = contour->ksp;
1243:   PetscFunctionReturn(PETSC_SUCCESS);
1244: }

1246: /*@C
1247:    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1248:    the CISS solver.

1250:    Not Collective

1252:    Input Parameter:
1253: .  eps - the eigenproblem solver solver

1255:    Output Parameters:
1256: +  nsolve - number of solver objects
1257: -  ksp - array of linear solver object

1259:    Notes:
1260:    The number of KSP solvers is equal to the number of integration points divided by
1261:    the number of partitions. This value is halved in the case of real matrices with
1262:    a region centered at the real axis.

1264:    Level: advanced

1266: .seealso: EPSCISSSetSizes()
1267: @*/
1268: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1269: {
1270:   PetscFunctionBegin;
1272:   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1273:   PetscFunctionReturn(PETSC_SUCCESS);
1274: }

1276: static PetscErrorCode EPSReset_CISS(EPS eps)
1277: {
1278:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1280:   PetscFunctionBegin;
1281:   PetscCall(BVDestroy(&ctx->S));
1282:   PetscCall(BVDestroy(&ctx->V));
1283:   PetscCall(BVDestroy(&ctx->Y));
1284:   if (!ctx->usest) PetscCall(SlepcContourDataReset(ctx->contour));
1285:   PetscCall(BVDestroy(&ctx->pV));
1286:   PetscFunctionReturn(PETSC_SUCCESS);
1287: }

1289: static PetscErrorCode EPSSetFromOptions_CISS(EPS eps,PetscOptionItems *PetscOptionsObject)
1290: {
1291:   PetscReal         r3,r4;
1292:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
1293:   PetscBool         b1,b2,flg,flg2,flg3,flg4,flg5,flg6;
1294:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
1295:   EPSCISSQuadRule   quad;
1296:   EPSCISSExtraction extraction;

1298:   PetscFunctionBegin;
1299:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS CISS Options");

1301:     PetscCall(EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1));
1302:     PetscCall(PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg));
1303:     PetscCall(PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2));
1304:     PetscCall(PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3));
1305:     PetscCall(PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4));
1306:     PetscCall(PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5));
1307:     PetscCall(PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6));
1308:     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1));

1310:     PetscCall(EPSCISSGetThreshold(eps,&r3,&r4));
1311:     PetscCall(PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg));
1312:     PetscCall(PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2));
1313:     if (flg || flg2) PetscCall(EPSCISSSetThreshold(eps,r3,r4));

1315:     PetscCall(EPSCISSGetRefinement(eps,&i6,&i7));
1316:     PetscCall(PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg));
1317:     PetscCall(PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2));
1318:     if (flg || flg2) PetscCall(EPSCISSSetRefinement(eps,i6,i7));

1320:     PetscCall(EPSCISSGetUseST(eps,&b2));
1321:     PetscCall(PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg));
1322:     if (flg) PetscCall(EPSCISSSetUseST(eps,b2));

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

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

1330:   PetscOptionsHeadEnd();

1332:   if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
1333:   PetscCall(RGSetFromOptions(eps->rg)); /* this is necessary here to set useconj */
1334:   if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
1335:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1336:   for (i=0;i<ctx->contour->npoints;i++) PetscCall(KSPSetFromOptions(ctx->contour->ksp[i]));
1337:   PetscCall(PetscSubcommSetFromOptions(ctx->contour->subcomm));
1338:   PetscFunctionReturn(PETSC_SUCCESS);
1339: }

1341: static PetscErrorCode EPSDestroy_CISS(EPS eps)
1342: {
1343:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1345:   PetscFunctionBegin;
1346:   PetscCall(SlepcContourDataDestroy(&ctx->contour));
1347:   PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
1348:   PetscCall(PetscFree(eps->data));
1349:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL));
1350:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL));
1351:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL));
1352:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL));
1353:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL));
1354:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL));
1355:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL));
1356:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL));
1357:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL));
1358:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL));
1359:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL));
1360:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL));
1361:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL));
1362:   PetscFunctionReturn(PETSC_SUCCESS);
1363: }

1365: static PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1366: {
1367:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1368:   PetscBool      isascii;
1369:   PetscViewer    sviewer;

1371:   PetscFunctionBegin;
1372:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1373:   if (isascii) {
1374:     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));
1375:     if (ctx->isreal) PetscCall(PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n"));
1376:     PetscCall(PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold));
1377:     PetscCall(PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize));
1378:     PetscCall(PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]));
1379:     PetscCall(PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]));
1380:     if (ctx->usest) PetscCall(PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n"));
1381:     else {
1382:       if (!ctx->contour || !ctx->contour->ksp) PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
1383:       PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1384:       PetscCall(PetscViewerASCIIPushTab(viewer));
1385:       if (ctx->npart>1 && ctx->contour->subcomm) {
1386:         PetscCall(PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1387:         if (!ctx->contour->subcomm->color) PetscCall(KSPView(ctx->contour->ksp[0],sviewer));
1388:         PetscCall(PetscViewerFlush(sviewer));
1389:         PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1390:         /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1391:         PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1392:       } else PetscCall(KSPView(ctx->contour->ksp[0],viewer));
1393:       PetscCall(PetscViewerASCIIPopTab(viewer));
1394:     }
1395:   }
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

1399: static PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1400: {
1401:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1402:   PetscBool      usest = ctx->usest;
1403:   KSP            ksp;
1404:   PC             pc;

1406:   PetscFunctionBegin;
1407:   if (!((PetscObject)eps->st)->type_name) {
1408:     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1409:     if (usest) PetscCall(STSetType(eps->st,STSINVERT));
1410:     else {
1411:       /* we are not going to use ST, so avoid factorizing the matrix */
1412:       PetscCall(STSetType(eps->st,STSHIFT));
1413:       if (eps->isgeneralized) {
1414:         PetscCall(STGetKSP(eps->st,&ksp));
1415:         PetscCall(KSPGetPC(ksp,&pc));
1416:         PetscCall(PCSetType(pc,PCNONE));
1417:       }
1418:     }
1419:   }
1420:   PetscFunctionReturn(PETSC_SUCCESS);
1421: }

1423: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1424: {
1425:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1427:   PetscFunctionBegin;
1428:   PetscCall(PetscNew(&ctx));
1429:   eps->data = ctx;

1431:   eps->useds = PETSC_TRUE;
1432:   eps->categ = EPS_CATEGORY_CONTOUR;

1434:   eps->ops->solve          = EPSSolve_CISS;
1435:   eps->ops->setup          = EPSSetUp_CISS;
1436:   eps->ops->setupsort      = EPSSetUpSort_CISS;
1437:   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1438:   eps->ops->destroy        = EPSDestroy_CISS;
1439:   eps->ops->reset          = EPSReset_CISS;
1440:   eps->ops->view           = EPSView_CISS;
1441:   eps->ops->computevectors = EPSComputeVectors_CISS;
1442:   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;

1444:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS));
1445:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS));
1446:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS));
1447:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS));
1448:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS));
1449:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS));
1450:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS));
1451:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS));
1452:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS));
1453:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS));
1454:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS));
1455:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS));
1456:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS));

1458:   /* set default values of parameters */
1459:   ctx->N                  = 32;
1460:   ctx->L                  = 16;
1461:   ctx->M                  = ctx->N/4;
1462:   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
1463:   ctx->L_max              = 64;
1464:   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1465:   ctx->usest              = PETSC_TRUE;
1466:   ctx->usest_set          = PETSC_FALSE;
1467:   ctx->isreal             = PETSC_FALSE;
1468:   ctx->refine_inner       = 0;
1469:   ctx->refine_blocksize   = 0;
1470:   ctx->npart              = 1;
1471:   ctx->quad               = (EPSCISSQuadRule)0;
1472:   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
1473:   PetscFunctionReturn(PETSC_SUCCESS);
1474: }