Actual source code: nciss.c

slepc-3.18.2 2023-01-26
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] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, K. Kimura, "A
 26:            numerical method for nonlinear eigenvalue problems using contour
 27:            integrals", JSIAM Lett. 1:52-55, 2009.

 29:        [2] S. Yokota and T. Sakurai, "A projection method for nonlinear
 30:            eigenvalue problems using contour integrals", JSIAM Lett.
 31:            5:41-44, 2013.
 32: */

 34: #include <slepc/private/nepimpl.h>
 35: #include <slepc/private/slepccontour.h>

 37: typedef struct _n_nep_ciss_project *NEP_CISS_PROJECT;
 38: typedef struct {
 39:   /* parameters */
 40:   PetscInt          N;             /* number of integration points (32) */
 41:   PetscInt          L;             /* block size (16) */
 42:   PetscInt          M;             /* moment degree (N/4 = 4) */
 43:   PetscReal         delta;         /* threshold of singular value (1e-12) */
 44:   PetscInt          L_max;         /* maximum number of columns of the source matrix V */
 45:   PetscReal         spurious_threshold; /* discard spurious eigenpairs */
 46:   PetscBool         isreal;        /* T(z) is real for real z */
 47:   PetscInt          npart;         /* number of partitions */
 48:   PetscInt          refine_inner;
 49:   PetscInt          refine_blocksize;
 50:   NEPCISSExtraction extraction;
 51:   /* private data */
 52:   SlepcContourData  contour;
 53:   PetscReal         *sigma;        /* threshold for numerical rank */
 54:   PetscScalar       *weight;
 55:   PetscScalar       *omega;
 56:   PetscScalar       *pp;
 57:   BV                V;
 58:   BV                S;
 59:   BV                Y;
 60:   PetscBool         useconj;
 61:   Mat               J;             /* auxiliary matrix when using subcomm */
 62:   BV                pV;
 63:   NEP_CISS_PROJECT  dsctxf;
 64:   PetscObjectId     rgid;
 65:   PetscObjectState  rgstate;
 66: } NEP_CISS;

 68: struct _n_nep_ciss_project {
 69:   NEP  nep;
 70:   BV   Q;
 71: };

 73: static PetscErrorCode NEPContourDSComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx)
 74: {
 75:   NEP_CISS_PROJECT proj = (NEP_CISS_PROJECT)ctx;
 76:   Mat              M,fun;

 78:   if (!deriv) {
 79:     NEPComputeFunction(proj->nep,lambda,proj->nep->function,proj->nep->function);
 80:     fun = proj->nep->function;
 81:   } else {
 82:     NEPComputeJacobian(proj->nep,lambda,proj->nep->jacobian);
 83:     fun = proj->nep->jacobian;
 84:   }
 85:   DSGetMat(ds,mat,&M);
 86:   BVMatProject(proj->Q,fun,proj->Q,M);
 87:   DSRestoreMat(ds,mat,&M);
 88:   return 0;
 89: }

 91: static PetscErrorCode NEPComputeFunctionSubcomm(NEP nep,PetscScalar lambda,Mat T,Mat P,PetscBool deriv)
 92: {
 93:   PetscInt       i;
 94:   PetscScalar    alpha;
 95:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

 97:   PetscAssert(nep->fui!=NEP_USER_INTERFACE_CALLBACK,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Should not arrive here with callbacks");
 98:   MatZeroEntries(T);
 99:   if (!deriv && T != P) MatZeroEntries(P);
100:   for (i=0;i<nep->nt;i++) {
101:     if (!deriv) FNEvaluateFunction(nep->f[i],lambda,&alpha);
102:     else FNEvaluateDerivative(nep->f[i],lambda,&alpha);
103:     MatAXPY(T,alpha,ctx->contour->pA[i],nep->mstr);
104:     if (!deriv && T != P) MatAXPY(P,alpha,ctx->contour->pP[i],nep->mstrp);
105:   }
106:   return 0;
107: }

109: /*
110:   Set up KSP solvers for every integration point
111: */
112: static PetscErrorCode NEPCISSSetUp(NEP nep,Mat T,Mat P)
113: {
114:   NEP_CISS         *ctx = (NEP_CISS*)nep->data;
115:   SlepcContourData contour;
116:   PetscInt         i,p_id;
117:   Mat              Amat,Pmat;

119:   if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
120:   contour = ctx->contour;
121:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
122:   for (i=0;i<contour->npoints;i++) {
123:     p_id = i*contour->subcomm->n + contour->subcomm->color;
124:     MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&Amat);
125:     if (T != P) MatDuplicate(P,MAT_DO_NOT_COPY_VALUES,&Pmat); else Pmat = Amat;
126:     if (contour->subcomm->n == 1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) NEPComputeFunction(nep,ctx->omega[p_id],Amat,Pmat);
127:     else NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],Amat,Pmat,PETSC_FALSE);
128:     NEP_KSPSetOperators(contour->ksp[i],Amat,Pmat);
129:     MatDestroy(&Amat);
130:     if (T != P) MatDestroy(&Pmat);
131:   }
132:   return 0;
133: }

135: /*
136:   Y_i = F(z_i)^{-1}Fp(z_i)V for every integration point, Y=[Y_i] is in the context
137: */
138: static PetscErrorCode NEPCISSSolve(NEP nep,Mat dT,BV V,PetscInt L_start,PetscInt L_end)
139: {
140:   NEP_CISS         *ctx = (NEP_CISS*)nep->data;
141:   SlepcContourData contour = ctx->contour;
142:   PetscInt         i,p_id;
143:   Mat              MV,BMV=NULL,MC;

145:   BVSetActiveColumns(V,L_start,L_end);
146:   BVGetMat(V,&MV);
147:   for (i=0;i<contour->npoints;i++) {
148:     p_id = i*contour->subcomm->n + contour->subcomm->color;
149:     if (contour->subcomm->n==1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) NEPComputeJacobian(nep,ctx->omega[p_id],dT);
150:     else NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],dT,NULL,PETSC_TRUE);
151:     BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end);
152:     BVGetMat(ctx->Y,&MC);
153:     if (!i) {
154:       MatProductCreate(dT,MV,NULL,&BMV);
155:       MatProductSetType(BMV,MATPRODUCT_AB);
156:       MatProductSetFromOptions(BMV);
157:       MatProductSymbolic(BMV);
158:     }
159:     MatProductNumeric(BMV);
160:     KSPMatSolve(contour->ksp[i],BMV,MC);
161:     BVRestoreMat(ctx->Y,&MC);
162:   }
163:   MatDestroy(&BMV);
164:   BVRestoreMat(V,&MV);
165:   return 0;
166: }

168: PetscErrorCode NEPSetUp_CISS(NEP nep)
169: {
170:   NEP_CISS         *ctx = (NEP_CISS*)nep->data;
171:   SlepcContourData contour;
172:   PetscInt         nwork;
173:   PetscBool        istrivial,isellipse,flg;
174:   NEP_CISS_PROJECT dsctxf;
175:   PetscObjectId    id;
176:   PetscObjectState state;
177:   Vec              v0;

179:   if (nep->ncv==PETSC_DEFAULT) nep->ncv = ctx->L_max*ctx->M;
180:   else {
181:     ctx->L_max = nep->ncv/ctx->M;
182:     if (!ctx->L_max) {
183:       ctx->L_max = 1;
184:       nep->ncv = ctx->L_max*ctx->M;
185:     }
186:   }
187:   ctx->L = PetscMin(ctx->L,ctx->L_max);
188:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = 5;
189:   if (nep->mpd==PETSC_DEFAULT) nep->mpd = nep->ncv;
190:   if (!nep->which) nep->which = NEP_ALL;
192:   NEPCheckUnsupported(nep,NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);

194:   /* check region */
195:   RGIsTrivial(nep->rg,&istrivial);
197:   RGGetComplement(nep->rg,&flg);
199:   PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);

202:   /* if the region has changed, then reset contour data */
203:   PetscObjectGetId((PetscObject)nep->rg,&id);
204:   PetscObjectStateGet((PetscObject)nep->rg,&state);
205:   if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
206:     SlepcContourDataDestroy(&ctx->contour);
207:     PetscInfo(nep,"Resetting the contour data structure due to a change of region\n");
208:     ctx->rgid = id; ctx->rgstate = state;
209:   }

211:   /* create contour data structure */
212:   if (!ctx->contour) {
213:     RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj);
214:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour);
215:   }

217:   NEPAllocateSolution(nep,0);
218:   if (ctx->weight) PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
219:   PetscMalloc4(ctx->N,&ctx->weight,ctx->N,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);

221:   /* allocate basis vectors */
222:   BVDestroy(&ctx->S);
223:   BVDuplicateResize(nep->V,ctx->L*ctx->M,&ctx->S);
224:   BVDestroy(&ctx->V);
225:   BVDuplicateResize(nep->V,ctx->L,&ctx->V);

227:   contour = ctx->contour;
228:   if (contour->subcomm && contour->subcomm->n != 1 && nep->fui==NEP_USER_INTERFACE_CALLBACK) {
229:     NEPComputeFunction(nep,0,nep->function,nep->function_pre);
230:     SlepcContourRedundantMat(contour,1,&nep->function,(nep->function!=nep->function_pre)?&nep->function_pre:NULL);
231:   } else SlepcContourRedundantMat(contour,nep->nt,nep->A,nep->P);
232:   if (contour->pA) {
233:     if (!ctx->J) MatDuplicate(contour->pA[0],MAT_DO_NOT_COPY_VALUES,&ctx->J);
234:     BVGetColumn(ctx->V,0,&v0);
235:     SlepcContourScatterCreate(contour,v0);
236:     BVRestoreColumn(ctx->V,0,&v0);
237:     BVDestroy(&ctx->pV);
238:     BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV);
239:     BVSetSizesFromVec(ctx->pV,contour->xsub,nep->n);
240:     BVSetFromOptions(ctx->pV);
241:     BVResize(ctx->pV,ctx->L,PETSC_FALSE);
242:   }

244:   BVDestroy(&ctx->Y);
245:   if (contour->pA) {
246:     BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y);
247:     BVSetSizesFromVec(ctx->Y,contour->xsub,nep->n);
248:     BVSetFromOptions(ctx->Y);
249:     BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE);
250:   } else BVDuplicateResize(nep->V,contour->npoints*ctx->L,&ctx->Y);

252:   if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) DSSetType(nep->ds,DSGNHEP);
253:   else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) DSSetType(nep->ds,DSNHEP);
254:   else {
255:     DSSetType(nep->ds,DSNEP);
256:     DSSetMethod(nep->ds,1);
257:     DSNEPSetRG(nep->ds,nep->rg);
258:     if (nep->fui==NEP_USER_INTERFACE_SPLIT) DSNEPSetFN(nep->ds,nep->nt,nep->f);
259:     else {
260:       PetscNew(&dsctxf);
261:       DSNEPSetComputeMatrixFunction(nep->ds,NEPContourDSComputeMatrix,dsctxf);
262:       dsctxf->nep = nep;
263:       ctx->dsctxf = dsctxf;
264:     }
265:   }
266:   DSAllocate(nep->ds,nep->ncv);
267:   nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
268:   NEPSetWorkVecs(nep,nwork);
269:   return 0;
270: }

272: PetscErrorCode NEPSolve_CISS(NEP nep)
273: {
274:   NEP_CISS         *ctx = (NEP_CISS*)nep->data;
275:   SlepcContourData contour = ctx->contour;
276:   Mat              X,M,E,T,P,J;
277:   BV               V;
278:   PetscInt         i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
279:   PetscScalar      *Mu,*H0,*H1,*rr,*temp,center;
280:   PetscReal        error,max_error,radius,rgscale,est_eig,eta;
281:   PetscBool        isellipse,*fl1;
282:   Vec              si;
283:   SlepcSC          sc;
284:   PetscRandom      rand;

286:   DSSetFromOptions(nep->ds);
287:   DSGetSlepcSC(nep->ds,&sc);
288:   sc->comparison    = SlepcCompareLargestMagnitude;
289:   sc->comparisonctx = NULL;
290:   sc->map           = NULL;
291:   sc->mapobj        = NULL;
292:   DSGetLeadingDimension(nep->ds,&ld);
293:   RGComputeQuadrature(nep->rg,RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight);
294:   if (contour->pA) {
295:     T = contour->pA[0];
296:     P = ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && contour->pP))? contour->pP[0]: T;
297:   } else {
298:     T = nep->function;
299:     P = nep->function_pre? nep->function_pre: nep->function;
300:   }
301:   NEPCISSSetUp(nep,T,P);
302:   BVSetActiveColumns(ctx->V,0,ctx->L);
303:   BVSetRandomSign(ctx->V);
304:   BVGetRandomContext(ctx->V,&rand);
305:   if (contour->pA) {
306:     J = ctx->J;
307:     V = ctx->pV;
308:     BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
309:   } else {
310:     J = nep->jacobian;
311:     V = ctx->V;
312:   }
313:   NEPCISSSolve(nep,J,V,0,ctx->L);
314:   PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
315:   if (isellipse) {
316:     BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig);
317:     PetscInfo(nep,"Estimated eigenvalue count: %f\n",(double)est_eig);
318:     eta = PetscPowReal(10.0,-PetscLog10Real(nep->tol)/ctx->N);
319:     L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
320:     if (L_add>ctx->L_max-ctx->L) {
321:       PetscInfo(nep,"Number of eigenvalues inside the contour path may be too large\n");
322:       L_add = ctx->L_max-ctx->L;
323:     }
324:   }
325:   /* Updates L after estimate the number of eigenvalue */
326:   if (L_add>0) {
327:     PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add);
328:     BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
329:     BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
330:     BVSetRandomSign(ctx->V);
331:     if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
332:     ctx->L += L_add;
333:     NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
334:   }

336:   PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
337:   for (i=0;i<ctx->refine_blocksize;i++) {
338:     BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
339:     CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
340:     PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
341:     SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
342:     PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
343:     if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
344:     L_add = L_base;
345:     if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
346:     PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add);
347:     BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
348:     BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
349:     BVSetRandomSign(ctx->V);
350:     if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
351:     ctx->L += L_add;
352:     NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
353:     if (L_add) {
354:       PetscFree2(Mu,H0);
355:       PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
356:     }
357:   }

359:   RGGetScale(nep->rg,&rgscale);
360:   RGEllipseGetParameters(nep->rg,&center,&radius,NULL);

362:   if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);

364:   while (nep->reason == NEP_CONVERGED_ITERATING) {
365:     nep->its++;
366:     for (inner=0;inner<=ctx->refine_inner;inner++) {
367:       if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
368:         BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
369:         CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
370:         PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
371:         SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
372:         PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
373:       } else {
374:         BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
375:         /* compute SVD of S */
376:         BVSVDAndRank(ctx->S,ctx->M,ctx->L,ctx->delta,(ctx->extraction==NEP_CISS_EXTRACTION_CAA)?BV_SVD_METHOD_QR_CAA:BV_SVD_METHOD_QR,H0,ctx->sigma,&nv);
377:       }
378:       PetscInfo(nep,"Estimated rank: nv = %" PetscInt_FMT "\n",nv);
379:       if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
380:         BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
381:         BVSetActiveColumns(ctx->S,0,ctx->L);
382:         BVSetActiveColumns(ctx->V,0,ctx->L);
383:         BVCopy(ctx->S,ctx->V);
384:         if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
385:         NEPCISSSolve(nep,J,V,0,ctx->L);
386:       } else break;
387:     }
388:     nep->nconv = 0;
389:     if (nv == 0) { nep->reason = NEP_CONVERGED_TOL; break; }
390:     else {
391:       /* Extracting eigenpairs */
392:       DSSetDimensions(nep->ds,nv,0,0);
393:       DSSetState(nep->ds,DS_STATE_RAW);
394:       if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
395:         CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
396:         CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1);
397:         DSGetArray(nep->ds,DS_MAT_A,&temp);
398:         for (j=0;j<nv;j++)
399:           for (i=0;i<nv;i++)
400:             temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
401:         DSRestoreArray(nep->ds,DS_MAT_A,&temp);
402:         DSGetArray(nep->ds,DS_MAT_B,&temp);
403:         for (j=0;j<nv;j++)
404:           for (i=0;i<nv;i++)
405:             temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
406:         DSRestoreArray(nep->ds,DS_MAT_B,&temp);
407:       } else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) {
408:         BVSetActiveColumns(ctx->S,0,nv);
409:         DSGetArray(nep->ds,DS_MAT_A,&temp);
410:         for (i=0;i<nv;i++) PetscArraycpy(temp+i*ld,H0+i*nv,nv);
411:         DSRestoreArray(nep->ds,DS_MAT_A,&temp);
412:       } else {
413:         BVSetActiveColumns(ctx->S,0,nv);
414:         if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
415:           for (i=0;i<nep->nt;i++) {
416:             DSGetMat(nep->ds,DSMatExtra[i],&E);
417:             BVMatProject(ctx->S,nep->A[i],ctx->S,E);
418:             DSRestoreMat(nep->ds,DSMatExtra[i],&E);
419:           }
420:         } else { ctx->dsctxf->Q = ctx->S; }
421:       }
422:       DSSolve(nep->ds,nep->eigr,nep->eigi);
423:       DSSynchronize(nep->ds,nep->eigr,nep->eigi);
424:       DSGetDimensions(nep->ds,NULL,NULL,NULL,&nv);
425:       if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
426:         for (i=0;i<nv;i++) {
427:           nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
428:         }
429:       }
430:       PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
431:       DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
432:       DSGetMat(nep->ds,DS_MAT_X,&X);
433:       SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1);
434:       DSRestoreMat(nep->ds,DS_MAT_X,&X);
435:       RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside);
436:       for (i=0;i<nv;i++) {
437:         if (fl1[i] && inside[i]>=0) {
438:           rr[i] = 1.0;
439:           nep->nconv++;
440:         } else rr[i] = 0.0;
441:       }
442:       DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv);
443:       DSSynchronize(nep->ds,nep->eigr,nep->eigi);
444:       if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
445:         for (i=0;i<nv;i++) nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
446:       }
447:       PetscFree3(fl1,inside,rr);
448:       BVSetActiveColumns(nep->V,0,nv);
449:       DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
450:       if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
451:         BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
452:         BVSetActiveColumns(ctx->S,0,nv);
453:         BVCopy(ctx->S,nep->V);
454:         DSGetMat(nep->ds,DS_MAT_X,&X);
455:         BVMultInPlace(ctx->S,X,0,nep->nconv);
456:         BVMultInPlace(nep->V,X,0,nep->nconv);
457:         DSRestoreMat(nep->ds,DS_MAT_X,&X);
458:       } else {
459:         DSGetMat(nep->ds,DS_MAT_X,&X);
460:         BVMultInPlace(ctx->S,X,0,nep->nconv);
461:         DSRestoreMat(nep->ds,DS_MAT_X,&X);
462:         BVSetActiveColumns(ctx->S,0,nep->nconv);
463:         BVCopy(ctx->S,nep->V);
464:       }
465:       max_error = 0.0;
466:       for (i=0;i<nep->nconv;i++) {
467:         BVGetColumn(nep->V,i,&si);
468:         VecNormalize(si,NULL);
469:         NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error);
470:         (*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx);
471:         BVRestoreColumn(nep->V,i,&si);
472:         max_error = PetscMax(max_error,error);
473:       }
474:       if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
475:       else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
476:       else {
477:         if (nep->nconv > ctx->L) nv = nep->nconv;
478:         else if (ctx->L > nv) nv = ctx->L;
479:         nv = PetscMin(nv,ctx->L*ctx->M);
480:         MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
481:         MatSetRandom(M,rand);
482:         BVSetActiveColumns(ctx->S,0,nv);
483:         BVMultInPlace(ctx->S,M,0,ctx->L);
484:         MatDestroy(&M);
485:         BVSetActiveColumns(ctx->S,0,ctx->L);
486:         BVSetActiveColumns(ctx->V,0,ctx->L);
487:         BVCopy(ctx->S,ctx->V);
488:         if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
489:         NEPCISSSolve(nep,J,V,0,ctx->L);
490:       }
491:     }
492:   }
493:   PetscFree2(Mu,H0);
494:   if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscFree(H1);
495:   return 0;
496: }

498: static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
499: {
500:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
501:   PetscInt       oN,oL,oM,oLmax,onpart;
502:   PetscMPIInt    size;

504:   oN = ctx->N;
505:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
506:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
507:   } else {
510:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
511:   }
512:   oL = ctx->L;
513:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
514:     ctx->L = 16;
515:   } else {
517:     ctx->L = bs;
518:   }
519:   oM = ctx->M;
520:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
521:     ctx->M = ctx->N/4;
522:   } else {
525:     ctx->M = PetscMax(ms,2);
526:   }
527:   onpart = ctx->npart;
528:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
529:     ctx->npart = 1;
530:   } else {
531:     MPI_Comm_size(PetscObjectComm((PetscObject)nep),&size);
533:     ctx->npart = npart;
534:   }
535:   oLmax = ctx->L_max;
536:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
537:     ctx->L_max = 64;
538:   } else {
540:     ctx->L_max = PetscMax(bsmax,ctx->L);
541:   }
542:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
543:     SlepcContourDataDestroy(&ctx->contour);
544:     PetscInfo(nep,"Resetting the contour data structure due to a change of parameters\n");
545:     nep->state = NEP_STATE_INITIAL;
546:   }
547:   ctx->isreal = realmats;
548:   if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) nep->state = NEP_STATE_INITIAL;
549:   return 0;
550: }

552: /*@
553:    NEPCISSSetSizes - Sets the values of various size parameters in the CISS solver.

555:    Logically Collective on nep

557:    Input Parameters:
558: +  nep   - the nonlinear eigensolver context
559: .  ip    - number of integration points
560: .  bs    - block size
561: .  ms    - moment size
562: .  npart - number of partitions when splitting the communicator
563: .  bsmax - max block size
564: -  realmats - T(z) is real for real z

566:    Options Database Keys:
567: +  -nep_ciss_integration_points - Sets the number of integration points
568: .  -nep_ciss_blocksize - Sets the block size
569: .  -nep_ciss_moments - Sets the moment size
570: .  -nep_ciss_partitions - Sets the number of partitions
571: .  -nep_ciss_maxblocksize - Sets the maximum block size
572: -  -nep_ciss_realmats - T(z) is real for real z

574:    Notes:
575:    The default number of partitions is 1. This means the internal KSP object is shared
576:    among all processes of the NEP communicator. Otherwise, the communicator is split
577:    into npart communicators, so that npart KSP solves proceed simultaneously.

579:    The realmats flag can be set to true when T(.) is guaranteed to be real
580:    when the argument is a real value, for example, when all matrices in
581:    the split form are real. When set to true, the solver avoids some computations.

583:    Level: advanced

585: .seealso: NEPCISSGetSizes()
586: @*/
587: PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
588: {
596:   PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
597:   return 0;
598: }

600: static PetscErrorCode NEPCISSGetSizes_CISS(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
601: {
602:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

604:   if (ip) *ip = ctx->N;
605:   if (bs) *bs = ctx->L;
606:   if (ms) *ms = ctx->M;
607:   if (npart) *npart = ctx->npart;
608:   if (bsmax) *bsmax = ctx->L_max;
609:   if (realmats) *realmats = ctx->isreal;
610:   return 0;
611: }

613: /*@
614:    NEPCISSGetSizes - Gets the values of various size parameters in the CISS solver.

616:    Not Collective

618:    Input Parameter:
619: .  nep - the nonlinear eigensolver context

621:    Output Parameters:
622: +  ip    - number of integration points
623: .  bs    - block size
624: .  ms    - moment size
625: .  npart - number of partitions when splitting the communicator
626: .  bsmax - max block size
627: -  realmats - T(z) is real for real z

629:    Level: advanced

631: .seealso: NEPCISSSetSizes()
632: @*/
633: PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
634: {
636:   PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
637:   return 0;
638: }

640: static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
641: {
642:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

644:   if (delta == PETSC_DEFAULT) {
645:     ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
646:   } else {
648:     ctx->delta = delta;
649:   }
650:   if (spur == PETSC_DEFAULT) {
651:     ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
652:   } else {
654:     ctx->spurious_threshold = spur;
655:   }
656:   return 0;
657: }

659: /*@
660:    NEPCISSSetThreshold - Sets the values of various threshold parameters in
661:    the CISS solver.

663:    Logically Collective on nep

665:    Input Parameters:
666: +  nep   - the nonlinear eigensolver context
667: .  delta - threshold for numerical rank
668: -  spur  - spurious threshold (to discard spurious eigenpairs)

670:    Options Database Keys:
671: +  -nep_ciss_delta - Sets the delta
672: -  -nep_ciss_spurious_threshold - Sets the spurious threshold

674:    Level: advanced

676: .seealso: NEPCISSGetThreshold()
677: @*/
678: PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
679: {
683:   PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
684:   return 0;
685: }

687: static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
688: {
689:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

691:   if (delta) *delta = ctx->delta;
692:   if (spur)  *spur = ctx->spurious_threshold;
693:   return 0;
694: }

696: /*@
697:    NEPCISSGetThreshold - Gets the values of various threshold parameters in
698:    the CISS solver.

700:    Not Collective

702:    Input Parameter:
703: .  nep - the nonlinear eigensolver context

705:    Output Parameters:
706: +  delta - threshold for numerical rank
707: -  spur  - spurious threshold (to discard spurious eigenpairs)

709:    Level: advanced

711: .seealso: NEPCISSSetThreshold()
712: @*/
713: PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
714: {
716:   PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
717:   return 0;
718: }

720: static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
721: {
722:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

724:   if (inner == PETSC_DEFAULT) {
725:     ctx->refine_inner = 0;
726:   } else {
728:     ctx->refine_inner = inner;
729:   }
730:   if (blsize == PETSC_DEFAULT) {
731:     ctx->refine_blocksize = 0;
732:   } else {
734:     ctx->refine_blocksize = blsize;
735:   }
736:   return 0;
737: }

739: /*@
740:    NEPCISSSetRefinement - Sets the values of various refinement parameters
741:    in the CISS solver.

743:    Logically Collective on nep

745:    Input Parameters:
746: +  nep    - the nonlinear eigensolver context
747: .  inner  - number of iterative refinement iterations (inner loop)
748: -  blsize - number of iterative refinement iterations (blocksize loop)

750:    Options Database Keys:
751: +  -nep_ciss_refine_inner - Sets number of inner iterations
752: -  -nep_ciss_refine_blocksize - Sets number of blocksize iterations

754:    Level: advanced

756: .seealso: NEPCISSGetRefinement()
757: @*/
758: PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
759: {
763:   PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
764:   return 0;
765: }

767: static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
768: {
769:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

771:   if (inner)  *inner = ctx->refine_inner;
772:   if (blsize) *blsize = ctx->refine_blocksize;
773:   return 0;
774: }

776: /*@
777:    NEPCISSGetRefinement - Gets the values of various refinement parameters
778:    in the CISS solver.

780:    Not Collective

782:    Input Parameter:
783: .  nep - the nonlinear eigensolver context

785:    Output Parameters:
786: +  inner  - number of iterative refinement iterations (inner loop)
787: -  blsize - number of iterative refinement iterations (blocksize loop)

789:    Level: advanced

791: .seealso: NEPCISSSetRefinement()
792: @*/
793: PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
794: {
796:   PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
797:   return 0;
798: }

800: static PetscErrorCode NEPCISSSetExtraction_CISS(NEP nep,NEPCISSExtraction extraction)
801: {
802:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

804:   if (ctx->extraction != extraction) {
805:     ctx->extraction = extraction;
806:     nep->state      = NEP_STATE_INITIAL;
807:   }
808:   return 0;
809: }

811: /*@
812:    NEPCISSSetExtraction - Sets the extraction technique used in the CISS solver.

814:    Logically Collective on nep

816:    Input Parameters:
817: +  nep        - the nonlinear eigensolver context
818: -  extraction - the extraction technique

820:    Options Database Key:
821: .  -nep_ciss_extraction - Sets the extraction technique (either 'ritz', 'hankel' or 'caa')

823:    Notes:
824:    By default, the Rayleigh-Ritz extraction is used (NEP_CISS_EXTRACTION_RITZ).

826:    If the 'hankel' or the 'caa' option is specified (NEP_CISS_EXTRACTION_HANKEL or
827:    NEP_CISS_EXTRACTION_CAA), then the Block Hankel method, or the Communication-avoiding
828:    Arnoldi method, respectively, is used for extracting eigenpairs.

830:    Level: advanced

832: .seealso: NEPCISSGetExtraction(), NEPCISSExtraction
833: @*/
834: PetscErrorCode NEPCISSSetExtraction(NEP nep,NEPCISSExtraction extraction)
835: {
838:   PetscTryMethod(nep,"NEPCISSSetExtraction_C",(NEP,NEPCISSExtraction),(nep,extraction));
839:   return 0;
840: }

842: static PetscErrorCode NEPCISSGetExtraction_CISS(NEP nep,NEPCISSExtraction *extraction)
843: {
844:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

846:   *extraction = ctx->extraction;
847:   return 0;
848: }

850: /*@
851:    NEPCISSGetExtraction - Gets the extraction technique used in the CISS solver.

853:    Not Collective

855:    Input Parameter:
856: .  nep - the nonlinear eigensolver context

858:    Output Parameters:
859: .  extraction - extraction technique

861:    Level: advanced

863: .seealso: NEPCISSSetExtraction() NEPCISSExtraction
864: @*/
865: PetscErrorCode NEPCISSGetExtraction(NEP nep,NEPCISSExtraction *extraction)
866: {
869:   PetscUseMethod(nep,"NEPCISSGetExtraction_C",(NEP,NEPCISSExtraction*),(nep,extraction));
870:   return 0;
871: }

873: static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
874: {
875:   NEP_CISS         *ctx = (NEP_CISS*)nep->data;
876:   SlepcContourData contour;
877:   PetscInt         i;
878:   PC               pc;
879:   MPI_Comm         child;

881:   if (!ctx->contour) {  /* initialize contour data structure first */
882:     RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj);
883:     SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour);
884:   }
885:   contour = ctx->contour;
886:   if (!contour->ksp) {
887:     PetscMalloc1(contour->npoints,&contour->ksp);
888:     PetscSubcommGetChild(contour->subcomm,&child);
889:     for (i=0;i<contour->npoints;i++) {
890:       KSPCreate(child,&contour->ksp[i]);
891:       PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)nep,1);
892:       KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)nep)->prefix);
893:       KSPAppendOptionsPrefix(contour->ksp[i],"nep_ciss_");
894:       PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)nep)->options);
895:       KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE);
896:       KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
897:       KSPGetPC(contour->ksp[i],&pc);
898:       if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
899:         KSPSetType(contour->ksp[i],KSPBCGS);
900:         PCSetType(pc,PCBJACOBI);
901:       } else {
902:         KSPSetType(contour->ksp[i],KSPPREONLY);
903:         PCSetType(pc,PCLU);
904:       }
905:     }
906:   }
907:   if (nsolve) *nsolve = contour->npoints;
908:   if (ksp)    *ksp    = contour->ksp;
909:   return 0;
910: }

912: /*@C
913:    NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
914:    the CISS solver.

916:    Not Collective

918:    Input Parameter:
919: .  nep - nonlinear eigenvalue solver

921:    Output Parameters:
922: +  nsolve - number of solver objects
923: -  ksp - array of linear solver object

925:    Notes:
926:    The number of KSP solvers is equal to the number of integration points divided by
927:    the number of partitions. This value is halved in the case of real matrices with
928:    a region centered at the real axis.

930:    Level: advanced

932: .seealso: NEPCISSSetSizes()
933: @*/
934: PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
935: {
937:   PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
938:   return 0;
939: }

941: PetscErrorCode NEPReset_CISS(NEP nep)
942: {
943:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

945:   BVDestroy(&ctx->S);
946:   BVDestroy(&ctx->V);
947:   BVDestroy(&ctx->Y);
948:   SlepcContourDataReset(ctx->contour);
949:   MatDestroy(&ctx->J);
950:   BVDestroy(&ctx->pV);
951:   if (ctx->extraction == NEP_CISS_EXTRACTION_RITZ && nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscFree(ctx->dsctxf);
952:   return 0;
953: }

955: PetscErrorCode NEPSetFromOptions_CISS(NEP nep,PetscOptionItems *PetscOptionsObject)
956: {
957:   NEP_CISS          *ctx = (NEP_CISS*)nep->data;
958:   PetscReal         r1,r2;
959:   PetscInt          i,i1,i2,i3,i4,i5,i6,i7;
960:   PetscBool         b1,flg,flg2,flg3,flg4,flg5,flg6;
961:   NEPCISSExtraction extraction;

963:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP CISS Options");

965:     NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1);
966:     PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,&flg);
967:     PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,&flg2);
968:     PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,&flg3);
969:     PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,&flg4);
970:     PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,&flg5);
971:     PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,&flg6);
972:     if (flg || flg2 || flg3 || flg4 || flg5 || flg6) NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1);

974:     NEPCISSGetThreshold(nep,&r1,&r2);
975:     PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,&flg);
976:     PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,&flg2);
977:     if (flg || flg2) NEPCISSSetThreshold(nep,r1,r2);

979:     NEPCISSGetRefinement(nep,&i6,&i7);
980:     PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,&flg);
981:     PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,&flg2);
982:     if (flg || flg2) NEPCISSSetRefinement(nep,i6,i7);

984:     PetscOptionsEnum("-nep_ciss_extraction","Extraction technique","NEPCISSSetExtraction",NEPCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
985:     if (flg) NEPCISSSetExtraction(nep,extraction);

987:   PetscOptionsHeadEnd();

989:   if (!nep->rg) NEPGetRG(nep,&nep->rg);
990:   RGSetFromOptions(nep->rg); /* this is necessary here to set useconj */
991:   if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
992:   PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
993:   for (i=0;i<ctx->contour->npoints;i++) KSPSetFromOptions(ctx->contour->ksp[i]);
994:   PetscSubcommSetFromOptions(ctx->contour->subcomm);
995:   return 0;
996: }

998: PetscErrorCode NEPDestroy_CISS(NEP nep)
999: {
1000:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

1002:   SlepcContourDataDestroy(&ctx->contour);
1003:   PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1004:   PetscFree(nep->data);
1005:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL);
1006:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL);
1007:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL);
1008:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL);
1009:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL);
1010:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL);
1011:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NULL);
1012:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NULL);
1013:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL);
1014:   return 0;
1015: }

1017: PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1018: {
1019:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
1020:   PetscBool      isascii;
1021:   PetscViewer    sviewer;

1023:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1024:   if (isascii) {
1025:     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);
1026:     if (ctx->isreal) PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n");
1027:     PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1028:     PetscViewerASCIIPrintf(viewer,"  iterative refinement  { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize);
1029:     PetscViewerASCIIPrintf(viewer,"  extraction: %s\n",NEPCISSExtractions[ctx->extraction]);
1030:     if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
1031:     PetscAssert(ctx->contour && ctx->contour->ksp,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Something went wrong with NEPCISSGetKSPs()");
1032:     PetscViewerASCIIPushTab(viewer);
1033:     if (ctx->npart>1 && ctx->contour->subcomm) {
1034:       PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1035:       if (!ctx->contour->subcomm->color) KSPView(ctx->contour->ksp[0],sviewer);
1036:       PetscViewerFlush(sviewer);
1037:       PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1038:       PetscViewerFlush(viewer);
1039:       /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1040:       PetscViewerASCIIPopSynchronized(viewer);
1041:     } else KSPView(ctx->contour->ksp[0],viewer);
1042:     PetscViewerASCIIPopTab(viewer);
1043:   }
1044:   return 0;
1045: }

1047: SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1048: {
1049:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

1051:   PetscNew(&ctx);
1052:   nep->data = ctx;
1053:   /* set default values of parameters */
1054:   ctx->N                  = 32;
1055:   ctx->L                  = 16;
1056:   ctx->M                  = ctx->N/4;
1057:   ctx->delta              = SLEPC_DEFAULT_TOL*1e-4;
1058:   ctx->L_max              = 64;
1059:   ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1060:   ctx->isreal             = PETSC_FALSE;
1061:   ctx->npart              = 1;

1063:   nep->useds = PETSC_TRUE;

1065:   nep->ops->solve          = NEPSolve_CISS;
1066:   nep->ops->setup          = NEPSetUp_CISS;
1067:   nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1068:   nep->ops->reset          = NEPReset_CISS;
1069:   nep->ops->destroy        = NEPDestroy_CISS;
1070:   nep->ops->view           = NEPView_CISS;

1072:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS);
1073:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS);
1074:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS);
1075:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS);
1076:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS);
1077:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS);
1078:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NEPCISSSetExtraction_CISS);
1079:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NEPCISSGetExtraction_CISS);
1080:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS);
1081:   return 0;
1082: }