Actual source code: ciss.c

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

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

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

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

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

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

222:   if (eps->ncv==PETSC_DEFAULT) {
223:     eps->ncv = ctx->L_max*ctx->M;
224:     if (eps->ncv>eps->n) {
225:       eps->ncv = eps->n;
226:       ctx->L_max = eps->ncv/ctx->M;
228:     }
229:   } else {
230:     EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
231:     ctx->L_max = eps->ncv/ctx->M;
232:     if (!ctx->L_max) {
233:       ctx->L_max = 1;
234:       eps->ncv = ctx->L_max*ctx->M;
235:     }
236:   }
237:   ctx->L = PetscMin(ctx->L,ctx->L_max);
238:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 5;
239:   if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
240:   if (!eps->which) eps->which = EPS_ALL;
242:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);

244:   /* check region */
245:   RGIsTrivial(eps->rg,&istrivial);
247:   RGGetComplement(eps->rg,&flg);
249:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
250:   PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
251:   PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);

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

263: #if !defined(PETSC_USE_COMPLEX)
265: #endif
266:   if (isinterval) {
267:     RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
268: #if !defined(PETSC_USE_COMPLEX)
270: #endif
271:     if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
272:   }
273:   if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;

275:   /* create contour data structure */
276:   if (!ctx->contour) {
277:     RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj);
278:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour);
279:   }

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

286:   /* allocate basis vectors */
287:   BVDestroy(&ctx->S);
288:   BVDuplicateResize(eps->V,ctx->L*ctx->M,&ctx->S);
289:   BVDestroy(&ctx->V);
290:   BVDuplicateResize(eps->V,ctx->L,&ctx->V);

292:   STGetMatrix(eps->st,0,&A[0]);
293:   MatIsShell(A[0],&flg);
295:   if (eps->isgeneralized) STGetMatrix(eps->st,1,&A[1]);

297:   if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;

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

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

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

323:   BVDestroy(&ctx->Y);
324:   if (contour->pA) {
325:     BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y);
326:     BVSetSizesFromVec(ctx->Y,contour->xsub,eps->n);
327:     BVSetFromOptions(ctx->Y);
328:     BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE);
329:   } else BVDuplicateResize(eps->V,contour->npoints*ctx->L,&ctx->Y);

331:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) DSSetType(eps->ds,DSGNHEP);
332:   else if (eps->isgeneralized) {
333:     if (eps->ishermitian && eps->ispositive) DSSetType(eps->ds,DSGHEP);
334:     else DSSetType(eps->ds,DSGNHEP);
335:   } else {
336:     if (eps->ishermitian) DSSetType(eps->ds,DSHEP);
337:     else DSSetType(eps->ds,DSNHEP);
338:   }
339:   DSAllocate(eps->ds,eps->ncv);

341: #if !defined(PETSC_USE_COMPLEX)
342:   EPSSetWorkVecs(eps,3);
343:   if (!eps->ishermitian) PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n");
344: #else
345:   EPSSetWorkVecs(eps,2);
346: #endif
347:   return 0;
348: }

350: PetscErrorCode EPSSetUpSort_CISS(EPS eps)
351: {
352:   SlepcSC        sc;

354:   /* fill sorting criterion context */
355:   eps->sc->comparison    = SlepcCompareSmallestReal;
356:   eps->sc->comparisonctx = NULL;
357:   eps->sc->map           = NULL;
358:   eps->sc->mapobj        = NULL;

360:   /* fill sorting criterion for DS */
361:   DSGetSlepcSC(eps->ds,&sc);
362:   sc->comparison    = SlepcCompareLargestMagnitude;
363:   sc->comparisonctx = NULL;
364:   sc->map           = NULL;
365:   sc->mapobj        = NULL;
366:   return 0;
367: }

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

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

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

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

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

529:       DSSolve(eps->ds,eps->eigr,eps->eigi);
530:       DSSynchronize(eps->ds,eps->eigr,eps->eigi);

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

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

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

615: PetscErrorCode EPSComputeVectors_CISS(EPS eps)
616: {
617:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
618:   PetscInt       n;
619:   Mat            Z,B=NULL;

621:   if (eps->ishermitian) {
622:     if (eps->isgeneralized && !eps->ispositive) EPSComputeVectors_Indefinite(eps);
623:     else EPSComputeVectors_Hermitian(eps);
624:     if (eps->isgeneralized && eps->ispositive && ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
625:       /* normalize to have unit B-norm */
626:       STGetMatrix(eps->st,1,&B);
627:       BVSetMatrix(eps->V,B,PETSC_FALSE);
628:       BVNormalize(eps->V,NULL);
629:       BVSetMatrix(eps->V,NULL,PETSC_FALSE);
630:     }
631:     return 0;
632:   }
633:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
634:   BVSetActiveColumns(eps->V,0,n);

636:   /* right eigenvectors */
637:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

639:   /* V = V * Z */
640:   DSGetMat(eps->ds,DS_MAT_X,&Z);
641:   BVMultInPlace(eps->V,Z,0,n);
642:   DSRestoreMat(eps->ds,DS_MAT_X,&Z);
643:   BVSetActiveColumns(eps->V,0,eps->nconv);

645:   /* normalize */
646:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) BVNormalize(eps->V,NULL);
647:   return 0;
648: }

650: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
651: {
652:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
653:   PetscInt       oN,oL,oM,oLmax,onpart;
654:   PetscMPIInt    size;

656:   oN = ctx->N;
657:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
658:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
659:   } else {
662:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
663:   }
664:   oL = ctx->L;
665:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
666:     ctx->L = 16;
667:   } else {
669:     ctx->L = bs;
670:   }
671:   oM = ctx->M;
672:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
673:     ctx->M = ctx->N/4;
674:   } else {
677:     ctx->M = ms;
678:   }
679:   onpart = ctx->npart;
680:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
681:     ctx->npart = 1;
682:   } else {
683:     MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
685:     ctx->npart = npart;
686:   }
687:   oLmax = ctx->L_max;
688:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
689:     ctx->L_max = 64;
690:   } else {
692:     ctx->L_max = PetscMax(bsmax,ctx->L);
693:   }
694:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
695:     SlepcContourDataDestroy(&ctx->contour);
696:     PetscInfo(eps,"Resetting the contour data structure due to a change of parameters\n");
697:     eps->state = EPS_STATE_INITIAL;
698:   }
699:   ctx->isreal = realmats;
700:   if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) eps->state = EPS_STATE_INITIAL;
701:   return 0;
702: }

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

707:    Logically Collective on eps

709:    Input Parameters:
710: +  eps   - the eigenproblem solver context
711: .  ip    - number of integration points
712: .  bs    - block size
713: .  ms    - moment size
714: .  npart - number of partitions when splitting the communicator
715: .  bsmax - max block size
716: -  realmats - A and B are real

718:    Options Database Keys:
719: +  -eps_ciss_integration_points - Sets the number of integration points
720: .  -eps_ciss_blocksize - Sets the block size
721: .  -eps_ciss_moments - Sets the moment size
722: .  -eps_ciss_partitions - Sets the number of partitions
723: .  -eps_ciss_maxblocksize - Sets the maximum block size
724: -  -eps_ciss_realmats - A and B are real

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

731:    Level: advanced

733: .seealso: EPSCISSGetSizes()
734: @*/
735: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
736: {
744:   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
745:   return 0;
746: }

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

752:   if (ip) *ip = ctx->N;
753:   if (bs) *bs = ctx->L;
754:   if (ms) *ms = ctx->M;
755:   if (npart) *npart = ctx->npart;
756:   if (bsmax) *bsmax = ctx->L_max;
757:   if (realmats) *realmats = ctx->isreal;
758:   return 0;
759: }

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

764:    Not Collective

766:    Input Parameter:
767: .  eps - the eigenproblem solver context

769:    Output Parameters:
770: +  ip    - number of integration points
771: .  bs    - block size
772: .  ms    - moment size
773: .  npart - number of partitions when splitting the communicator
774: .  bsmax - max block size
775: -  realmats - A and B are real

777:    Level: advanced

779: .seealso: EPSCISSSetSizes()
780: @*/
781: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
782: {
784:   PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
785:   return 0;
786: }

788: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
789: {
790:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

792:   if (delta == PETSC_DEFAULT) {
793:     ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
794:   } else {
796:     ctx->delta = delta;
797:   }
798:   if (spur == PETSC_DEFAULT) {
799:     ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
800:   } else {
802:     ctx->spurious_threshold = spur;
803:   }
804:   return 0;
805: }

807: /*@
808:    EPSCISSSetThreshold - Sets the values of various threshold parameters in
809:    the CISS solver.

811:    Logically Collective on eps

813:    Input Parameters:
814: +  eps   - the eigenproblem solver context
815: .  delta - threshold for numerical rank
816: -  spur  - spurious threshold (to discard spurious eigenpairs)

818:    Options Database Keys:
819: +  -eps_ciss_delta - Sets the delta
820: -  -eps_ciss_spurious_threshold - Sets the spurious threshold

822:    Level: advanced

824: .seealso: EPSCISSGetThreshold()
825: @*/
826: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
827: {
831:   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
832:   return 0;
833: }

835: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
836: {
837:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

839:   if (delta) *delta = ctx->delta;
840:   if (spur)  *spur = ctx->spurious_threshold;
841:   return 0;
842: }

844: /*@
845:    EPSCISSGetThreshold - Gets the values of various threshold parameters
846:    in the CISS solver.

848:    Not Collective

850:    Input Parameter:
851: .  eps - the eigenproblem solver context

853:    Output Parameters:
854: +  delta - threshold for numerical rank
855: -  spur  - spurious threshold (to discard spurious eigenpairs)

857:    Level: advanced

859: .seealso: EPSCISSSetThreshold()
860: @*/
861: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
862: {
864:   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
865:   return 0;
866: }

868: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
869: {
870:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

872:   if (inner == PETSC_DEFAULT) {
873:     ctx->refine_inner = 0;
874:   } else {
876:     ctx->refine_inner = inner;
877:   }
878:   if (blsize == PETSC_DEFAULT) {
879:     ctx->refine_blocksize = 0;
880:   } else {
882:     ctx->refine_blocksize = blsize;
883:   }
884:   return 0;
885: }

887: /*@
888:    EPSCISSSetRefinement - Sets the values of various refinement parameters
889:    in the CISS solver.

891:    Logically Collective on eps

893:    Input Parameters:
894: +  eps    - the eigenproblem solver context
895: .  inner  - number of iterative refinement iterations (inner loop)
896: -  blsize - number of iterative refinement iterations (blocksize loop)

898:    Options Database Keys:
899: +  -eps_ciss_refine_inner - Sets number of inner iterations
900: -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations

902:    Level: advanced

904: .seealso: EPSCISSGetRefinement()
905: @*/
906: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
907: {
911:   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
912:   return 0;
913: }

915: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
916: {
917:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

919:   if (inner)  *inner = ctx->refine_inner;
920:   if (blsize) *blsize = ctx->refine_blocksize;
921:   return 0;
922: }

924: /*@
925:    EPSCISSGetRefinement - Gets the values of various refinement parameters
926:    in the CISS solver.

928:    Not Collective

930:    Input Parameter:
931: .  eps - the eigenproblem solver context

933:    Output Parameters:
934: +  inner  - number of iterative refinement iterations (inner loop)
935: -  blsize - number of iterative refinement iterations (blocksize loop)

937:    Level: advanced

939: .seealso: EPSCISSSetRefinement()
940: @*/
941: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
942: {
944:   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
945:   return 0;
946: }

948: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
949: {
950:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

952:   ctx->usest     = usest;
953:   ctx->usest_set = PETSC_TRUE;
954:   eps->state     = EPS_STATE_INITIAL;
955:   return 0;
956: }

958: /*@
959:    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
960:    use the ST object for the linear solves.

962:    Logically Collective on eps

964:    Input Parameters:
965: +  eps    - the eigenproblem solver context
966: -  usest  - boolean flag to use the ST object or not

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

971:    Level: advanced

973: .seealso: EPSCISSGetUseST()
974: @*/
975: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
976: {
979:   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
980:   return 0;
981: }

983: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
984: {
985:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

987:   *usest = ctx->usest;
988:   return 0;
989: }

991: /*@
992:    EPSCISSGetUseST - Gets the flag for using the ST object
993:    in the CISS solver.

995:    Not Collective

997:    Input Parameter:
998: .  eps - the eigenproblem solver context

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

1003:    Level: advanced

1005: .seealso: EPSCISSSetUseST()
1006: @*/
1007: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1008: {
1011:   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1012:   return 0;
1013: }

1015: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1016: {
1017:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1019:   if (ctx->quad != quad) {
1020:     ctx->quad  = quad;
1021:     eps->state = EPS_STATE_INITIAL;
1022:   }
1023:   return 0;
1024: }

1026: /*@
1027:    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.

1029:    Logically Collective on eps

1031:    Input Parameters:
1032: +  eps  - the eigenproblem solver context
1033: -  quad - the quadrature rule

1035:    Options Database Key:
1036: .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1037:                            'chebyshev')

1039:    Notes:
1040:    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).

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

1045:    Level: advanced

1047: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1048: @*/
1049: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1050: {
1053:   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1054:   return 0;
1055: }

1057: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1058: {
1059:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1061:   *quad = ctx->quad;
1062:   return 0;
1063: }

1065: /*@
1066:    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.

1068:    Not Collective

1070:    Input Parameter:
1071: .  eps - the eigenproblem solver context

1073:    Output Parameters:
1074: .  quad - quadrature rule

1076:    Level: advanced

1078: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1079: @*/
1080: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1081: {
1084:   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1085:   return 0;
1086: }

1088: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1089: {
1090:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1092:   if (ctx->extraction != extraction) {
1093:     ctx->extraction = extraction;
1094:     eps->state      = EPS_STATE_INITIAL;
1095:   }
1096:   return 0;
1097: }

1099: /*@
1100:    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.

1102:    Logically Collective on eps

1104:    Input Parameters:
1105: +  eps        - the eigenproblem solver context
1106: -  extraction - the extraction technique

1108:    Options Database Key:
1109: .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1110:                            'hankel')

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

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

1118:    Level: advanced

1120: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1121: @*/
1122: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1123: {
1126:   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1127:   return 0;
1128: }

1130: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1131: {
1132:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1134:   *extraction = ctx->extraction;
1135:   return 0;
1136: }

1138: /*@
1139:    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.

1141:    Not Collective

1143:    Input Parameter:
1144: .  eps - the eigenproblem solver context

1146:    Output Parameters:
1147: .  extraction - extraction technique

1149:    Level: advanced

1151: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1152: @*/
1153: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1154: {
1157:   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1158:   return 0;
1159: }

1161: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1162: {
1163:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
1164:   SlepcContourData contour;
1165:   PetscInt         i,nsplit;
1166:   PC               pc;
1167:   MPI_Comm         child;

1169:   if (!ctx->contour) {  /* initialize contour data structure first */
1170:     RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj);
1171:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour);
1172:   }
1173:   contour = ctx->contour;
1174:   if (!contour->ksp) {
1175:     PetscMalloc1(contour->npoints,&contour->ksp);
1176:     EPSGetST(eps,&eps->st);
1177:     STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL);
1178:     PetscSubcommGetChild(contour->subcomm,&child);
1179:     for (i=0;i<contour->npoints;i++) {
1180:       KSPCreate(child,&contour->ksp[i]);
1181:       PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1);
1182:       KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix);
1183:       KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_");
1184:       PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options);
1185:       KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE);
1186:       KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1187:       KSPGetPC(contour->ksp[i],&pc);
1188:       if (nsplit) {
1189:         KSPSetType(contour->ksp[i],KSPBCGS);
1190:         PCSetType(pc,PCBJACOBI);
1191:       } else {
1192:         KSPSetType(contour->ksp[i],KSPPREONLY);
1193:         PCSetType(pc,PCLU);
1194:       }
1195:     }
1196:   }
1197:   if (nsolve) *nsolve = contour->npoints;
1198:   if (ksp)    *ksp    = contour->ksp;
1199:   return 0;
1200: }

1202: /*@C
1203:    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1204:    the CISS solver.

1206:    Not Collective

1208:    Input Parameter:
1209: .  eps - the eigenproblem solver solver

1211:    Output Parameters:
1212: +  nsolve - number of solver objects
1213: -  ksp - array of linear solver object

1215:    Notes:
1216:    The number of KSP solvers is equal to the number of integration points divided by
1217:    the number of partitions. This value is halved in the case of real matrices with
1218:    a region centered at the real axis.

1220:    Level: advanced

1222: .seealso: EPSCISSSetSizes()
1223: @*/
1224: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1225: {
1227:   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1228:   return 0;
1229: }

1231: PetscErrorCode EPSReset_CISS(EPS eps)
1232: {
1233:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1235:   BVDestroy(&ctx->S);
1236:   BVDestroy(&ctx->V);
1237:   BVDestroy(&ctx->Y);
1238:   if (!ctx->usest) SlepcContourDataReset(ctx->contour);
1239:   BVDestroy(&ctx->pV);
1240:   return 0;
1241: }

1243: PetscErrorCode EPSSetFromOptions_CISS(EPS eps,PetscOptionItems *PetscOptionsObject)
1244: {
1245:   PetscReal         r3,r4;
1246:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
1247:   PetscBool         b1,b2,flg,flg2,flg3,flg4,flg5,flg6;
1248:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
1249:   EPSCISSQuadRule   quad;
1250:   EPSCISSExtraction extraction;

1252:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS CISS Options");

1254:     EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1);
1255:     PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg);
1256:     PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2);
1257:     PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3);
1258:     PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4);
1259:     PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5);
1260:     PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6);
1261:     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1);

1263:     EPSCISSGetThreshold(eps,&r3,&r4);
1264:     PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg);
1265:     PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2);
1266:     if (flg || flg2) EPSCISSSetThreshold(eps,r3,r4);

1268:     EPSCISSGetRefinement(eps,&i6,&i7);
1269:     PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg);
1270:     PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2);
1271:     if (flg || flg2) EPSCISSSetRefinement(eps,i6,i7);

1273:     EPSCISSGetUseST(eps,&b2);
1274:     PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg);
1275:     if (flg) EPSCISSSetUseST(eps,b2);

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

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

1283:   PetscOptionsHeadEnd();

1285:   if (!eps->rg) EPSGetRG(eps,&eps->rg);
1286:   RGSetFromOptions(eps->rg); /* this is necessary here to set useconj */
1287:   if (!ctx->contour || !ctx->contour->ksp) EPSCISSGetKSPs(eps,NULL,NULL);
1288:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1289:   for (i=0;i<ctx->contour->npoints;i++) KSPSetFromOptions(ctx->contour->ksp[i]);
1290:   PetscSubcommSetFromOptions(ctx->contour->subcomm);
1291:   return 0;
1292: }

1294: PetscErrorCode EPSDestroy_CISS(EPS eps)
1295: {
1296:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1298:   SlepcContourDataDestroy(&ctx->contour);
1299:   PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1300:   PetscFree(eps->data);
1301:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL);
1302:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL);
1303:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL);
1304:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL);
1305:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL);
1306:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL);
1307:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL);
1308:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL);
1309:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL);
1310:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL);
1311:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL);
1312:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL);
1313:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL);
1314:   return 0;
1315: }

1317: PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1318: {
1319:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1320:   PetscBool      isascii;
1321:   PetscViewer    sviewer;

1323:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1324:   if (isascii) {
1325:     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);
1326:     if (ctx->isreal) PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n");
1327:     PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1328:     PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize);
1329:     PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]);
1330:     PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]);
1331:     if (ctx->usest) PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n");
1332:     else {
1333:       if (!ctx->contour || !ctx->contour->ksp) EPSCISSGetKSPs(eps,NULL,NULL);
1334:       PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Something went wrong with EPSCISSGetKSPs()");
1335:       PetscViewerASCIIPushTab(viewer);
1336:       if (ctx->npart>1 && ctx->contour->subcomm) {
1337:         PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1338:         if (!ctx->contour->subcomm->color) KSPView(ctx->contour->ksp[0],sviewer);
1339:         PetscViewerFlush(sviewer);
1340:         PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1341:         PetscViewerFlush(viewer);
1342:         /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1343:         PetscViewerASCIIPopSynchronized(viewer);
1344:       } else KSPView(ctx->contour->ksp[0],viewer);
1345:       PetscViewerASCIIPopTab(viewer);
1346:     }
1347:   }
1348:   return 0;
1349: }

1351: PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1352: {
1353:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1354:   PetscBool      usest = ctx->usest;
1355:   KSP            ksp;
1356:   PC             pc;

1358:   if (!((PetscObject)eps->st)->type_name) {
1359:     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1360:     if (usest) STSetType(eps->st,STSINVERT);
1361:     else {
1362:       /* we are not going to use ST, so avoid factorizing the matrix */
1363:       STSetType(eps->st,STSHIFT);
1364:       if (eps->isgeneralized) {
1365:         STGetKSP(eps->st,&ksp);
1366:         KSPGetPC(ksp,&pc);
1367:         PCSetType(pc,PCNONE);
1368:       }
1369:     }
1370:   }
1371:   return 0;
1372: }

1374: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1375: {
1376:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1378:   PetscNew(&ctx);
1379:   eps->data = ctx;

1381:   eps->useds = PETSC_TRUE;
1382:   eps->categ = EPS_CATEGORY_CONTOUR;

1384:   eps->ops->solve          = EPSSolve_CISS;
1385:   eps->ops->setup          = EPSSetUp_CISS;
1386:   eps->ops->setupsort      = EPSSetUpSort_CISS;
1387:   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1388:   eps->ops->destroy        = EPSDestroy_CISS;
1389:   eps->ops->reset          = EPSReset_CISS;
1390:   eps->ops->view           = EPSView_CISS;
1391:   eps->ops->computevectors = EPSComputeVectors_CISS;
1392:   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;

1394:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS);
1395:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS);
1396:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS);
1397:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS);
1398:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS);
1399:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS);
1400:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS);
1401:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS);
1402:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS);
1403:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS);
1404:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS);
1405:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS);
1406:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS);

1408:   /* set default values of parameters */
1409:   ctx->N                  = 32;
1410:   ctx->L                  = 16;
1411:   ctx->M                  = ctx->N/4;
1412:   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
1413:   ctx->L_max              = 64;
1414:   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1415:   ctx->usest              = PETSC_TRUE;
1416:   ctx->usest_set          = PETSC_FALSE;
1417:   ctx->isreal             = PETSC_FALSE;
1418:   ctx->refine_inner       = 0;
1419:   ctx->refine_blocksize   = 0;
1420:   ctx->npart              = 1;
1421:   ctx->quad               = (EPSCISSQuadRule)0;
1422:   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
1423:   return 0;
1424: }