Actual source code: nciss.c
slepc-3.11.2 2019-07-30
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,¢er,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,¢er,&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,¢er,&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: }