Actual source code: nciss.c
slepc-3.18.0 2022-10-01
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: for (i=0;i<contour->npoints;i++) {
122: p_id = i*contour->subcomm->n + contour->subcomm->color;
123: MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&Amat);
124: if (T != P) MatDuplicate(P,MAT_DO_NOT_COPY_VALUES,&Pmat); else Pmat = Amat;
125: if (contour->subcomm->n == 1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) NEPComputeFunction(nep,ctx->omega[p_id],Amat,Pmat);
126: else NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],Amat,Pmat,PETSC_FALSE);
127: NEP_KSPSetOperators(contour->ksp[i],Amat,Pmat);
128: MatDestroy(&Amat);
129: if (T != P) MatDestroy(&Pmat);
130: }
131: return 0;
132: }
134: /*
135: Y_i = F(z_i)^{-1}Fp(z_i)V for every integration point, Y=[Y_i] is in the context
136: */
137: static PetscErrorCode NEPCISSSolve(NEP nep,Mat dT,BV V,PetscInt L_start,PetscInt L_end)
138: {
139: NEP_CISS *ctx = (NEP_CISS*)nep->data;
140: SlepcContourData contour = ctx->contour;
141: PetscInt i,p_id;
142: Mat MV,BMV=NULL,MC;
144: BVSetActiveColumns(V,L_start,L_end);
145: BVGetMat(V,&MV);
146: for (i=0;i<contour->npoints;i++) {
147: p_id = i*contour->subcomm->n + contour->subcomm->color;
148: if (contour->subcomm->n==1 || nep->fui==NEP_USER_INTERFACE_CALLBACK) NEPComputeJacobian(nep,ctx->omega[p_id],dT);
149: else NEPComputeFunctionSubcomm(nep,ctx->omega[p_id],dT,NULL,PETSC_TRUE);
150: BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end);
151: BVGetMat(ctx->Y,&MC);
152: if (!i) {
153: MatProductCreate(dT,MV,NULL,&BMV);
154: MatProductSetType(BMV,MATPRODUCT_AB);
155: MatProductSetFromOptions(BMV);
156: MatProductSymbolic(BMV);
157: }
158: MatProductNumeric(BMV);
159: KSPMatSolve(contour->ksp[i],BMV,MC);
160: BVRestoreMat(ctx->Y,&MC);
161: }
162: MatDestroy(&BMV);
163: BVRestoreMat(V,&MV);
164: return 0;
165: }
167: PetscErrorCode NEPSetUp_CISS(NEP nep)
168: {
169: NEP_CISS *ctx = (NEP_CISS*)nep->data;
170: SlepcContourData contour;
171: PetscInt nwork;
172: PetscBool istrivial,isellipse,flg;
173: NEP_CISS_PROJECT dsctxf;
174: PetscObjectId id;
175: PetscObjectState state;
176: Vec v0;
178: if (nep->ncv==PETSC_DEFAULT) nep->ncv = ctx->L_max*ctx->M;
179: else {
180: ctx->L_max = nep->ncv/ctx->M;
181: if (!ctx->L_max) {
182: ctx->L_max = 1;
183: nep->ncv = ctx->L_max*ctx->M;
184: }
185: }
186: ctx->L = PetscMin(ctx->L,ctx->L_max);
187: if (nep->max_it==PETSC_DEFAULT) nep->max_it = 5;
188: if (nep->mpd==PETSC_DEFAULT) nep->mpd = nep->ncv;
189: if (!nep->which) nep->which = NEP_ALL;
191: NEPCheckUnsupported(nep,NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);
193: /* check region */
194: RGIsTrivial(nep->rg,&istrivial);
196: RGGetComplement(nep->rg,&flg);
198: PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
201: /* if the region has changed, then reset contour data */
202: PetscObjectGetId((PetscObject)nep->rg,&id);
203: PetscObjectStateGet((PetscObject)nep->rg,&state);
204: if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
205: SlepcContourDataDestroy(&ctx->contour);
206: PetscInfo(nep,"Resetting the contour data structure due to a change of region\n");
207: ctx->rgid = id; ctx->rgstate = state;
208: }
210: /* create contour data structure */
211: if (!ctx->contour) {
212: RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj);
213: SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour);
214: }
216: NEPAllocateSolution(nep,0);
217: if (ctx->weight) PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
218: PetscMalloc4(ctx->N,&ctx->weight,ctx->N,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
220: /* allocate basis vectors */
221: BVDestroy(&ctx->S);
222: BVDuplicateResize(nep->V,ctx->L*ctx->M,&ctx->S);
223: BVDestroy(&ctx->V);
224: BVDuplicateResize(nep->V,ctx->L,&ctx->V);
226: contour = ctx->contour;
227: if (contour->subcomm && contour->subcomm->n != 1 && nep->fui==NEP_USER_INTERFACE_CALLBACK) {
228: NEPComputeFunction(nep,0,nep->function,nep->function_pre);
229: SlepcContourRedundantMat(contour,1,&nep->function,(nep->function!=nep->function_pre)?&nep->function_pre:NULL);
230: } else SlepcContourRedundantMat(contour,nep->nt,nep->A,nep->P);
231: if (contour->pA) {
232: if (!ctx->J) MatDuplicate(contour->pA[0],MAT_DO_NOT_COPY_VALUES,&ctx->J);
233: BVGetColumn(ctx->V,0,&v0);
234: SlepcContourScatterCreate(contour,v0);
235: BVRestoreColumn(ctx->V,0,&v0);
236: BVDestroy(&ctx->pV);
237: BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV);
238: BVSetSizesFromVec(ctx->pV,contour->xsub,nep->n);
239: BVSetFromOptions(ctx->pV);
240: BVResize(ctx->pV,ctx->L,PETSC_FALSE);
241: }
243: BVDestroy(&ctx->Y);
244: if (contour->pA) {
245: BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y);
246: BVSetSizesFromVec(ctx->Y,contour->xsub,nep->n);
247: BVSetFromOptions(ctx->Y);
248: BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE);
249: } else BVDuplicateResize(nep->V,contour->npoints*ctx->L,&ctx->Y);
251: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) DSSetType(nep->ds,DSGNHEP);
252: else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) DSSetType(nep->ds,DSNHEP);
253: else {
254: DSSetType(nep->ds,DSNEP);
255: DSSetMethod(nep->ds,1);
256: DSNEPSetRG(nep->ds,nep->rg);
257: if (nep->fui==NEP_USER_INTERFACE_SPLIT) DSNEPSetFN(nep->ds,nep->nt,nep->f);
258: else {
259: PetscNew(&dsctxf);
260: DSNEPSetComputeMatrixFunction(nep->ds,NEPContourDSComputeMatrix,dsctxf);
261: dsctxf->nep = nep;
262: ctx->dsctxf = dsctxf;
263: }
264: }
265: DSAllocate(nep->ds,nep->ncv);
266: nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
267: NEPSetWorkVecs(nep,nwork);
268: return 0;
269: }
271: PetscErrorCode NEPSolve_CISS(NEP nep)
272: {
273: NEP_CISS *ctx = (NEP_CISS*)nep->data;
274: SlepcContourData contour = ctx->contour;
275: Mat X,M,E,T,P,J;
276: BV V;
277: PetscInt i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
278: PetscScalar *Mu,*H0,*H1,*rr,*temp,center;
279: PetscReal error,max_error,radius,rgscale,est_eig,eta;
280: PetscBool isellipse,*fl1;
281: Vec si;
282: SlepcSC sc;
283: PetscRandom rand;
285: DSSetFromOptions(nep->ds);
286: DSGetSlepcSC(nep->ds,&sc);
287: sc->comparison = SlepcCompareLargestMagnitude;
288: sc->comparisonctx = NULL;
289: sc->map = NULL;
290: sc->mapobj = NULL;
291: DSGetLeadingDimension(nep->ds,&ld);
292: RGComputeQuadrature(nep->rg,RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight);
293: if (contour->pA) {
294: T = contour->pA[0];
295: P = ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && contour->pP))? contour->pP[0]: T;
296: } else {
297: T = nep->function;
298: P = nep->function_pre? nep->function_pre: nep->function;
299: }
300: NEPCISSSetUp(nep,T,P);
301: BVSetActiveColumns(ctx->V,0,ctx->L);
302: BVSetRandomSign(ctx->V);
303: BVGetRandomContext(ctx->V,&rand);
304: if (contour->pA) {
305: J = ctx->J;
306: V = ctx->pV;
307: BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
308: } else {
309: J = nep->jacobian;
310: V = ctx->V;
311: }
312: NEPCISSSolve(nep,J,V,0,ctx->L);
313: PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
314: if (isellipse) {
315: BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig);
316: PetscInfo(nep,"Estimated eigenvalue count: %f\n",(double)est_eig);
317: eta = PetscPowReal(10.0,-PetscLog10Real(nep->tol)/ctx->N);
318: L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
319: if (L_add>ctx->L_max-ctx->L) {
320: PetscInfo(nep,"Number of eigenvalues inside the contour path may be too large\n");
321: L_add = ctx->L_max-ctx->L;
322: }
323: }
324: /* Updates L after estimate the number of eigenvalue */
325: if (L_add>0) {
326: PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add);
327: BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
328: BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
329: BVSetRandomSign(ctx->V);
330: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
331: ctx->L += L_add;
332: NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
333: }
335: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
336: for (i=0;i<ctx->refine_blocksize;i++) {
337: BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
338: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
339: PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
340: SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
341: PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
342: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
343: L_add = L_base;
344: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
345: PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add);
346: BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
347: BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
348: BVSetRandomSign(ctx->V);
349: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
350: ctx->L += L_add;
351: NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
352: if (L_add) {
353: PetscFree2(Mu,H0);
354: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
355: }
356: }
358: RGGetScale(nep->rg,&rgscale);
359: RGEllipseGetParameters(nep->rg,¢er,&radius,NULL);
361: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);
363: while (nep->reason == NEP_CONVERGED_ITERATING) {
364: nep->its++;
365: for (inner=0;inner<=ctx->refine_inner;inner++) {
366: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
367: BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
368: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
369: PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
370: SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
371: PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
372: } else {
373: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
374: /* compute SVD of S */
375: 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);
376: }
377: PetscInfo(nep,"Estimated rank: nv = %" PetscInt_FMT "\n",nv);
378: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
379: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
380: BVSetActiveColumns(ctx->S,0,ctx->L);
381: BVSetActiveColumns(ctx->V,0,ctx->L);
382: BVCopy(ctx->S,ctx->V);
383: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
384: NEPCISSSolve(nep,J,V,0,ctx->L);
385: } else break;
386: }
387: nep->nconv = 0;
388: if (nv == 0) { nep->reason = NEP_CONVERGED_TOL; break; }
389: else {
390: /* Extracting eigenpairs */
391: DSSetDimensions(nep->ds,nv,0,0);
392: DSSetState(nep->ds,DS_STATE_RAW);
393: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
394: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
395: CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1);
396: DSGetArray(nep->ds,DS_MAT_A,&temp);
397: for (j=0;j<nv;j++)
398: for (i=0;i<nv;i++)
399: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
400: DSRestoreArray(nep->ds,DS_MAT_A,&temp);
401: DSGetArray(nep->ds,DS_MAT_B,&temp);
402: for (j=0;j<nv;j++)
403: for (i=0;i<nv;i++)
404: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
405: DSRestoreArray(nep->ds,DS_MAT_B,&temp);
406: } else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) {
407: BVSetActiveColumns(ctx->S,0,nv);
408: DSGetArray(nep->ds,DS_MAT_A,&temp);
409: for (i=0;i<nv;i++) PetscArraycpy(temp+i*ld,H0+i*nv,nv);
410: DSRestoreArray(nep->ds,DS_MAT_A,&temp);
411: } else {
412: BVSetActiveColumns(ctx->S,0,nv);
413: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
414: for (i=0;i<nep->nt;i++) {
415: DSGetMat(nep->ds,DSMatExtra[i],&E);
416: BVMatProject(ctx->S,nep->A[i],ctx->S,E);
417: DSRestoreMat(nep->ds,DSMatExtra[i],&E);
418: }
419: } else { ctx->dsctxf->Q = ctx->S; }
420: }
421: DSSolve(nep->ds,nep->eigr,nep->eigi);
422: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
423: DSGetDimensions(nep->ds,NULL,NULL,NULL,&nv);
424: if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
425: for (i=0;i<nv;i++) {
426: nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
427: }
428: }
429: PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
430: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
431: DSGetMat(nep->ds,DS_MAT_X,&X);
432: SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1);
433: DSRestoreMat(nep->ds,DS_MAT_X,&X);
434: RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside);
435: for (i=0;i<nv;i++) {
436: if (fl1[i] && inside[i]>=0) {
437: rr[i] = 1.0;
438: nep->nconv++;
439: } else rr[i] = 0.0;
440: }
441: DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv);
442: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
443: if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
444: for (i=0;i<nv;i++) nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
445: }
446: PetscFree3(fl1,inside,rr);
447: BVSetActiveColumns(nep->V,0,nv);
448: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
449: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
450: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
451: BVSetActiveColumns(ctx->S,0,nv);
452: BVCopy(ctx->S,nep->V);
453: DSGetMat(nep->ds,DS_MAT_X,&X);
454: BVMultInPlace(ctx->S,X,0,nep->nconv);
455: BVMultInPlace(nep->V,X,0,nep->nconv);
456: DSRestoreMat(nep->ds,DS_MAT_X,&X);
457: } else {
458: DSGetMat(nep->ds,DS_MAT_X,&X);
459: BVMultInPlace(ctx->S,X,0,nep->nconv);
460: DSRestoreMat(nep->ds,DS_MAT_X,&X);
461: BVSetActiveColumns(ctx->S,0,nep->nconv);
462: BVCopy(ctx->S,nep->V);
463: }
464: max_error = 0.0;
465: for (i=0;i<nep->nconv;i++) {
466: BVGetColumn(nep->V,i,&si);
467: VecNormalize(si,NULL);
468: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error);
469: (*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx);
470: BVRestoreColumn(nep->V,i,&si);
471: max_error = PetscMax(max_error,error);
472: }
473: if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
474: else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
475: else {
476: if (nep->nconv > ctx->L) nv = nep->nconv;
477: else if (ctx->L > nv) nv = ctx->L;
478: nv = PetscMin(nv,ctx->L*ctx->M);
479: MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
480: MatSetRandom(M,rand);
481: BVSetActiveColumns(ctx->S,0,nv);
482: BVMultInPlace(ctx->S,M,0,ctx->L);
483: MatDestroy(&M);
484: BVSetActiveColumns(ctx->S,0,ctx->L);
485: BVSetActiveColumns(ctx->V,0,ctx->L);
486: BVCopy(ctx->S,ctx->V);
487: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
488: NEPCISSSolve(nep,J,V,0,ctx->L);
489: }
490: }
491: }
492: PetscFree2(Mu,H0);
493: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscFree(H1);
494: return 0;
495: }
497: static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
498: {
499: NEP_CISS *ctx = (NEP_CISS*)nep->data;
500: PetscInt oN,oL,oM,oLmax,onpart;
501: PetscMPIInt size;
503: oN = ctx->N;
504: if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
505: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
506: } else {
509: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
510: }
511: oL = ctx->L;
512: if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
513: ctx->L = 16;
514: } else {
516: ctx->L = bs;
517: }
518: oM = ctx->M;
519: if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
520: ctx->M = ctx->N/4;
521: } else {
524: ctx->M = PetscMax(ms,2);
525: }
526: onpart = ctx->npart;
527: if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
528: ctx->npart = 1;
529: } else {
530: MPI_Comm_size(PetscObjectComm((PetscObject)nep),&size);
532: ctx->npart = npart;
533: }
534: oLmax = ctx->L_max;
535: if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
536: ctx->L_max = 64;
537: } else {
539: ctx->L_max = PetscMax(bsmax,ctx->L);
540: }
541: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
542: SlepcContourDataDestroy(&ctx->contour);
543: PetscInfo(nep,"Resetting the contour data structure due to a change of parameters\n");
544: nep->state = NEP_STATE_INITIAL;
545: }
546: ctx->isreal = realmats;
547: if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) nep->state = NEP_STATE_INITIAL;
548: return 0;
549: }
551: /*@
552: NEPCISSSetSizes - Sets the values of various size parameters in the CISS solver.
554: Logically Collective on nep
556: Input Parameters:
557: + nep - the nonlinear eigensolver context
558: . ip - number of integration points
559: . bs - block size
560: . ms - moment size
561: . npart - number of partitions when splitting the communicator
562: . bsmax - max block size
563: - realmats - T(z) is real for real z
565: Options Database Keys:
566: + -nep_ciss_integration_points - Sets the number of integration points
567: . -nep_ciss_blocksize - Sets the block size
568: . -nep_ciss_moments - Sets the moment size
569: . -nep_ciss_partitions - Sets the number of partitions
570: . -nep_ciss_maxblocksize - Sets the maximum block size
571: - -nep_ciss_realmats - T(z) is real for real z
573: Notes:
574: The default number of partitions is 1. This means the internal KSP object is shared
575: among all processes of the NEP communicator. Otherwise, the communicator is split
576: into npart communicators, so that npart KSP solves proceed simultaneously.
578: The realmats flag can be set to true when T(.) is guaranteed to be real
579: when the argument is a real value, for example, when all matrices in
580: the split form are real. When set to true, the solver avoids some computations.
582: Level: advanced
584: .seealso: NEPCISSGetSizes()
585: @*/
586: PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
587: {
595: PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
596: return 0;
597: }
599: static PetscErrorCode NEPCISSGetSizes_CISS(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
600: {
601: NEP_CISS *ctx = (NEP_CISS*)nep->data;
603: if (ip) *ip = ctx->N;
604: if (bs) *bs = ctx->L;
605: if (ms) *ms = ctx->M;
606: if (npart) *npart = ctx->npart;
607: if (bsmax) *bsmax = ctx->L_max;
608: if (realmats) *realmats = ctx->isreal;
609: return 0;
610: }
612: /*@
613: NEPCISSGetSizes - Gets the values of various size parameters in the CISS solver.
615: Not Collective
617: Input Parameter:
618: . nep - the nonlinear eigensolver context
620: Output Parameters:
621: + ip - number of integration points
622: . bs - block size
623: . ms - moment size
624: . npart - number of partitions when splitting the communicator
625: . bsmax - max block size
626: - realmats - T(z) is real for real z
628: Level: advanced
630: .seealso: NEPCISSSetSizes()
631: @*/
632: PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
633: {
635: PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
636: return 0;
637: }
639: static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
640: {
641: NEP_CISS *ctx = (NEP_CISS*)nep->data;
643: if (delta == PETSC_DEFAULT) {
644: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
645: } else {
647: ctx->delta = delta;
648: }
649: if (spur == PETSC_DEFAULT) {
650: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
651: } else {
653: ctx->spurious_threshold = spur;
654: }
655: return 0;
656: }
658: /*@
659: NEPCISSSetThreshold - Sets the values of various threshold parameters in
660: the CISS solver.
662: Logically Collective on nep
664: Input Parameters:
665: + nep - the nonlinear eigensolver context
666: . delta - threshold for numerical rank
667: - spur - spurious threshold (to discard spurious eigenpairs)
669: Options Database Keys:
670: + -nep_ciss_delta - Sets the delta
671: - -nep_ciss_spurious_threshold - Sets the spurious threshold
673: Level: advanced
675: .seealso: NEPCISSGetThreshold()
676: @*/
677: PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
678: {
682: PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
683: return 0;
684: }
686: static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
687: {
688: NEP_CISS *ctx = (NEP_CISS*)nep->data;
690: if (delta) *delta = ctx->delta;
691: if (spur) *spur = ctx->spurious_threshold;
692: return 0;
693: }
695: /*@
696: NEPCISSGetThreshold - Gets the values of various threshold parameters in
697: the CISS solver.
699: Not Collective
701: Input Parameter:
702: . nep - the nonlinear eigensolver context
704: Output Parameters:
705: + delta - threshold for numerical rank
706: - spur - spurious threshold (to discard spurious eigenpairs)
708: Level: advanced
710: .seealso: NEPCISSSetThreshold()
711: @*/
712: PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
713: {
715: PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
716: return 0;
717: }
719: static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
720: {
721: NEP_CISS *ctx = (NEP_CISS*)nep->data;
723: if (inner == PETSC_DEFAULT) {
724: ctx->refine_inner = 0;
725: } else {
727: ctx->refine_inner = inner;
728: }
729: if (blsize == PETSC_DEFAULT) {
730: ctx->refine_blocksize = 0;
731: } else {
733: ctx->refine_blocksize = blsize;
734: }
735: return 0;
736: }
738: /*@
739: NEPCISSSetRefinement - Sets the values of various refinement parameters
740: in the CISS solver.
742: Logically Collective on nep
744: Input Parameters:
745: + nep - the nonlinear eigensolver context
746: . inner - number of iterative refinement iterations (inner loop)
747: - blsize - number of iterative refinement iterations (blocksize loop)
749: Options Database Keys:
750: + -nep_ciss_refine_inner - Sets number of inner iterations
751: - -nep_ciss_refine_blocksize - Sets number of blocksize iterations
753: Level: advanced
755: .seealso: NEPCISSGetRefinement()
756: @*/
757: PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
758: {
762: PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
763: return 0;
764: }
766: static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
767: {
768: NEP_CISS *ctx = (NEP_CISS*)nep->data;
770: if (inner) *inner = ctx->refine_inner;
771: if (blsize) *blsize = ctx->refine_blocksize;
772: return 0;
773: }
775: /*@
776: NEPCISSGetRefinement - Gets the values of various refinement parameters
777: in the CISS solver.
779: Not Collective
781: Input Parameter:
782: . nep - the nonlinear eigensolver context
784: Output Parameters:
785: + inner - number of iterative refinement iterations (inner loop)
786: - blsize - number of iterative refinement iterations (blocksize loop)
788: Level: advanced
790: .seealso: NEPCISSSetRefinement()
791: @*/
792: PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
793: {
795: PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
796: return 0;
797: }
799: static PetscErrorCode NEPCISSSetExtraction_CISS(NEP nep,NEPCISSExtraction extraction)
800: {
801: NEP_CISS *ctx = (NEP_CISS*)nep->data;
803: if (ctx->extraction != extraction) {
804: ctx->extraction = extraction;
805: nep->state = NEP_STATE_INITIAL;
806: }
807: return 0;
808: }
810: /*@
811: NEPCISSSetExtraction - Sets the extraction technique used in the CISS solver.
813: Logically Collective on nep
815: Input Parameters:
816: + nep - the nonlinear eigensolver context
817: - extraction - the extraction technique
819: Options Database Key:
820: . -nep_ciss_extraction - Sets the extraction technique (either 'ritz', 'hankel' or 'caa')
822: Notes:
823: By default, the Rayleigh-Ritz extraction is used (NEP_CISS_EXTRACTION_RITZ).
825: If the 'hankel' or the 'caa' option is specified (NEP_CISS_EXTRACTION_HANKEL or
826: NEP_CISS_EXTRACTION_CAA), then the Block Hankel method, or the Communication-avoiding
827: Arnoldi method, respectively, is used for extracting eigenpairs.
829: Level: advanced
831: .seealso: NEPCISSGetExtraction(), NEPCISSExtraction
832: @*/
833: PetscErrorCode NEPCISSSetExtraction(NEP nep,NEPCISSExtraction extraction)
834: {
837: PetscTryMethod(nep,"NEPCISSSetExtraction_C",(NEP,NEPCISSExtraction),(nep,extraction));
838: return 0;
839: }
841: static PetscErrorCode NEPCISSGetExtraction_CISS(NEP nep,NEPCISSExtraction *extraction)
842: {
843: NEP_CISS *ctx = (NEP_CISS*)nep->data;
845: *extraction = ctx->extraction;
846: return 0;
847: }
849: /*@
850: NEPCISSGetExtraction - Gets the extraction technique used in the CISS solver.
852: Not Collective
854: Input Parameter:
855: . nep - the nonlinear eigensolver context
857: Output Parameters:
858: . extraction - extraction technique
860: Level: advanced
862: .seealso: NEPCISSSetExtraction() NEPCISSExtraction
863: @*/
864: PetscErrorCode NEPCISSGetExtraction(NEP nep,NEPCISSExtraction *extraction)
865: {
868: PetscUseMethod(nep,"NEPCISSGetExtraction_C",(NEP,NEPCISSExtraction*),(nep,extraction));
869: return 0;
870: }
872: static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
873: {
874: NEP_CISS *ctx = (NEP_CISS*)nep->data;
875: SlepcContourData contour;
876: PetscInt i;
877: PC pc;
878: MPI_Comm child;
880: if (!ctx->contour) { /* initialize contour data structure first */
881: RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj);
882: SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour);
883: }
884: contour = ctx->contour;
885: if (!contour->ksp) {
886: PetscMalloc1(contour->npoints,&contour->ksp);
887: PetscSubcommGetChild(contour->subcomm,&child);
888: for (i=0;i<contour->npoints;i++) {
889: KSPCreate(child,&contour->ksp[i]);
890: PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)nep,1);
891: KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)nep)->prefix);
892: KSPAppendOptionsPrefix(contour->ksp[i],"nep_ciss_");
893: PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)nep)->options);
894: KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE);
895: KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
896: KSPGetPC(contour->ksp[i],&pc);
897: if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
898: KSPSetType(contour->ksp[i],KSPBCGS);
899: PCSetType(pc,PCBJACOBI);
900: } else {
901: KSPSetType(contour->ksp[i],KSPPREONLY);
902: PCSetType(pc,PCLU);
903: }
904: }
905: }
906: if (nsolve) *nsolve = contour->npoints;
907: if (ksp) *ksp = contour->ksp;
908: return 0;
909: }
911: /*@C
912: NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
913: the CISS solver.
915: Not Collective
917: Input Parameter:
918: . nep - nonlinear eigenvalue solver
920: Output Parameters:
921: + nsolve - number of solver objects
922: - ksp - array of linear solver object
924: Notes:
925: The number of KSP solvers is equal to the number of integration points divided by
926: the number of partitions. This value is halved in the case of real matrices with
927: a region centered at the real axis.
929: Level: advanced
931: .seealso: NEPCISSSetSizes()
932: @*/
933: PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
934: {
936: PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
937: return 0;
938: }
940: PetscErrorCode NEPReset_CISS(NEP nep)
941: {
942: NEP_CISS *ctx = (NEP_CISS*)nep->data;
944: BVDestroy(&ctx->S);
945: BVDestroy(&ctx->V);
946: BVDestroy(&ctx->Y);
947: SlepcContourDataReset(ctx->contour);
948: MatDestroy(&ctx->J);
949: BVDestroy(&ctx->pV);
950: if (ctx->extraction == NEP_CISS_EXTRACTION_RITZ && nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscFree(ctx->dsctxf);
951: return 0;
952: }
954: PetscErrorCode NEPSetFromOptions_CISS(NEP nep,PetscOptionItems *PetscOptionsObject)
955: {
956: NEP_CISS *ctx = (NEP_CISS*)nep->data;
957: PetscReal r1,r2;
958: PetscInt i,i1,i2,i3,i4,i5,i6,i7;
959: PetscBool b1,flg,flg2,flg3,flg4,flg5,flg6;
960: NEPCISSExtraction extraction;
962: PetscOptionsHeadBegin(PetscOptionsObject,"NEP CISS Options");
964: NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1);
965: PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,&flg);
966: PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,&flg2);
967: PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,&flg3);
968: PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,&flg4);
969: PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,&flg5);
970: PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,&flg6);
971: if (flg || flg2 || flg3 || flg4 || flg5 || flg6) NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1);
973: NEPCISSGetThreshold(nep,&r1,&r2);
974: PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,&flg);
975: PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,&flg2);
976: if (flg || flg2) NEPCISSSetThreshold(nep,r1,r2);
978: NEPCISSGetRefinement(nep,&i6,&i7);
979: PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,&flg);
980: PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,&flg2);
981: if (flg || flg2) NEPCISSSetRefinement(nep,i6,i7);
983: PetscOptionsEnum("-nep_ciss_extraction","Extraction technique","NEPCISSSetExtraction",NEPCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
984: if (flg) NEPCISSSetExtraction(nep,extraction);
986: PetscOptionsHeadEnd();
988: if (!nep->rg) NEPGetRG(nep,&nep->rg);
989: RGSetFromOptions(nep->rg); /* this is necessary here to set useconj */
990: if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
991: for (i=0;i<ctx->contour->npoints;i++) KSPSetFromOptions(ctx->contour->ksp[i]);
992: PetscSubcommSetFromOptions(ctx->contour->subcomm);
993: return 0;
994: }
996: PetscErrorCode NEPDestroy_CISS(NEP nep)
997: {
998: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1000: SlepcContourDataDestroy(&ctx->contour);
1001: PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1002: PetscFree(nep->data);
1003: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL);
1004: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL);
1005: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL);
1006: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL);
1007: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL);
1008: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL);
1009: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NULL);
1010: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NULL);
1011: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL);
1012: return 0;
1013: }
1015: PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1016: {
1017: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1018: PetscBool isascii;
1019: PetscViewer sviewer;
1021: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1022: if (isascii) {
1023: 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);
1024: if (ctx->isreal) PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n");
1025: PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1026: PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize);
1027: PetscViewerASCIIPrintf(viewer," extraction: %s\n",NEPCISSExtractions[ctx->extraction]);
1028: if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
1029: PetscViewerASCIIPushTab(viewer);
1030: if (ctx->npart>1 && ctx->contour->subcomm) {
1031: PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1032: if (!ctx->contour->subcomm->color) KSPView(ctx->contour->ksp[0],sviewer);
1033: PetscViewerFlush(sviewer);
1034: PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1035: PetscViewerFlush(viewer);
1036: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1037: PetscViewerASCIIPopSynchronized(viewer);
1038: } else KSPView(ctx->contour->ksp[0],viewer);
1039: PetscViewerASCIIPopTab(viewer);
1040: }
1041: return 0;
1042: }
1044: SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1045: {
1046: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1048: PetscNew(&ctx);
1049: nep->data = ctx;
1050: /* set default values of parameters */
1051: ctx->N = 32;
1052: ctx->L = 16;
1053: ctx->M = ctx->N/4;
1054: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
1055: ctx->L_max = 64;
1056: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1057: ctx->isreal = PETSC_FALSE;
1058: ctx->npart = 1;
1060: nep->useds = PETSC_TRUE;
1062: nep->ops->solve = NEPSolve_CISS;
1063: nep->ops->setup = NEPSetUp_CISS;
1064: nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1065: nep->ops->reset = NEPReset_CISS;
1066: nep->ops->destroy = NEPDestroy_CISS;
1067: nep->ops->view = NEPView_CISS;
1069: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS);
1070: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS);
1071: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS);
1072: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS);
1073: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS);
1074: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS);
1075: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NEPCISSSetExtraction_CISS);
1076: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NEPCISSGetExtraction_CISS);
1077: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS);
1078: return 0;
1079: }