Actual source code: ciss.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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: static PetscErrorCode EPSCISSSolveSystem(EPS eps,Mat A,Mat B,BV V,PetscInt L_start,PetscInt L_end,PetscBool initksp)
 69: {
 70:   PetscErrorCode   ierr;
 71:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
 72:   SlepcContourData contour;
 73:   PetscInt         i,p_id;
 74:   Mat              Fz,kspMat,MV,BMV=NULL,MC;
 75:   KSP              ksp;
 76:   const char       *prefix;

 79:   if (!ctx->contour || !ctx->contour->ksp) { EPSCISSGetKSPs(eps,NULL,NULL); }
 80:   contour = ctx->contour;
 81:   if (ctx->usest) {
 82:     MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Fz);
 83:   }
 84:   BVSetActiveColumns(V,L_start,L_end);
 85:   BVGetMat(V,&MV);
 86:   if (B) {
 87:     MatProductCreate(B,MV,NULL,&BMV);
 88:     MatProductSetType(BMV,MATPRODUCT_AB);
 89:     MatProductSetFromOptions(BMV);
 90:     MatProductSymbolic(BMV);
 91:   }
 92:   for (i=0;i<contour->npoints;i++) {
 93:     p_id = i*contour->subcomm->n + contour->subcomm->color;
 94:     if (!ctx->usest && initksp) {
 95:       MatDuplicate(A,MAT_COPY_VALUES,&kspMat);
 96:       if (B) {
 97:         MatAXPY(kspMat,-ctx->omega[p_id],B,UNKNOWN_NONZERO_PATTERN);
 98:       } else {
 99:         MatShift(kspMat,-ctx->omega[p_id]);
100:       }
101:       KSPSetOperators(contour->ksp[i],kspMat,kspMat);
102:       /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS) */
103:       KSPGetOptionsPrefix(contour->ksp[i],&prefix);
104:       MatSetOptionsPrefix(kspMat,prefix);
105:       MatDestroy(&kspMat);
106:     } else if (ctx->usest) {
107:       STSetShift(eps->st,ctx->omega[p_id]);
108:       STGetKSP(eps->st,&ksp);
109:     }
110:     BVSetActiveColumns(ctx->Y,i*ctx->L_max+L_start,i*ctx->L_max+L_end);
111:     BVGetMat(ctx->Y,&MC);
112:     if (B) {
113:       MatProductNumeric(BMV);
114:       if (ctx->usest) {
115:         KSPMatSolve(ksp,BMV,MC);
116:       } else {
117:         KSPMatSolve(contour->ksp[i],BMV,MC);
118:       }
119:     } else {
120:       if (ctx->usest) {
121:         KSPMatSolve(ksp,MV,MC);
122:       } else {
123:         KSPMatSolve(contour->ksp[i],MV,MC);
124:       }
125:     }
126:     if (ctx->usest && i<contour->npoints-1) { KSPReset(ksp); }
127:     BVRestoreMat(ctx->Y,&MC);
128:   }
129:   MatDestroy(&BMV);
130:   BVRestoreMat(V,&MV);
131:   if (ctx->usest) { MatDestroy(&Fz); }
132:   return(0);
133: }

135: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
136: {
138:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
139:   PetscInt       i;
140:   PetscScalar    center;
141:   PetscReal      radius,a,b,c,d,rgscale;
142: #if defined(PETSC_USE_COMPLEX)
143:   PetscReal      start_ang,end_ang,vscale,theta;
144: #endif
145:   PetscBool      isring,isellipse,isinterval;

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

203: PetscErrorCode EPSSetUp_CISS(EPS eps)
204: {
205:   PetscErrorCode   ierr;
206:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
207:   SlepcContourData contour;
208:   PetscBool        istrivial,isring,isellipse,isinterval,flg;
209:   PetscReal        c,d;
210:   PetscRandom      rand;
211:   PetscObjectId    id;
212:   PetscObjectState state;
213:   Mat              A[2];
214:   Vec              v0;

217:   if (eps->ncv==PETSC_DEFAULT) {
218:     eps->ncv = ctx->L_max*ctx->M;
219:     if (eps->ncv>eps->n) {
220:       eps->ncv = eps->n;
221:       ctx->L_max = eps->ncv/ctx->M;
222:       if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
223:     }
224:   } else {
225:     EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
226:     ctx->L_max = eps->ncv/ctx->M;
227:     if (!ctx->L_max) {
228:       ctx->L_max = 1;
229:       eps->ncv = ctx->L_max*ctx->M;
230:     }
231:   }
232:   ctx->L = PetscMin(ctx->L,ctx->L_max);
233:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 5;
234:   if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
235:   if (!eps->which) eps->which = EPS_ALL;
236:   if (eps->which!=EPS_ALL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
237:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);

239:   /* check region */
240:   RGIsTrivial(eps->rg,&istrivial);
241:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
242:   RGGetComplement(eps->rg,&flg);
243:   if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
244:   PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
245:   PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
246:   PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
247:   if (!isellipse && !isring && !isinterval) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");

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

258: #if !defined(PETSC_USE_COMPLEX)
259:   if (isring) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
260: #endif
261:   if (isinterval) {
262:     RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
263: #if !defined(PETSC_USE_COMPLEX)
264:     if (c!=d || c!=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"In real scalars, endpoints of the imaginary axis must be both zero");
265: #endif
266:     if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
267:   }
268:   if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;

270:   /* create contour data structure */
271:   if (!ctx->contour) {
272:     RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj);
273:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour);
274:   }

276:   EPSAllocateSolution(eps,0);
277:   BVGetRandomContext(eps->V,&rand);  /* make sure the random context is available when duplicating */
278:   if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
279:   PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
280:   PetscLogObjectMemory((PetscObject)eps,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));

282:   /* allocate basis vectors */
283:   BVDestroy(&ctx->S);
284:   BVDuplicateResize(eps->V,ctx->L_max*ctx->M,&ctx->S);
285:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->S);
286:   BVDestroy(&ctx->V);
287:   BVDuplicateResize(eps->V,ctx->L_max,&ctx->V);
288:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->V);

290:   STGetMatrix(eps->st,0,&A[0]);
291:   MatIsShell(A[0],&flg);
292:   if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
293:   if (eps->isgeneralized) { STGetMatrix(eps->st,1,&A[1]); }

295:   if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
296:   if (ctx->usest && ctx->npart>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");

298:   contour = ctx->contour;
299:   SlepcContourRedundantMat(contour,eps->isgeneralized?2:1,A);
300:   if (contour->pA) {
301:     BVGetColumn(ctx->V,0,&v0);
302:     SlepcContourScatterCreate(contour,v0);
303:     BVRestoreColumn(ctx->V,0,&v0);
304:     BVDestroy(&ctx->pV);
305:     BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV);
306:     BVSetSizesFromVec(ctx->pV,contour->xsub,eps->n);
307:     BVSetFromOptions(ctx->pV);
308:     BVResize(ctx->pV,ctx->L_max,PETSC_FALSE);
309:     PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->pV);
310:   }

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

315:   BVDestroy(&ctx->Y);
316:   if (contour->pA) {
317:     BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y);
318:     BVSetSizesFromVec(ctx->Y,contour->xsub,eps->n);
319:     BVSetFromOptions(ctx->Y);
320:     BVResize(ctx->Y,contour->npoints*ctx->L_max,PETSC_FALSE);
321:   } else {
322:     BVDuplicateResize(eps->V,contour->npoints*ctx->L_max,&ctx->Y);
323:   }
324:   PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->Y);

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

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

352: PetscErrorCode EPSSetUpSort_CISS(EPS eps)
353: {
355:   SlepcSC        sc;

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

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

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

393:   w[0] = eps->work[0];
394: #if defined(PETSC_USE_COMPLEX)
395:   w[1] = NULL;
396: #else
397:   w[1] = eps->work[2];
398: #endif
399:   w[2] = eps->work[1];
400:   VecGetLocalSize(w[0],&nlocal);
401:   DSGetLeadingDimension(eps->ds,&ld);
402:   STGetNumMatrices(eps->st,&nmat);
403:   STGetMatrix(eps->st,0,&A);
404:   if (nmat>1) { STGetMatrix(eps->st,1,&B); }
405:   else B = NULL;
406:   RGComputeQuadrature(eps->rg,ctx->quad==EPS_CISS_QUADRULE_CHEBYSHEV?RG_QUADRULE_CHEBYSHEV:RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight);
407:   BVSetActiveColumns(ctx->V,0,ctx->L);
408:   BVSetRandomSign(ctx->V);
409:   BVGetRandomContext(ctx->V,&rand);

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

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

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

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

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

563:       DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
564:       DSGetMat(eps->ds,DS_MAT_X,&X);
565:       BVMultInPlace(ctx->S,X,0,eps->nconv);
566:       if (eps->ishermitian) {
567:         BVMultInPlace(eps->V,X,0,eps->nconv);
568:       }
569:       MatDestroy(&X);
570:       max_error = 0.0;
571:       for (i=0;i<eps->nconv;i++) {
572:         BVGetColumn(ctx->S,i,&si);
573: #if !defined(PETSC_USE_COMPLEX)
574:         if (eps->eigi[i]!=0.0) { BVGetColumn(ctx->S,i+1,&si1); }
575: #endif
576:         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:           VecNorm(si,NORM_2,&norm);
579: #if !defined(PETSC_USE_COMPLEX)
580:           if (eps->eigi[i]!=0.0) {
581:             VecNorm(si1,NORM_2,&normi);
582:             norm = SlepcAbsEigenvalue(norm,normi);
583:           }
584: #endif
585:           error /= norm;
586:         }
587:         (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx);
588:         BVRestoreColumn(ctx->S,i,&si);
589: #if !defined(PETSC_USE_COMPLEX)
590:         if (eps->eigi[i]!=0.0) {
591:           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:         MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
604:         MatSetRandom(M,rand);
605:         BVSetActiveColumns(ctx->S,0,nv);
606:         BVMultInPlace(ctx->S,M,0,ctx->L);
607:         MatDestroy(&M);
608:         BVSetActiveColumns(ctx->S,0,ctx->L);
609:         BVSetActiveColumns(ctx->V,0,ctx->L);
610:         BVCopy(ctx->S,ctx->V);
611:         if (contour->pA) {
612:           BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
613:           EPSCISSSolveSystem(eps,contour->pA[0],contour->pA[1],ctx->pV,0,ctx->L,PETSC_FALSE);
614:         } else {
615:           EPSCISSSolveSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
616:         }
617:       }
618:     }
619:   }
620:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
621:     PetscFree(H1);
622:   }
623:   PetscFree2(Mu,H0);
624:   return(0);
625: }

627: PetscErrorCode EPSComputeVectors_CISS(EPS eps)
628: {
630:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
631:   PetscInt       n;
632:   Mat            Z,B=NULL;

635:   if (eps->ishermitian) {
636:     if (eps->isgeneralized && !eps->ispositive) {
637:       EPSComputeVectors_Indefinite(eps);
638:     } else {
639:       EPSComputeVectors_Hermitian(eps);
640:     }
641:     if (eps->isgeneralized && eps->ispositive && ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
642:       /* normalize to have unit B-norm */
643:       STGetMatrix(eps->st,1,&B);
644:       BVSetMatrix(eps->V,B,PETSC_FALSE);
645:       BVNormalize(eps->V,NULL);
646:       BVSetMatrix(eps->V,NULL,PETSC_FALSE);
647:     }
648:     return(0);
649:   }
650:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
651:   BVSetActiveColumns(eps->V,0,n);

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

656:   /* V = V * Z */
657:   DSGetMat(eps->ds,DS_MAT_X,&Z);
658:   BVMultInPlace(eps->V,Z,0,n);
659:   MatDestroy(&Z);
660:   BVSetActiveColumns(eps->V,0,eps->nconv);

662:   /* normalize */
663:   if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
664:     BVNormalize(eps->V,NULL);
665:   }
666:   return(0);
667: }

669: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
670: {
672:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
673:   PetscInt       oN,oL,oM,oLmax,onpart;

676:   oN = ctx->N;
677:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
678:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
679:   } else {
680:     if (ip<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
681:     if (ip%2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
682:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
683:   }
684:   oL = ctx->L;
685:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
686:     ctx->L = 16;
687:   } else {
688:     if (bs<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
689:     ctx->L = bs;
690:   }
691:   oM = ctx->M;
692:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
693:     ctx->M = ctx->N/4;
694:   } else {
695:     if (ms<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
696:     if (ms>ctx->N) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
697:     ctx->M = ms;
698:   }
699:   onpart = ctx->npart;
700:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
701:     ctx->npart = 1;
702:   } else {
703:     if (npart<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
704:     ctx->npart = npart;
705:   }
706:   oLmax = ctx->L_max;
707:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
708:     ctx->L_max = 64;
709:   } else {
710:     if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
711:     ctx->L_max = PetscMax(bsmax,ctx->L);
712:   }
713:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
714:     SlepcContourDataDestroy(&ctx->contour);
715:     PetscInfo(eps,"Resetting the contour data structure due to a change of parameters\n");
716:     eps->state = EPS_STATE_INITIAL;
717:   }
718:   ctx->isreal = realmats;
719:   if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) eps->state = EPS_STATE_INITIAL;
720:   return(0);
721: }

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

726:    Logically Collective on eps

728:    Input Parameters:
729: +  eps   - the eigenproblem solver context
730: .  ip    - number of integration points
731: .  bs    - block size
732: .  ms    - moment size
733: .  npart - number of partitions when splitting the communicator
734: .  bsmax - max block size
735: -  realmats - A and B are real

737:    Options Database Keys:
738: +  -eps_ciss_integration_points - Sets the number of integration points
739: .  -eps_ciss_blocksize - Sets the block size
740: .  -eps_ciss_moments - Sets the moment size
741: .  -eps_ciss_partitions - Sets the number of partitions
742: .  -eps_ciss_maxblocksize - Sets the maximum block size
743: -  -eps_ciss_realmats - A and B are real

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

750:    Level: advanced

752: .seealso: EPSCISSGetSizes()
753: @*/
754: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
755: {

766:   PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
767:   return(0);
768: }

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

775:   if (ip) *ip = ctx->N;
776:   if (bs) *bs = ctx->L;
777:   if (ms) *ms = ctx->M;
778:   if (npart) *npart = ctx->npart;
779:   if (bsmax) *bsmax = ctx->L_max;
780:   if (realmats) *realmats = ctx->isreal;
781:   return(0);
782: }

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

787:    Not Collective

789:    Input Parameter:
790: .  eps - the eigenproblem solver context

792:    Output Parameters:
793: +  ip    - number of integration points
794: .  bs    - block size
795: .  ms    - moment size
796: .  npart - number of partitions when splitting the communicator
797: .  bsmax - max block size
798: -  realmats - A and B are real

800:    Level: advanced

802: .seealso: EPSCISSSetSizes()
803: @*/
804: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
805: {

810:   PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
811:   return(0);
812: }

814: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
815: {
816:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

819:   if (delta == PETSC_DEFAULT) {
820:     ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
821:   } else {
822:     if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
823:     ctx->delta = delta;
824:   }
825:   if (spur == PETSC_DEFAULT) {
826:     ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
827:   } else {
828:     if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
829:     ctx->spurious_threshold = spur;
830:   }
831:   return(0);
832: }

834: /*@
835:    EPSCISSSetThreshold - Sets the values of various threshold parameters in
836:    the CISS solver.

838:    Logically Collective on eps

840:    Input Parameters:
841: +  eps   - the eigenproblem solver context
842: .  delta - threshold for numerical rank
843: -  spur  - spurious threshold (to discard spurious eigenpairs)

845:    Options Database Keys:
846: +  -eps_ciss_delta - Sets the delta
847: -  -eps_ciss_spurious_threshold - Sets the spurious threshold

849:    Level: advanced

851: .seealso: EPSCISSGetThreshold()
852: @*/
853: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
854: {

861:   PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
862:   return(0);
863: }

865: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
866: {
867:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

870:   if (delta) *delta = ctx->delta;
871:   if (spur)  *spur = ctx->spurious_threshold;
872:   return(0);
873: }

875: /*@
876:    EPSCISSGetThreshold - Gets the values of various threshold parameters
877:    in the CISS solver.

879:    Not Collective

881:    Input Parameter:
882: .  eps - the eigenproblem solver context

884:    Output Parameters:
885: +  delta - threshold for numerical rank
886: -  spur  - spurious threshold (to discard spurious eigenpairs)

888:    Level: advanced

890: .seealso: EPSCISSSetThreshold()
891: @*/
892: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
893: {

898:   PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
899:   return(0);
900: }

902: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
903: {
904:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

907:   if (inner == PETSC_DEFAULT) {
908:     ctx->refine_inner = 0;
909:   } else {
910:     if (inner<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
911:     ctx->refine_inner = inner;
912:   }
913:   if (blsize == PETSC_DEFAULT) {
914:     ctx->refine_blocksize = 0;
915:   } else {
916:     if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
917:     ctx->refine_blocksize = blsize;
918:   }
919:   return(0);
920: }

922: /*@
923:    EPSCISSSetRefinement - Sets the values of various refinement parameters
924:    in the CISS solver.

926:    Logically Collective on eps

928:    Input Parameters:
929: +  eps    - the eigenproblem solver context
930: .  inner  - number of iterative refinement iterations (inner loop)
931: -  blsize - number of iterative refinement iterations (blocksize loop)

933:    Options Database Keys:
934: +  -eps_ciss_refine_inner - Sets number of inner iterations
935: -  -eps_ciss_refine_blocksize - Sets number of blocksize iterations

937:    Level: advanced

939: .seealso: EPSCISSGetRefinement()
940: @*/
941: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
942: {

949:   PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
950:   return(0);
951: }

953: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
954: {
955:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

958:   if (inner)  *inner = ctx->refine_inner;
959:   if (blsize) *blsize = ctx->refine_blocksize;
960:   return(0);
961: }

963: /*@
964:    EPSCISSGetRefinement - Gets the values of various refinement parameters
965:    in the CISS solver.

967:    Not Collective

969:    Input Parameter:
970: .  eps - the eigenproblem solver context

972:    Output Parameters:
973: +  inner  - number of iterative refinement iterations (inner loop)
974: -  blsize - number of iterative refinement iterations (blocksize loop)

976:    Level: advanced

978: .seealso: EPSCISSSetRefinement()
979: @*/
980: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
981: {

986:   PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
987:   return(0);
988: }

990: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
991: {
992:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

995:   ctx->usest     = usest;
996:   ctx->usest_set = PETSC_TRUE;
997:   eps->state     = EPS_STATE_INITIAL;
998:   return(0);
999: }

1001: /*@
1002:    EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
1003:    use the ST object for the linear solves.

1005:    Logically Collective on eps

1007:    Input Parameters:
1008: +  eps    - the eigenproblem solver context
1009: -  usest  - boolean flag to use the ST object or not

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

1014:    Level: advanced

1016: .seealso: EPSCISSGetUseST()
1017: @*/
1018: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
1019: {

1025:   PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1026:   return(0);
1027: }

1029: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1030: {
1031:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1034:   *usest = ctx->usest;
1035:   return(0);
1036: }

1038: /*@
1039:    EPSCISSGetUseST - Gets the flag for using the ST object
1040:    in the CISS solver.

1042:    Not Collective

1044:    Input Parameter:
1045: .  eps - the eigenproblem solver context

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

1050:    Level: advanced

1052: .seealso: EPSCISSSetUseST()
1053: @*/
1054: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1055: {

1061:   PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1062:   return(0);
1063: }

1065: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1066: {
1067:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1070:   if (ctx->quad != quad) {
1071:     ctx->quad  = quad;
1072:     eps->state = EPS_STATE_INITIAL;
1073:   }
1074:   return(0);
1075: }

1077: /*@
1078:    EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.

1080:    Logically Collective on eps

1082:    Input Parameters:
1083: +  eps  - the eigenproblem solver context
1084: -  quad - the quadrature rule

1086:    Options Database Key:
1087: .  -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1088:                            'chebyshev')

1090:    Notes:
1091:    By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).

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

1096:    Level: advanced

1098: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1099: @*/
1100: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1101: {

1107:   PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1108:   return(0);
1109: }

1111: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1112: {
1113:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1116:   *quad = ctx->quad;
1117:   return(0);
1118: }

1120: /*@
1121:    EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.

1123:    Not Collective

1125:    Input Parameter:
1126: .  eps - the eigenproblem solver context

1128:    Output Parameters:
1129: .  quad - quadrature rule

1131:    Level: advanced

1133: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1134: @*/
1135: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1136: {

1142:   PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1143:   return(0);
1144: }

1146: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1147: {
1148:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1151:   if (ctx->extraction != extraction) {
1152:     ctx->extraction = extraction;
1153:     eps->state      = EPS_STATE_INITIAL;
1154:   }
1155:   return(0);
1156: }

1158: /*@
1159:    EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.

1161:    Logically Collective on eps

1163:    Input Parameters:
1164: +  eps        - the eigenproblem solver context
1165: -  extraction - the extraction technique

1167:    Options Database Key:
1168: .  -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1169:                            'hankel')

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

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

1177:    Level: advanced

1179: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1180: @*/
1181: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1182: {

1188:   PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1189:   return(0);
1190: }

1192: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1193: {
1194:   EPS_CISS *ctx = (EPS_CISS*)eps->data;

1197:   *extraction = ctx->extraction;
1198:   return(0);
1199: }

1201: /*@
1202:    EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.

1204:    Not Collective

1206:    Input Parameter:
1207: .  eps - the eigenproblem solver context

1209:    Output Parameters:
1210: .  extraction - extraction technique

1212:    Level: advanced

1214: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1215: @*/
1216: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1217: {

1223:   PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1224:   return(0);
1225: }

1227: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1228: {
1229:   PetscErrorCode   ierr;
1230:   EPS_CISS         *ctx = (EPS_CISS*)eps->data;
1231:   SlepcContourData contour;
1232:   PetscInt         i;
1233:   PC               pc;

1236:   if (!ctx->contour) {  /* initialize contour data structure first */
1237:     RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj);
1238:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour);
1239:   }
1240:   contour = ctx->contour;
1241:   if (!contour->ksp) {
1242:     PetscMalloc1(contour->npoints,&contour->ksp);
1243:     for (i=0;i<contour->npoints;i++) {
1244:       KSPCreate(PetscSubcommChild(contour->subcomm),&contour->ksp[i]);
1245:       PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1);
1246:       KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix);
1247:       KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_");
1248:       PetscLogObjectParent((PetscObject)eps,(PetscObject)contour->ksp[i]);
1249:       PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options);
1250:       KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE);
1251:       KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1252:       KSPGetPC(contour->ksp[i],&pc);
1253:       KSPSetType(contour->ksp[i],KSPPREONLY);
1254:       PCSetType(pc,PCLU);
1255:     }
1256:   }
1257:   if (nsolve) *nsolve = contour->npoints;
1258:   if (ksp)    *ksp    = contour->ksp;
1259:   return(0);
1260: }

1262: /*@C
1263:    EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1264:    the CISS solver.

1266:    Not Collective

1268:    Input Parameter:
1269: .  eps - the eigenproblem solver solver

1271:    Output Parameters:
1272: +  nsolve - number of solver objects
1273: -  ksp - array of linear solver object

1275:    Notes:
1276:    The number of KSP solvers is equal to the number of integration points divided by
1277:    the number of partitions. This value is halved in the case of real matrices with
1278:    a region centered at the real axis.

1280:    Level: advanced

1282: .seealso: EPSCISSSetSizes()
1283: @*/
1284: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1285: {

1290:   PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1291:   return(0);
1292: }

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

1300:   BVDestroy(&ctx->S);
1301:   BVDestroy(&ctx->V);
1302:   BVDestroy(&ctx->Y);
1303:   if (!ctx->usest) {
1304:     SlepcContourDataReset(ctx->contour);
1305:   }
1306:   BVDestroy(&ctx->pV);
1307:   return(0);
1308: }

1310: PetscErrorCode EPSSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,EPS eps)
1311: {
1312:   PetscErrorCode    ierr;
1313:   PetscReal         r3,r4;
1314:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
1315:   PetscBool         b1,b2,flg,flg2,flg3,flg4,flg5,flg6;
1316:   EPS_CISS          *ctx = (EPS_CISS*)eps->data;
1317:   EPSCISSQuadRule   quad;
1318:   EPSCISSExtraction extraction;

1321:   PetscOptionsHead(PetscOptionsObject,"EPS CISS Options");

1323:     EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1);
1324:     PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg);
1325:     PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2);
1326:     PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3);
1327:     PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4);
1328:     PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5);
1329:     PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6);
1330:     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) { EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1); }

1332:     EPSCISSGetThreshold(eps,&r3,&r4);
1333:     PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg);
1334:     PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2);
1335:     if (flg || flg2) { EPSCISSSetThreshold(eps,r3,r4); }

1337:     EPSCISSGetRefinement(eps,&i6,&i7);
1338:     PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg);
1339:     PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2);
1340:     if (flg || flg2) { EPSCISSSetRefinement(eps,i6,i7); }

1342:     EPSCISSGetUseST(eps,&b2);
1343:     PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg);
1344:     if (flg) { EPSCISSSetUseST(eps,b2); }

1346:     PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg);
1347:     if (flg) { EPSCISSSetQuadRule(eps,quad); }

1349:     PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
1350:     if (flg) { EPSCISSSetExtraction(eps,extraction); }

1352:   PetscOptionsTail();

1354:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
1355:   RGSetFromOptions(eps->rg); /* this is necessary here to set useconj */
1356:   if (!ctx->contour || !ctx->contour->ksp) { EPSCISSGetKSPs(eps,NULL,NULL); }
1357:   for (i=0;i<ctx->contour->npoints;i++) {
1358:     KSPSetFromOptions(ctx->contour->ksp[i]);
1359:   }
1360:   PetscSubcommSetFromOptions(ctx->contour->subcomm);
1361:   return(0);
1362: }

1364: PetscErrorCode EPSDestroy_CISS(EPS eps)
1365: {
1367:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1370:   SlepcContourDataDestroy(&ctx->contour);
1371:   PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1372:   PetscFree(eps->data);
1373:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL);
1374:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL);
1375:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL);
1376:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL);
1377:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL);
1378:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL);
1379:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL);
1380:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL);
1381:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL);
1382:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL);
1383:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL);
1384:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL);
1385:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL);
1386:   return(0);
1387: }

1389: PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1390: {
1392:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1393:   PetscBool      isascii;
1394:   PetscViewer    sviewer;

1397:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1398:   if (isascii) {
1399:     PetscViewerASCIIPrintf(viewer,"  sizes { integration points: %D, block size: %D, moment size: %D, partitions: %D, maximum block size: %D }\n",ctx->N,ctx->L,ctx->M,ctx->npart,ctx->L_max);
1400:     if (ctx->isreal) {
1401:       PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n");
1402:     }
1403:     PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1404:     PetscViewerASCIIPrintf(viewer,"  iterative refinement { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1405:     PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",EPSCISSExtractions[ctx->extraction]);
1406:     PetscViewerASCIIPrintf(viewer,"  quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]);
1407:     if (ctx->usest) {
1408:       PetscViewerASCIIPrintf(viewer,"  using ST for linear solves\n");
1409:     } else {
1410:       if (!ctx->contour || !ctx->contour->ksp) { EPSCISSGetKSPs(eps,NULL,NULL); }
1411:       PetscViewerASCIIPushTab(viewer);
1412:       if (ctx->npart>1 && ctx->contour->subcomm) {
1413:         PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1414:         if (!ctx->contour->subcomm->color) {
1415:           KSPView(ctx->contour->ksp[0],sviewer);
1416:         }
1417:         PetscViewerFlush(sviewer);
1418:         PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1419:         PetscViewerFlush(viewer);
1420:         /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1421:         PetscViewerASCIIPopSynchronized(viewer);
1422:       } else {
1423:         KSPView(ctx->contour->ksp[0],viewer);
1424:       }
1425:       PetscViewerASCIIPopTab(viewer);
1426:     }
1427:   }
1428:   return(0);
1429: }

1431: PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1432: {
1434:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;
1435:   PetscBool      usest = ctx->usest;
1436:   KSP            ksp;
1437:   PC             pc;

1440:   if (!((PetscObject)eps->st)->type_name) {
1441:     if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1442:     if (usest) {
1443:       STSetType(eps->st,STSINVERT);
1444:     } else {
1445:       /* we are not going to use ST, so avoid factorizing the matrix */
1446:       STSetType(eps->st,STSHIFT);
1447:       if (eps->isgeneralized) {
1448:         STGetKSP(eps->st,&ksp);
1449:         KSPGetPC(ksp,&pc);
1450:         PCSetType(pc,PCNONE);
1451:       }
1452:     }
1453:   }
1454:   return(0);
1455: }

1457: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1458: {
1460:   EPS_CISS       *ctx = (EPS_CISS*)eps->data;

1463:   PetscNewLog(eps,&ctx);
1464:   eps->data = ctx;

1466:   eps->useds = PETSC_TRUE;
1467:   eps->categ = EPS_CATEGORY_CONTOUR;

1469:   eps->ops->solve          = EPSSolve_CISS;
1470:   eps->ops->setup          = EPSSetUp_CISS;
1471:   eps->ops->setupsort      = EPSSetUpSort_CISS;
1472:   eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1473:   eps->ops->destroy        = EPSDestroy_CISS;
1474:   eps->ops->reset          = EPSReset_CISS;
1475:   eps->ops->view           = EPSView_CISS;
1476:   eps->ops->computevectors = EPSComputeVectors_CISS;
1477:   eps->ops->setdefaultst   = EPSSetDefaultST_CISS;

1479:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS);
1480:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS);
1481:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS);
1482:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS);
1483:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS);
1484:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS);
1485:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS);
1486:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS);
1487:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS);
1488:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS);
1489:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS);
1490:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS);
1491:   PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS);

1493:   /* set default values of parameters */
1494:   ctx->N                  = 32;
1495:   ctx->L                  = 16;
1496:   ctx->M                  = ctx->N/4;
1497:   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
1498:   ctx->L_max              = 64;
1499:   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1500:   ctx->usest              = PETSC_TRUE;
1501:   ctx->usest_set          = PETSC_FALSE;
1502:   ctx->isreal             = PETSC_FALSE;
1503:   ctx->refine_inner       = 0;
1504:   ctx->refine_blocksize   = 0;
1505:   ctx->npart              = 1;
1506:   ctx->quad               = (EPSCISSQuadRule)0;
1507:   ctx->extraction         = EPS_CISS_EXTRACTION_RITZ;
1508:   return(0);
1509: }