Actual source code: nciss.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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/nepimpl.h>         /*I "slepcnep.h" I*/
 34: #include <slepcblaslapack.h>

 36: typedef struct {
 37:   /* parameters */
 38:   PetscInt     N;          /* number of integration points (32) */
 39:   PetscInt     L;          /* block size (16) */
 40:   PetscInt     M;          /* moment degree (N/4 = 4) */
 41:   PetscReal    delta;      /* threshold of singular value (1e-12) */
 42:   PetscInt     L_max;      /* maximum number of columns of the source matrix V */
 43:   PetscReal    spurious_threshold; /* discard spurious eigenpairs */
 44:   PetscBool    isreal;     /* T(z) is real for real z */
 45:   PetscInt     npart;      /* number of partitions */
 46:   PetscInt     refine_inner;
 47:   PetscInt     refine_blocksize;
 48:   /* private data */
 49:   PetscReal    *sigma;     /* threshold for numerical rank */
 50:   PetscInt     subcomm_id;
 51:   PetscInt     num_solve_point;
 52:   PetscScalar  *weight;
 53:   PetscScalar  *omega;
 54:   PetscScalar  *pp;
 55:   BV           V;
 56:   BV           S;
 57:   BV           Y;
 58:   KSP          *ksp;       /* ksp array for storing factorizations at integration points */
 59:   PetscBool    useconj;
 60:   PetscReal    est_eig;
 61:   PetscSubcomm subcomm;
 62: } NEP_CISS;

 64: /* destroy KSP objects when the number of solve points changes */
 65: PETSC_STATIC_INLINE PetscErrorCode NEPCISSResetSolvers(NEP nep)
 66: {
 68:   PetscInt       i;
 69:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

 72:   if (ctx->ksp) {
 73:     for (i=0;i<ctx->num_solve_point;i++) {
 74:       KSPDestroy(&ctx->ksp[i]);
 75:     }
 76:     PetscFree(ctx->ksp);
 77:   }
 78:   return(0);
 79: }

 81: /* clean PetscSubcomm object when the number of partitions changes */
 82: PETSC_STATIC_INLINE PetscErrorCode NEPCISSResetSubcomm(NEP nep)
 83: {
 85:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

 88:   NEPCISSResetSolvers(nep);
 89:   PetscSubcommDestroy(&ctx->subcomm);
 90:   return(0);
 91: }

 93: /* determine whether half of integration points can be avoided (use its conjugates);
 94:    depends on isreal and the center of the region */
 95: PETSC_STATIC_INLINE PetscErrorCode NEPCISSSetUseConj(NEP nep,PetscBool *useconj)
 96: {
 98:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
 99:   PetscScalar    center;
100:   PetscBool      isellipse=PETSC_FALSE;

103:   *useconj = PETSC_FALSE;
104:   PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
105:   if (isellipse) {
106:     RGEllipseGetParameters(nep->rg,&center,NULL,NULL);
107:     *useconj = (ctx->isreal && PetscImaginaryPart(center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
108:   }
109:   return(0);
110: }

112: /* create PetscSubcomm object and determine num_solve_point (depends on npart, N, useconj) */
113: PETSC_STATIC_INLINE PetscErrorCode NEPCISSSetUpSubComm(NEP nep,PetscInt *num_solve_point)
114: {
116:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
117:   PetscInt       N = ctx->N;

120:   PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&ctx->subcomm);
121:   PetscSubcommSetNumber(ctx->subcomm,ctx->npart);
122:   PetscSubcommSetType(ctx->subcomm,PETSC_SUBCOMM_INTERLACED);
123:   PetscLogObjectMemory((PetscObject)nep,sizeof(PetscSubcomm));
124:   ctx->subcomm_id = ctx->subcomm->color;
125:   NEPCISSSetUseConj(nep,&ctx->useconj);
126:   if (ctx->useconj) N = N/2;
127:   *num_solve_point = N / ctx->npart;
128:   if (N%ctx->npart > ctx->subcomm_id) (*num_solve_point)++;
129:   return(0);
130: }

132: static PetscErrorCode SetPathParameter(NEP nep)
133: {
135:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
136:   PetscInt       i;
137:   PetscScalar    center;
138:   PetscReal      theta,radius,vscale,rgscale;
139:   PetscBool      isellipse=PETSC_FALSE;

142:   PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
143:   RGGetScale(nep->rg,&rgscale);
144:   if (isellipse) {
145:     RGEllipseGetParameters(nep->rg,&center,&radius,&vscale);
146:   } else SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Region must be Ellipse");
147:   for (i=0;i<ctx->N;i++) {
148:     theta = ((2*PETSC_PI)/ctx->N)*(i+0.5);
149:     ctx->pp[i] = PetscCosReal(theta) + PETSC_i*vscale*PetscSinReal(theta);
150:     ctx->weight[i] = radius*(vscale*PetscCosReal(theta) + PETSC_i*PetscSinReal(theta))/(PetscReal)ctx->N;
151:     ctx->omega[i] = rgscale*(center + radius*ctx->pp[i]);
152:   }
153:   return(0);
154: }

156: static PetscErrorCode CISSVecSetRandom(BV V,PetscInt i0,PetscInt i1)
157: {
159:   PetscInt       i,j,nlocal;
160:   PetscScalar    *vdata;
161:   Vec            x;

164:   BVGetSizes(V,&nlocal,NULL,NULL);
165:   for (i=i0;i<i1;i++) {
166:     BVSetRandomColumn(V,i);
167:     BVGetColumn(V,i,&x);
168:     VecGetArray(x,&vdata);
169:     for (j=0;j<nlocal;j++) {
170:       vdata[j] = PetscRealPart(vdata[j]);
171:       if (PetscRealPart(vdata[j]) < 0.5) vdata[j] = -1.0;
172:       else vdata[j] = 1.0;
173:     }
174:     VecRestoreArray(x,&vdata);
175:     BVRestoreColumn(V,i,&x);
176:   }
177:   return(0);
178: }

180: static PetscErrorCode SolveLinearSystem(NEP nep,Mat T,Mat dT,BV V,PetscInt L_start,PetscInt L_end,PetscBool initksp)
181: {
183:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
184:   PetscInt       i,j,p_id;
185:   Mat            kspMat;
186:   Vec            Bvj,vj,yj;

189:   if (!ctx->ksp) { NEPCISSGetKSPs(nep,&ctx->num_solve_point,&ctx->ksp); }
190:   BVCreateVec(V,&Bvj);
191:   for (i=0;i<ctx->num_solve_point;i++) {
192:     p_id = i*ctx->subcomm->n + ctx->subcomm_id;
193:     if (initksp) {
194:       NEPComputeFunction(nep,ctx->omega[p_id],T,T);
195:       MatDuplicate(T,MAT_COPY_VALUES,&kspMat);
196:       KSPSetOperators(ctx->ksp[i],kspMat,kspMat);
197:       MatDestroy(&kspMat);
198:     }
199:     NEPComputeJacobian(nep,ctx->omega[p_id],dT);
200:     for (j=L_start;j<L_end;j++) {
201:       BVGetColumn(V,j,&vj);
202:       BVGetColumn(ctx->Y,i*ctx->L_max+j,&yj);
203:       MatMult(dT,vj,Bvj);
204:       KSPSolve(ctx->ksp[i],Bvj,yj);
205:       BVRestoreColumn(V,j,&vj);
206:       BVRestoreColumn(ctx->Y,i*ctx->L_max+j,&yj);
207:     }
208:   }
209:   VecDestroy(&Bvj);
210:   return(0);
211: }

213: static PetscErrorCode EstimateNumberEigs(NEP nep,PetscInt *L_add)
214: {
216:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
217:   PetscInt       i,j,p_id;
218:   PetscScalar    tmp,m = 1,sum = 0.0;
219:   PetscReal      eta;
220:   Vec            v,vtemp,vj,yj;

223:   BVGetColumn(ctx->Y,0,&yj);
224:   VecDuplicate(yj,&v);
225:   BVRestoreColumn(ctx->Y,0,&yj);
226:   BVCreateVec(ctx->V,&vtemp);
227:   for (j=0;j<ctx->L;j++) {
228:     VecSet(v,0);
229:     for (i=0;i<ctx->num_solve_point; i++) {
230:       p_id = i*ctx->subcomm->n + ctx->subcomm_id;
231:       BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
232:       BVMultVec(ctx->Y,ctx->weight[p_id],1,v,&m);
233:     }
234:     BVGetColumn(ctx->V,j,&vj);
235:     VecDot(vj,v,&tmp);
236:     BVRestoreColumn(ctx->V,j,&vj);
237:     if (ctx->useconj) sum += PetscRealPart(tmp)*2;
238:     else sum += tmp;
239:   }
240:   ctx->est_eig = PetscAbsScalar(sum/(PetscReal)ctx->L);
241:   eta = PetscPowReal(10,-PetscLog10Real(nep->tol)/ctx->N);
242:   PetscInfo1(nep,"Estimation_#Eig %f\n",(double)ctx->est_eig);
243:   *L_add = (PetscInt)PetscCeilReal((ctx->est_eig*eta)/ctx->M) - ctx->L;
244:   if (*L_add < 0) *L_add = 0;
245:   if (*L_add>ctx->L_max-ctx->L) {
246:     PetscInfo(nep,"Number of eigenvalues around the contour path may be too large\n");
247:     *L_add = ctx->L_max-ctx->L;
248:   }
249:   VecDestroy(&v);
250:   VecDestroy(&vtemp);
251:   return(0);
252: }

254: static PetscErrorCode CalcMu(NEP nep, PetscScalar *Mu)
255: {
257:   PetscMPIInt    sub_size,len;
258:   PetscInt       i,j,k,s;
259:   PetscScalar    *m,*temp,*temp2,*ppk,alp;
260:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
261:   Mat            M;

264:   MPI_Comm_size(PetscSubcommChild(ctx->subcomm),&sub_size);
265:   PetscMalloc3(ctx->num_solve_point*ctx->L*(ctx->L+1),&temp,2*ctx->M*ctx->L*ctx->L,&temp2,ctx->num_solve_point,&ppk);
266:   MatCreateSeqDense(PETSC_COMM_SELF,ctx->L,ctx->L_max*ctx->num_solve_point,NULL,&M);
267:   for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] = 0;
268:   BVSetActiveColumns(ctx->Y,0,ctx->L_max*ctx->num_solve_point);
269:   BVSetActiveColumns(ctx->V,0,ctx->L);
270:   BVDot(ctx->Y,ctx->V,M);
271:   MatDenseGetArray(M,&m);
272:   for (i=0;i<ctx->num_solve_point;i++) {
273:     for (j=0;j<ctx->L;j++) {
274:       for (k=0;k<ctx->L;k++) {
275:         temp[k+j*ctx->L+i*ctx->L*ctx->L]=m[k+j*ctx->L+i*ctx->L*ctx->L_max];
276:       }
277:     }
278:   }
279:   MatDenseRestoreArray(M,&m);
280:   for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
281:   for (k=0;k<2*ctx->M;k++) {
282:     for (j=0;j<ctx->L;j++) {
283:       for (i=0;i<ctx->num_solve_point;i++) {
284:         alp = ppk[i]*ctx->weight[i*ctx->subcomm->n + ctx->subcomm_id];
285:         for (s=0;s<ctx->L;s++) {
286:           if (ctx->useconj) temp2[s+(j+k*ctx->L)*ctx->L] += PetscRealPart(alp*temp[s+(j+i*ctx->L)*ctx->L])*2;
287:           else temp2[s+(j+k*ctx->L)*ctx->L] += alp*temp[s+(j+i*ctx->L)*ctx->L];
288:         }
289:       }
290:     }
291:     for (i=0;i<ctx->num_solve_point;i++)
292:       ppk[i] *= ctx->pp[i*ctx->subcomm->n + ctx->subcomm_id];
293:   }
294:   for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] /= sub_size;
295:   PetscMPIIntCast(2*ctx->M*ctx->L*ctx->L,&len);
296:   MPI_Allreduce(temp2,Mu,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)nep));
297:   PetscFree3(temp,temp2,ppk);
298:   MatDestroy(&M);
299:   return(0);
300: }

302: static PetscErrorCode BlockHankel(NEP nep,PetscScalar *Mu,PetscInt s,PetscScalar *H)
303: {
304:   NEP_CISS *ctx = (NEP_CISS*)nep->data;
305:   PetscInt  i,j,k,L=ctx->L,M=ctx->M;

308:   for (k=0;k<L*M;k++)
309:     for (j=0;j<M;j++)
310:       for (i=0;i<L;i++)
311:         H[j*L+i+k*L*M] = Mu[i+k*L+(j+s)*L*L];
312:   return(0);
313: }

315: static PetscErrorCode SVD_H0(NEP nep,PetscScalar *S,PetscInt *K)
316: {
317: #if defined(PETSC_MISSING_LAPACK_GESVD)
319:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable");
320: #else
322:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
323:   PetscInt       i,ml=ctx->L*ctx->M;
324:   PetscBLASInt   m,n,lda,ldu,ldvt,lwork,info;
325:   PetscScalar    *work;
326: #if defined(PETSC_USE_COMPLEX)
327:   PetscReal      *rwork;
328: #endif

331:   PetscMalloc1(5*ml,&work);
332: #if defined(PETSC_USE_COMPLEX)
333:   PetscMalloc1(5*ml,&rwork);
334: #endif
335:   PetscBLASIntCast(ml,&m);
336:   n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
337:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
338: #if defined(PETSC_USE_COMPLEX)
339:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
340: #else
341:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
342: #endif
343:   SlepcCheckLapackInfo("gesvd",info);
344:   PetscFPTrapPop();
345:   (*K) = 0;
346:   for (i=0;i<ml;i++) {
347:     if (ctx->sigma[i]/PetscMax(ctx->sigma[0],1)>ctx->delta) (*K)++;
348:   }
349:   PetscFree(work);
350: #if defined(PETSC_USE_COMPLEX)
351:   PetscFree(rwork);
352: #endif
353:   return(0);
354: #endif
355: }

357: static PetscErrorCode ConstructS(NEP nep)
358: {
360:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
361:   PetscInt       i,j,k,vec_local_size,p_id;
362:   Vec            v,sj,yj;
363:   PetscScalar    *ppk, *v_data, m = 1;

366:   BVGetSizes(ctx->Y,&vec_local_size,NULL,NULL);
367:   PetscMalloc1(ctx->num_solve_point,&ppk);
368:   for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
369:   BVGetColumn(ctx->Y,0,&yj);
370:   VecDuplicate(yj,&v);
371:   BVRestoreColumn(ctx->Y,0,&yj);
372:   for (k=0;k<ctx->M;k++) {
373:     for (j=0;j<ctx->L;j++) {
374:       VecSet(v,0);
375:       for (i=0;i<ctx->num_solve_point;i++) {
376:         p_id = i*ctx->subcomm->n + ctx->subcomm_id;
377:         BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
378:         BVMultVec(ctx->Y,ppk[i]*ctx->weight[p_id],1,v,&m);
379:       }
380:       if (ctx->useconj) {
381:         VecGetArray(v,&v_data);
382:         for (i=0;i<vec_local_size;i++) v_data[i] = PetscRealPart(v_data[i])*2;
383:         VecRestoreArray(v,&v_data);
384:       }
385:       BVGetColumn(ctx->S,k*ctx->L+j,&sj);
386:       VecCopy(v,sj);
387:       BVRestoreColumn(ctx->S,k*ctx->L+j,&sj);
388:     }
389:     for (i=0;i<ctx->num_solve_point;i++) {
390:       p_id = i*ctx->subcomm->n + ctx->subcomm_id;
391:       ppk[i] *= ctx->pp[p_id];
392:     }
393:   }
394:   PetscFree(ppk);
395:   VecDestroy(&v);
396:   return(0);
397: }

399: static PetscErrorCode isGhost(NEP nep,PetscInt ld,PetscInt nv,PetscBool *fl)
400: {
402:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
403:   PetscInt       i,j;
404:   PetscScalar    *pX;
405:   PetscReal      *tau,s1,s2,tau_max=0.0;

408:   PetscMalloc1(nv,&tau);
409:   DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
410:   DSGetArray(nep->ds,DS_MAT_X,&pX);

412:   for (i=0;i<nv;i++) {
413:     s1 = 0;
414:     s2 = 0;
415:     for (j=0;j<nv;j++) {
416:       s1 += PetscAbsScalar(PetscPowScalar(pX[i*ld+j],2));
417:       s2 += PetscPowReal(PetscAbsScalar(pX[i*ld+j]),2)/ctx->sigma[j];
418:     }
419:     tau[i] = s1/s2;
420:     tau_max = PetscMax(tau_max,tau[i]);
421:   }
422:   DSRestoreArray(nep->ds,DS_MAT_X,&pX);
423:   for (i=0;i<nv;i++) {
424:     tau[i] /= tau_max;
425:   }
426:   for (i=0;i<nv;i++) {
427:     if (tau[i]>=ctx->spurious_threshold) fl[i] = PETSC_TRUE;
428:     else fl[i] = PETSC_FALSE;
429:   }
430:   PetscFree(tau);
431:   return(0);
432: }

434: PetscErrorCode NEPSetUp_CISS(NEP nep)
435: {
437:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
438:   PetscInt       nwork;
439:   PetscBool      istrivial,isellipse,flg,useconj;

442:   if (!nep->ncv) {
443:     nep->ncv = ctx->L_max*ctx->M;
444:     if (nep->ncv>nep->n) {
445:       nep->ncv = nep->n;
446:       ctx->L_max = nep->ncv/ctx->M;
447:       if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
448:     }
449:   } else {
450:     ctx->L_max = nep->ncv/ctx->M;
451:     if (!ctx->L_max) {
452:       ctx->L_max = 1;
453:       nep->ncv = ctx->L_max*ctx->M;
454:     }
455:   }
456:   ctx->L = PetscMin(ctx->L,ctx->L_max);
457:   if (!nep->max_it) nep->max_it = 1;
458:   if (!nep->mpd) nep->mpd = nep->ncv;
459:   if (!nep->which) nep->which = NEP_ALL;
460:   if (nep->stopping!=NEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support user-defined stopping test");

462:   /* check region */
463:   RGIsTrivial(nep->rg,&istrivial);
464:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
465:   RGGetComplement(nep->rg,&flg);
466:   if (flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
467:   PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
468:   if (!isellipse) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Currently only implemented for elliptic regions");
469:   /* if useconj has changed, then reset subcomm data */
470:   NEPCISSSetUseConj(nep,&useconj);
471:   if (useconj!=ctx->useconj) { NEPCISSResetSubcomm(nep); }

473:   /* create split comm */
474:   if (!ctx->subcomm) { NEPCISSSetUpSubComm(nep,&ctx->num_solve_point); }

476:   NEPAllocateSolution(nep,0);
477:   if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
478:   PetscMalloc4(ctx->N,&ctx->weight,ctx->N,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
479:   PetscLogObjectMemory((PetscObject)nep,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));

481:   /* allocate basis vectors */
482:   BVDestroy(&ctx->S);
483:   BVDuplicateResize(nep->V,ctx->L_max*ctx->M,&ctx->S);
484:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->S);
485:   BVDestroy(&ctx->V);
486:   BVDuplicateResize(nep->V,ctx->L_max,&ctx->V);
487:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->V);

489:   BVDestroy(&ctx->Y);
490:   BVDuplicateResize(nep->V,ctx->num_solve_point*ctx->L_max,&ctx->Y);

492:   DSSetType(nep->ds,DSGNHEP);
493:   DSAllocate(nep->ds,nep->ncv);
494:   nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
495:   NEPSetWorkVecs(nep,nwork);
496:   return(0);
497: }

499: PetscErrorCode NEPSolve_CISS(NEP nep)
500: {
502:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
503:   Mat            X,M;
504:   PetscInt       i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
505:   PetscScalar    *Mu,*H0,*H1,*rr,*temp,center;
506:   PetscReal      error,max_error,radius,rgscale;
507:   PetscBool      *fl1;
508:   Vec            si;
509:   SlepcSC        sc;
510:   PetscRandom    rand;

513:   DSGetSlepcSC(nep->ds,&sc);
514:   sc->comparison    = SlepcCompareLargestMagnitude;
515:   sc->comparisonctx = NULL;
516:   sc->map           = NULL;
517:   sc->mapobj        = NULL;
518:   DSGetLeadingDimension(nep->ds,&ld);
519:   SetPathParameter(nep);
520:   CISSVecSetRandom(ctx->V,0,ctx->L);
521:   BVGetRandomContext(ctx->V,&rand);

523:   SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_TRUE);
524:   EstimateNumberEigs(nep,&L_add);
525:   if (L_add>0) {
526:     PetscInfo2(nep,"Changing L %D -> %D by Estimate #Eig\n",ctx->L,ctx->L+L_add);
527:     CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
528:     SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
529:     ctx->L += L_add;
530:   }
531:   PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
532:   for (i=0;i<ctx->refine_blocksize;i++) {
533:     CalcMu(nep,Mu);
534:     BlockHankel(nep,Mu,0,H0);
535:     SVD_H0(nep,H0,&nv);
536:     if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
537:     L_add = L_base;
538:     if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
539:     PetscInfo2(nep,"Changing L %D -> %D by SVD(H0)\n",ctx->L,ctx->L+L_add);
540:     CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
541:     SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
542:     ctx->L += L_add;
543:   }
544:   PetscFree2(Mu,H0);

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

549:   PetscMalloc3(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0,ctx->L*ctx->M*ctx->L*ctx->M,&H1);
550:   while (nep->reason == NEP_CONVERGED_ITERATING) {
551:     nep->its++;
552:     for (inner=0;inner<=ctx->refine_inner;inner++) {
553:       CalcMu(nep,Mu);
554:       BlockHankel(nep,Mu,0,H0);
555:       SVD_H0(nep,H0,&nv);
556:       if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
557:         ConstructS(nep);
558:         BVSetActiveColumns(ctx->S,0,ctx->L);
559:         BVCopy(ctx->S,ctx->V);
560:         SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_FALSE);
561:       } else break;
562:     }

564:     nep->nconv = 0;
565:     if (nv == 0) break;
566:     BlockHankel(nep,Mu,0,H0);
567:     BlockHankel(nep,Mu,1,H1);
568:     DSSetDimensions(nep->ds,nv,0,0,0);
569:     DSSetState(nep->ds,DS_STATE_RAW);
570:     DSGetArray(nep->ds,DS_MAT_A,&temp);
571:     for (j=0;j<nv;j++)
572:       for (i=0;i<nv;i++)
573:         temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
574:     DSRestoreArray(nep->ds,DS_MAT_A,&temp);
575:     DSGetArray(nep->ds,DS_MAT_B,&temp);
576:     for (j=0;j<nv;j++)
577:       for (i=0;i<nv;i++)
578:         temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
579:     DSRestoreArray(nep->ds,DS_MAT_B,&temp);
580:     DSSolve(nep->ds,nep->eigr,nep->eigi);
581:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);
582:     for (i=0;i<nv;i++){
583:       nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
584: #if !defined(PETSC_USE_COMPLEX)
585:       nep->eigi[i] = nep->eigi[i]*radius*rgscale;
586: #endif
587:     }
588:     PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
589:     isGhost(nep,ld,nv,fl1);
590:     RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside);
591:     for (i=0;i<nv;i++) {
592:       if (fl1[i] && inside[i]>0) {
593:         rr[i] = 1.0;
594:         nep->nconv++;
595:       } else rr[i] = 0.0;
596:     }
597:     DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv);
598:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);
599:     for (i=0;i<nv;i++){
600:       nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
601: #if !defined(PETSC_USE_COMPLEX)
602:       nep->eigi[i] = nep->eigi[i]*radius*rgscale;
603: #endif
604:     }
605:     PetscFree3(fl1,inside,rr);
606:     BVSetActiveColumns(nep->V,0,nv);
607:     ConstructS(nep);
608:     BVSetActiveColumns(ctx->S,0,nv);
609:     BVCopy(ctx->S,nep->V);

611:     DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
612:     DSGetMat(nep->ds,DS_MAT_X,&X);
613:     BVMultInPlace(ctx->S,X,0,nep->nconv);
614:     BVMultInPlace(nep->V,X,0,nep->nconv);
615:     MatDestroy(&X);
616:     max_error = 0.0;
617:     for (i=0;i<nep->nconv;i++) {
618:       BVGetColumn(nep->V,i,&si);
619:       VecNormalize(si,NULL);
620:       NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error);
621:       (*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx);
622:       BVRestoreColumn(nep->V,i,&si);
623:       max_error = PetscMax(max_error,error);
624:     }
625:     if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
626:     else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
627:     else {
628:       if (nep->nconv > ctx->L) nv = nep->nconv;
629:       else if (ctx->L > nv) nv = ctx->L;
630:       MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
631:       MatDenseGetArray(M,&temp);
632:       for (i=0;i<ctx->L*nv;i++) {
633:         PetscRandomGetValue(rand,&temp[i]);
634:         temp[i] = PetscRealPart(temp[i]);
635:       }
636:       MatDenseRestoreArray(M,&temp);
637:       BVSetActiveColumns(ctx->S,0,nv);
638:       BVMultInPlace(ctx->S,M,0,ctx->L);
639:       MatDestroy(&M);
640:       BVSetActiveColumns(ctx->S,0,ctx->L);
641:       BVCopy(ctx->S,ctx->V);
642:       SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_FALSE);
643:     }
644:   }
645:   PetscFree3(Mu,H0,H1);
646:   return(0);
647: }

649: static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
650: {
652:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
653:   PetscInt       oN,onpart;

656:   oN = ctx->N;
657:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
658:     if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
659:   } else {
660:     if (ip<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
661:     if (ip%2) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
662:     if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
663:   }
664:   if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
665:     ctx->L = 16;
666:   } else {
667:     if (bs<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
668:     ctx->L = bs;
669:   }
670:   if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
671:     ctx->M = ctx->N/4;
672:   } else {
673:     if (ms<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
674:     if (ms>ctx->N) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
675:     ctx->M = ms;
676:   }
677:   onpart = ctx->npart;
678:   if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
679:     ctx->npart = 1;
680:   } else {
681:     if (npart<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
682:     ctx->npart = npart;
683:     if (npart>1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Current implementation does not allow partitions");
684:   }
685:   if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
686:     ctx->L_max = 64;
687:   } else {
688:     if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
689:     ctx->L_max = PetscMax(bsmax,ctx->L);
690:   }
691:   if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) { NEPCISSResetSubcomm(nep); }
692:   ctx->isreal = realmats;
693:   nep->state = NEP_STATE_INITIAL;
694:   return(0);
695: }

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

700:    Logically Collective on NEP

702:    Input Parameters:
703: +  nep   - the eigenproblem solver context
704: .  ip    - number of integration points
705: .  bs    - block size
706: .  ms    - moment size
707: .  npart - number of partitions when splitting the communicator
708: .  bsmax - max block size
709: -  realmats - T(z) is real for real z

711:    Options Database Keys:
712: +  -nep_ciss_integration_points - Sets the number of integration points
713: .  -nep_ciss_blocksize - Sets the block size
714: .  -nep_ciss_moments - Sets the moment size
715: .  -nep_ciss_partitions - Sets the number of partitions
716: .  -nep_ciss_maxblocksize - Sets the maximum block size
717: -  -nep_ciss_realmats - T(z) is real for real z

719:    Note:
720:    The default number of partitions is 1. This means the internal KSP object is shared
721:    among all processes of the NEP communicator. Otherwise, the communicator is split
722:    into npart communicators, so that npart KSP solves proceed simultaneously.

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

728:    Level: advanced

730: .seealso: NEPCISSGetSizes()
731: @*/
732: PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
733: {

744:   PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
745:   return(0);
746: }

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

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

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

765:    Not Collective

767:    Input Parameter:
768: .  nep - the eigenproblem solver context

770:    Output Parameters:
771: +  ip    - number of integration points
772: .  bs    - block size
773: .  ms    - moment size
774: .  npart - number of partitions when splitting the communicator
775: .  bsmax - max block size
776: -  realmats - T(z) is real for real z

778:    Level: advanced

780: .seealso: NEPCISSSetSizes()
781: @*/
782: PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
783: {

788:   PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
789:   return(0);
790: }

792: static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
793: {
794:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

797:   if (delta == PETSC_DEFAULT) {
798:     ctx->delta = 1e-12;
799:   } else {
800:     if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
801:     ctx->delta = delta;
802:   }
803:   if (spur == PETSC_DEFAULT) {
804:     ctx->spurious_threshold = 1e-4;
805:   } else {
806:     if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
807:     ctx->spurious_threshold = spur;
808:   }
809:   return(0);
810: }

812: /*@
813:    NEPCISSSetThreshold - Sets the values of various threshold parameters in
814:    the CISS solver.

816:    Logically Collective on NEP

818:    Input Parameters:
819: +  nep   - the eigenproblem solver context
820: .  delta - threshold for numerical rank
821: -  spur  - spurious threshold (to discard spurious eigenpairs)

823:    Options Database Keys:
824: +  -nep_ciss_delta - Sets the delta
825: -  -nep_ciss_spurious_threshold - Sets the spurious threshold

827:    Level: advanced

829: .seealso: NEPCISSGetThreshold()
830: @*/
831: PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
832: {

839:   PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
840:   return(0);
841: }

843: static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
844: {
845:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

848:   if (delta) *delta = ctx->delta;
849:   if (spur)  *spur = ctx->spurious_threshold;
850:   return(0);
851: }

853: /*@
854:    NEPCISSGetThreshold - Gets the values of various threshold parameters in
855:    the CISS solver.

857:    Not Collective

859:    Input Parameter:
860: .  nep - the eigenproblem solver context

862:    Output Parameters:
863: +  delta - threshold for numerical rank
864: -  spur  - spurious threshold (to discard spurious eigenpairs)

866:    Level: advanced

868: .seealso: NEPCISSSetThreshold()
869: @*/
870: PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
871: {

876:   PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
877:   return(0);
878: }

880: static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
881: {
882:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

885:   if (inner == PETSC_DEFAULT) {
886:     ctx->refine_inner = 0;
887:   } else {
888:     if (inner<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
889:     ctx->refine_inner = inner;
890:   }
891:   if (blsize == PETSC_DEFAULT) {
892:     ctx->refine_blocksize = 0;
893:   } else {
894:     if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
895:     ctx->refine_blocksize = blsize;
896:   }
897:   return(0);
898: }

900: /*@
901:    NEPCISSSetRefinement - Sets the values of various refinement parameters
902:    in the CISS solver.

904:    Logically Collective on NEP

906:    Input Parameters:
907: +  nep    - the eigenproblem solver context
908: .  inner  - number of iterative refinement iterations (inner loop)
909: -  blsize - number of iterative refinement iterations (blocksize loop)

911:    Options Database Keys:
912: +  -nep_ciss_refine_inner - Sets number of inner iterations
913: -  -nep_ciss_refine_blocksize - Sets number of blocksize iterations

915:    Level: advanced

917: .seealso: NEPCISSGetRefinement()
918: @*/
919: PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
920: {

927:   PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
928:   return(0);
929: }

931: static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
932: {
933:   NEP_CISS *ctx = (NEP_CISS*)nep->data;

936:   if (inner)  *inner = ctx->refine_inner;
937:   if (blsize) *blsize = ctx->refine_blocksize;
938:   return(0);
939: }

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

945:    Not Collective

947:    Input Parameter:
948: .  nep - the eigenproblem solver context

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

954:    Level: advanced

956: .seealso: NEPCISSSetRefinement()
957: @*/
958: PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
959: {

964:   PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
965:   return(0);
966: }

968: static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
969: {
971:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
972:   PetscInt       i;
973:   PC             pc;

976:   if (!ctx->ksp) {
977:     if (!ctx->subcomm) {  /* initialize subcomm first */
978:       NEPCISSSetUseConj(nep,&ctx->useconj);
979:       NEPCISSSetUpSubComm(nep,&ctx->num_solve_point);
980:     }
981:     PetscMalloc1(ctx->num_solve_point,&ctx->ksp);
982:     for (i=0;i<ctx->num_solve_point;i++) {
983:       KSPCreate(PetscSubcommChild(ctx->subcomm),&ctx->ksp[i]);
984:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
985:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
986:       KSPAppendOptionsPrefix(ctx->ksp[i],"nep_ciss_");
987:       PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
988:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
989:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
990:       KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
991:       KSPGetPC(ctx->ksp[i],&pc);
992:       KSPSetType(ctx->ksp[i],KSPPREONLY);
993:       PCSetType(pc,PCLU);
994:     }
995:   }
996:   if (nsolve) *nsolve = ctx->num_solve_point;
997:   if (ksp)    *ksp    = ctx->ksp;
998:   return(0);
999: }

1001: /*@C
1002:    NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
1003:    the CISS solver.

1005:    Not Collective

1007:    Input Parameter:
1008: .  nep - nonlinear eigenvalue solver

1010:    Output Parameters:
1011: +  nsolve - number of solver objects
1012: -  ksp - array of linear solver object

1014:    Notes:
1015:    The number of KSP solvers is equal to the number of integration points divided by
1016:    the number of partitions. This value is halved in the case of real matrices with
1017:    a region centered at the real axis.

1019:    Level: advanced

1021: .seealso: NEPCISSSetSizes()
1022: @*/
1023: PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1024: {

1029:   PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1030:   return(0);
1031: }

1033: PetscErrorCode NEPReset_CISS(NEP nep)
1034: {
1036:   PetscInt       i;
1037:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

1040:   BVDestroy(&ctx->S);
1041:   BVDestroy(&ctx->V);
1042:   BVDestroy(&ctx->Y);
1043:   for (i=0;i<ctx->num_solve_point;i++) {
1044:     KSPReset(ctx->ksp[i]);
1045:   }
1046:   return(0);
1047: }

1049: PetscErrorCode NEPSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,NEP nep)
1050: {
1052:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
1053:   PetscReal      r1,r2;
1054:   PetscInt       i,i1,i2,i3,i4,i5,i6,i7;
1055:   PetscBool      b1;

1058:   PetscOptionsHead(PetscOptionsObject,"NEP CISS Options");

1060:     NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1);
1061:     PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,NULL);
1062:     PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,NULL);
1063:     PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,NULL);
1064:     PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,NULL);
1065:     PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,NULL);
1066:     PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,NULL);
1067:     NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1);

1069:     NEPCISSGetThreshold(nep,&r1,&r2);
1070:     PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,NULL);
1071:     PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,NULL);
1072:     NEPCISSSetThreshold(nep,r1,r2);

1074:     NEPCISSGetRefinement(nep,&i6,&i7);
1075:     PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,NULL);
1076:     PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,NULL);
1077:     NEPCISSSetRefinement(nep,i6,i7);

1079:   PetscOptionsTail();

1081:   if (!nep->rg) { NEPGetRG(nep,&nep->rg); }
1082:   RGSetFromOptions(nep->rg); /* this is necessary here to set useconj */
1083:   if (!ctx->ksp) { NEPCISSGetKSPs(nep,&ctx->num_solve_point,&ctx->ksp); }
1084:   for (i=0;i<ctx->num_solve_point;i++) {
1085:     KSPSetFromOptions(ctx->ksp[i]);
1086:   }
1087:   PetscSubcommSetFromOptions(ctx->subcomm);
1088:   return(0);
1089: }

1091: PetscErrorCode NEPDestroy_CISS(NEP nep)
1092: {
1094:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

1097:   NEPCISSResetSubcomm(nep);
1098:   PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1099:   PetscFree(nep->data);
1100:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL);
1101:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL);
1102:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL);
1103:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL);
1104:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL);
1105:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL);
1106:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL);
1107:   return(0);
1108: }

1110: PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1111: {
1113:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;
1114:   PetscBool      isascii;

1117:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1118:   if (isascii) {
1119:     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);
1120:     if (ctx->isreal) {
1121:       PetscViewerASCIIPrintf(viewer,"  exploiting symmetry of integration points\n");
1122:     }
1123:     PetscViewerASCIIPrintf(viewer,"  threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1124:     PetscViewerASCIIPrintf(viewer,"  iterative refinement  { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1125:     if (!ctx->ksp) { NEPCISSGetKSPs(nep,&ctx->num_solve_point,&ctx->ksp); }
1126:     PetscViewerASCIIPushTab(viewer);
1127:     KSPView(ctx->ksp[0],viewer);
1128:     PetscViewerASCIIPopTab(viewer);
1129:   }
1130:   return(0);
1131: }

1133: SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1134: {
1136:   NEP_CISS       *ctx = (NEP_CISS*)nep->data;

1139:   PetscNewLog(nep,&ctx);
1140:   nep->data = ctx;
1141:   /* set default values of parameters */
1142:   ctx->N                  = 32;
1143:   ctx->L                  = 16;
1144:   ctx->M                  = ctx->N/4;
1145:   ctx->delta              = 1e-12;
1146:   ctx->L_max              = 64;
1147:   ctx->spurious_threshold = 1e-4;
1148:   ctx->isreal             = PETSC_FALSE;
1149:   ctx->npart              = 1;

1151:   nep->useds = PETSC_TRUE;

1153:   nep->ops->solve          = NEPSolve_CISS;
1154:   nep->ops->setup          = NEPSetUp_CISS;
1155:   nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1156:   nep->ops->reset          = NEPReset_CISS;
1157:   nep->ops->destroy        = NEPDestroy_CISS;
1158:   nep->ops->view           = NEPView_CISS;

1160:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS);
1161:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS);
1162:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS);
1163:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS);
1164:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS);
1165:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS);
1166:   PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS);
1167:   return(0);
1168: }