Actual source code: nciss.c
slepc-3.17.0 2022-03-31
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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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);
219: PetscLogObjectMemory((PetscObject)nep,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));
221: /* allocate basis vectors */
222: BVDestroy(&ctx->S);
223: BVDuplicateResize(nep->V,ctx->L*ctx->M,&ctx->S);
224: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->S);
225: BVDestroy(&ctx->V);
226: BVDuplicateResize(nep->V,ctx->L,&ctx->V);
227: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->V);
229: contour = ctx->contour;
230: if (contour->subcomm && contour->subcomm->n != 1 && nep->fui==NEP_USER_INTERFACE_CALLBACK) {
231: NEPComputeFunction(nep,0,nep->function,nep->function_pre);
232: SlepcContourRedundantMat(contour,1,&nep->function,(nep->function!=nep->function_pre)?&nep->function_pre:NULL);
233: } else SlepcContourRedundantMat(contour,nep->nt,nep->A,nep->P);
234: if (contour->pA) {
235: if (!ctx->J) {
236: MatDuplicate(contour->pA[0],MAT_DO_NOT_COPY_VALUES,&ctx->J);
237: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->J);
238: }
239: BVGetColumn(ctx->V,0,&v0);
240: SlepcContourScatterCreate(contour,v0);
241: BVRestoreColumn(ctx->V,0,&v0);
242: BVDestroy(&ctx->pV);
243: BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV);
244: BVSetSizesFromVec(ctx->pV,contour->xsub,nep->n);
245: BVSetFromOptions(ctx->pV);
246: BVResize(ctx->pV,ctx->L,PETSC_FALSE);
247: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->pV);
248: }
250: BVDestroy(&ctx->Y);
251: if (contour->pA) {
252: BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y);
253: BVSetSizesFromVec(ctx->Y,contour->xsub,nep->n);
254: BVSetFromOptions(ctx->Y);
255: BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE);
256: } else BVDuplicateResize(nep->V,contour->npoints*ctx->L,&ctx->Y);
258: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) DSSetType(nep->ds,DSGNHEP);
259: else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) DSSetType(nep->ds,DSNHEP);
260: else {
261: DSSetType(nep->ds,DSNEP);
262: DSSetMethod(nep->ds,1);
263: DSNEPSetRG(nep->ds,nep->rg);
264: if (nep->fui==NEP_USER_INTERFACE_SPLIT) DSNEPSetFN(nep->ds,nep->nt,nep->f);
265: else {
266: PetscNew(&dsctxf);
267: DSNEPSetComputeMatrixFunction(nep->ds,NEPContourDSComputeMatrix,dsctxf);
268: dsctxf->nep = nep;
269: ctx->dsctxf = dsctxf;
270: }
271: }
272: DSAllocate(nep->ds,nep->ncv);
273: nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
274: NEPSetWorkVecs(nep,nwork);
275: PetscFunctionReturn(0);
276: }
278: PetscErrorCode NEPSolve_CISS(NEP nep)
279: {
280: NEP_CISS *ctx = (NEP_CISS*)nep->data;
281: SlepcContourData contour = ctx->contour;
282: Mat X,M,E,T,P,J;
283: BV V;
284: PetscInt i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
285: PetscScalar *Mu,*H0,*H1,*rr,*temp,center;
286: PetscReal error,max_error,radius,rgscale,est_eig,eta;
287: PetscBool isellipse,*fl1;
288: Vec si;
289: SlepcSC sc;
290: PetscRandom rand;
292: DSSetFromOptions(nep->ds);
293: DSGetSlepcSC(nep->ds,&sc);
294: sc->comparison = SlepcCompareLargestMagnitude;
295: sc->comparisonctx = NULL;
296: sc->map = NULL;
297: sc->mapobj = NULL;
298: DSGetLeadingDimension(nep->ds,&ld);
299: RGComputeQuadrature(nep->rg,RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight);
300: if (contour->pA) {
301: T = contour->pA[0];
302: P = ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && contour->pP))? contour->pP[0]: T;
303: } else {
304: T = nep->function;
305: P = nep->function_pre? nep->function_pre: nep->function;
306: }
307: NEPCISSSetUp(nep,T,P);
308: BVSetActiveColumns(ctx->V,0,ctx->L);
309: BVSetRandomSign(ctx->V);
310: BVGetRandomContext(ctx->V,&rand);
311: if (contour->pA) {
312: J = ctx->J;
313: V = ctx->pV;
314: BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
315: } else {
316: J = nep->jacobian;
317: V = ctx->V;
318: }
319: NEPCISSSolve(nep,J,V,0,ctx->L);
320: PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
321: if (isellipse) {
322: BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig);
323: PetscInfo(nep,"Estimated eigenvalue count: %f\n",(double)est_eig);
324: eta = PetscPowReal(10.0,-PetscLog10Real(nep->tol)/ctx->N);
325: L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
326: if (L_add>ctx->L_max-ctx->L) {
327: PetscInfo(nep,"Number of eigenvalues inside the contour path may be too large\n");
328: L_add = ctx->L_max-ctx->L;
329: }
330: }
331: /* Updates L after estimate the number of eigenvalue */
332: if (L_add>0) {
333: PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add);
334: BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
335: BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
336: BVSetRandomSign(ctx->V);
337: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
338: ctx->L += L_add;
339: NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
340: }
342: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
343: for (i=0;i<ctx->refine_blocksize;i++) {
344: BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
345: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
346: PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
347: SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
348: PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
349: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
350: L_add = L_base;
351: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
352: PetscInfo(nep,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add);
353: BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints);
354: BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add);
355: BVSetRandomSign(ctx->V);
356: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
357: ctx->L += L_add;
358: NEPCISSSolve(nep,J,V,ctx->L-L_add,ctx->L);
359: if (L_add) {
360: PetscFree2(Mu,H0);
361: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
362: }
363: }
365: RGGetScale(nep->rg,&rgscale);
366: RGEllipseGetParameters(nep->rg,¢er,&radius,NULL);
368: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);
370: while (nep->reason == NEP_CONVERGED_ITERATING) {
371: nep->its++;
372: for (inner=0;inner<=ctx->refine_inner;inner++) {
373: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
374: BVDotQuadrature(ctx->Y,V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj);
375: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
376: PetscLogEventBegin(NEP_CISS_SVD,nep,0,0,0);
377: SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
378: PetscLogEventEnd(NEP_CISS_SVD,nep,0,0,0);
379: } else {
380: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
381: /* compute SVD of S */
382: 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);
383: }
384: PetscInfo(nep,"Estimated rank: nv = %" PetscInt_FMT "\n",nv);
385: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
386: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
387: BVSetActiveColumns(ctx->S,0,ctx->L);
388: BVSetActiveColumns(ctx->V,0,ctx->L);
389: BVCopy(ctx->S,ctx->V);
390: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
391: NEPCISSSolve(nep,J,V,0,ctx->L);
392: } else break;
393: }
394: nep->nconv = 0;
395: if (nv == 0) { nep->reason = NEP_CONVERGED_TOL; break; }
396: else {
397: /* Extracting eigenpairs */
398: DSSetDimensions(nep->ds,nv,0,0);
399: DSSetState(nep->ds,DS_STATE_RAW);
400: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
401: CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0);
402: CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1);
403: DSGetArray(nep->ds,DS_MAT_A,&temp);
404: for (j=0;j<nv;j++)
405: for (i=0;i<nv;i++)
406: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
407: DSRestoreArray(nep->ds,DS_MAT_A,&temp);
408: DSGetArray(nep->ds,DS_MAT_B,&temp);
409: for (j=0;j<nv;j++)
410: for (i=0;i<nv;i++)
411: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
412: DSRestoreArray(nep->ds,DS_MAT_B,&temp);
413: } else if (ctx->extraction == NEP_CISS_EXTRACTION_CAA) {
414: BVSetActiveColumns(ctx->S,0,nv);
415: DSGetArray(nep->ds,DS_MAT_A,&temp);
416: for (i=0;i<nv;i++) PetscArraycpy(temp+i*ld,H0+i*nv,nv);
417: DSRestoreArray(nep->ds,DS_MAT_A,&temp);
418: } else {
419: BVSetActiveColumns(ctx->S,0,nv);
420: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
421: for (i=0;i<nep->nt;i++) {
422: DSGetMat(nep->ds,DSMatExtra[i],&E);
423: BVMatProject(ctx->S,nep->A[i],ctx->S,E);
424: DSRestoreMat(nep->ds,DSMatExtra[i],&E);
425: }
426: } else { ctx->dsctxf->Q = ctx->S; }
427: }
428: DSSolve(nep->ds,nep->eigr,nep->eigi);
429: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
430: DSGetDimensions(nep->ds,NULL,NULL,NULL,&nv);
431: if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
432: for (i=0;i<nv;i++) {
433: nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
434: }
435: }
436: PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
437: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
438: DSGetMat(nep->ds,DS_MAT_X,&X);
439: SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1);
440: MatDestroy(&X);
441: RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside);
442: for (i=0;i<nv;i++) {
443: if (fl1[i] && inside[i]>=0) {
444: rr[i] = 1.0;
445: nep->nconv++;
446: } else rr[i] = 0.0;
447: }
448: DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv);
449: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
450: if (ctx->extraction == NEP_CISS_EXTRACTION_CAA || ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
451: for (i=0;i<nv;i++) nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
452: }
453: PetscFree3(fl1,inside,rr);
454: BVSetActiveColumns(nep->V,0,nv);
455: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
456: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) {
457: BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj);
458: BVSetActiveColumns(ctx->S,0,nv);
459: BVCopy(ctx->S,nep->V);
460: DSGetMat(nep->ds,DS_MAT_X,&X);
461: BVMultInPlace(ctx->S,X,0,nep->nconv);
462: BVMultInPlace(nep->V,X,0,nep->nconv);
463: MatDestroy(&X);
464: } else {
465: DSGetMat(nep->ds,DS_MAT_X,&X);
466: BVMultInPlace(ctx->S,X,0,nep->nconv);
467: MatDestroy(&X);
468: BVSetActiveColumns(ctx->S,0,nep->nconv);
469: BVCopy(ctx->S,nep->V);
470: }
471: max_error = 0.0;
472: for (i=0;i<nep->nconv;i++) {
473: BVGetColumn(nep->V,i,&si);
474: VecNormalize(si,NULL);
475: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error);
476: (*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx);
477: BVRestoreColumn(nep->V,i,&si);
478: max_error = PetscMax(max_error,error);
479: }
480: if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
481: else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
482: else {
483: if (nep->nconv > ctx->L) nv = nep->nconv;
484: else if (ctx->L > nv) nv = ctx->L;
485: nv = PetscMin(nv,ctx->L*ctx->M);
486: MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
487: MatSetRandom(M,rand);
488: BVSetActiveColumns(ctx->S,0,nv);
489: BVMultInPlace(ctx->S,M,0,ctx->L);
490: MatDestroy(&M);
491: BVSetActiveColumns(ctx->S,0,ctx->L);
492: BVSetActiveColumns(ctx->V,0,ctx->L);
493: BVCopy(ctx->S,ctx->V);
494: if (contour->pA) BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup);
495: NEPCISSSolve(nep,J,V,0,ctx->L);
496: }
497: }
498: }
499: PetscFree2(Mu,H0);
500: if (ctx->extraction == NEP_CISS_EXTRACTION_HANKEL) PetscFree(H1);
501: PetscFunctionReturn(0);
502: }
504: static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
505: {
506: NEP_CISS *ctx = (NEP_CISS*)nep->data;
507: PetscInt oN,oL,oM,oLmax,onpart;
508: PetscMPIInt size;
510: oN = ctx->N;
511: if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
512: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
513: } else {
516: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
517: }
518: oL = ctx->L;
519: if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
520: ctx->L = 16;
521: } else {
523: ctx->L = bs;
524: }
525: oM = ctx->M;
526: if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
527: ctx->M = ctx->N/4;
528: } else {
531: ctx->M = PetscMax(ms,2);
532: }
533: onpart = ctx->npart;
534: if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
535: ctx->npart = 1;
536: } else {
537: MPI_Comm_size(PetscObjectComm((PetscObject)nep),&size);
539: ctx->npart = npart;
540: }
541: oLmax = ctx->L_max;
542: if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
543: ctx->L_max = 64;
544: } else {
546: ctx->L_max = PetscMax(bsmax,ctx->L);
547: }
548: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
549: SlepcContourDataDestroy(&ctx->contour);
550: PetscInfo(nep,"Resetting the contour data structure due to a change of parameters\n");
551: nep->state = NEP_STATE_INITIAL;
552: }
553: ctx->isreal = realmats;
554: if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) nep->state = NEP_STATE_INITIAL;
555: PetscFunctionReturn(0);
556: }
558: /*@
559: NEPCISSSetSizes - Sets the values of various size parameters in the CISS solver.
561: Logically Collective on nep
563: Input Parameters:
564: + nep - the nonlinear eigensolver context
565: . ip - number of integration points
566: . bs - block size
567: . ms - moment size
568: . npart - number of partitions when splitting the communicator
569: . bsmax - max block size
570: - realmats - T(z) is real for real z
572: Options Database Keys:
573: + -nep_ciss_integration_points - Sets the number of integration points
574: . -nep_ciss_blocksize - Sets the block size
575: . -nep_ciss_moments - Sets the moment size
576: . -nep_ciss_partitions - Sets the number of partitions
577: . -nep_ciss_maxblocksize - Sets the maximum block size
578: - -nep_ciss_realmats - T(z) is real for real z
580: Notes:
581: The default number of partitions is 1. This means the internal KSP object is shared
582: among all processes of the NEP communicator. Otherwise, the communicator is split
583: into npart communicators, so that npart KSP solves proceed simultaneously.
585: The realmats flag can be set to true when T(.) is guaranteed to be real
586: when the argument is a real value, for example, when all matrices in
587: the split form are real. When set to true, the solver avoids some computations.
589: Level: advanced
591: .seealso: NEPCISSGetSizes()
592: @*/
593: PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
594: {
602: PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
603: PetscFunctionReturn(0);
604: }
606: static PetscErrorCode NEPCISSGetSizes_CISS(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
607: {
608: NEP_CISS *ctx = (NEP_CISS*)nep->data;
610: if (ip) *ip = ctx->N;
611: if (bs) *bs = ctx->L;
612: if (ms) *ms = ctx->M;
613: if (npart) *npart = ctx->npart;
614: if (bsmax) *bsmax = ctx->L_max;
615: if (realmats) *realmats = ctx->isreal;
616: PetscFunctionReturn(0);
617: }
619: /*@
620: NEPCISSGetSizes - Gets the values of various size parameters in the CISS solver.
622: Not Collective
624: Input Parameter:
625: . nep - the nonlinear eigensolver context
627: Output Parameters:
628: + ip - number of integration points
629: . bs - block size
630: . ms - moment size
631: . npart - number of partitions when splitting the communicator
632: . bsmax - max block size
633: - realmats - T(z) is real for real z
635: Level: advanced
637: .seealso: NEPCISSSetSizes()
638: @*/
639: PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
640: {
642: PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
643: PetscFunctionReturn(0);
644: }
646: static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
647: {
648: NEP_CISS *ctx = (NEP_CISS*)nep->data;
650: if (delta == PETSC_DEFAULT) {
651: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
652: } else {
654: ctx->delta = delta;
655: }
656: if (spur == PETSC_DEFAULT) {
657: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
658: } else {
660: ctx->spurious_threshold = spur;
661: }
662: PetscFunctionReturn(0);
663: }
665: /*@
666: NEPCISSSetThreshold - Sets the values of various threshold parameters in
667: the CISS solver.
669: Logically Collective on nep
671: Input Parameters:
672: + nep - the nonlinear eigensolver context
673: . delta - threshold for numerical rank
674: - spur - spurious threshold (to discard spurious eigenpairs)
676: Options Database Keys:
677: + -nep_ciss_delta - Sets the delta
678: - -nep_ciss_spurious_threshold - Sets the spurious threshold
680: Level: advanced
682: .seealso: NEPCISSGetThreshold()
683: @*/
684: PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
685: {
689: PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
690: PetscFunctionReturn(0);
691: }
693: static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
694: {
695: NEP_CISS *ctx = (NEP_CISS*)nep->data;
697: if (delta) *delta = ctx->delta;
698: if (spur) *spur = ctx->spurious_threshold;
699: PetscFunctionReturn(0);
700: }
702: /*@
703: NEPCISSGetThreshold - Gets the values of various threshold parameters in
704: the CISS solver.
706: Not Collective
708: Input Parameter:
709: . nep - the nonlinear eigensolver context
711: Output Parameters:
712: + delta - threshold for numerical rank
713: - spur - spurious threshold (to discard spurious eigenpairs)
715: Level: advanced
717: .seealso: NEPCISSSetThreshold()
718: @*/
719: PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
720: {
722: PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
723: PetscFunctionReturn(0);
724: }
726: static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
727: {
728: NEP_CISS *ctx = (NEP_CISS*)nep->data;
730: if (inner == PETSC_DEFAULT) {
731: ctx->refine_inner = 0;
732: } else {
734: ctx->refine_inner = inner;
735: }
736: if (blsize == PETSC_DEFAULT) {
737: ctx->refine_blocksize = 0;
738: } else {
740: ctx->refine_blocksize = blsize;
741: }
742: PetscFunctionReturn(0);
743: }
745: /*@
746: NEPCISSSetRefinement - Sets the values of various refinement parameters
747: in the CISS solver.
749: Logically Collective on nep
751: Input Parameters:
752: + nep - the nonlinear eigensolver context
753: . inner - number of iterative refinement iterations (inner loop)
754: - blsize - number of iterative refinement iterations (blocksize loop)
756: Options Database Keys:
757: + -nep_ciss_refine_inner - Sets number of inner iterations
758: - -nep_ciss_refine_blocksize - Sets number of blocksize iterations
760: Level: advanced
762: .seealso: NEPCISSGetRefinement()
763: @*/
764: PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
765: {
769: PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
770: PetscFunctionReturn(0);
771: }
773: static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
774: {
775: NEP_CISS *ctx = (NEP_CISS*)nep->data;
777: if (inner) *inner = ctx->refine_inner;
778: if (blsize) *blsize = ctx->refine_blocksize;
779: PetscFunctionReturn(0);
780: }
782: /*@
783: NEPCISSGetRefinement - Gets the values of various refinement parameters
784: in the CISS solver.
786: Not Collective
788: Input Parameter:
789: . nep - the nonlinear eigensolver context
791: Output Parameters:
792: + inner - number of iterative refinement iterations (inner loop)
793: - blsize - number of iterative refinement iterations (blocksize loop)
795: Level: advanced
797: .seealso: NEPCISSSetRefinement()
798: @*/
799: PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
800: {
802: PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
803: PetscFunctionReturn(0);
804: }
806: static PetscErrorCode NEPCISSSetExtraction_CISS(NEP nep,NEPCISSExtraction extraction)
807: {
808: NEP_CISS *ctx = (NEP_CISS*)nep->data;
810: if (ctx->extraction != extraction) {
811: ctx->extraction = extraction;
812: nep->state = NEP_STATE_INITIAL;
813: }
814: PetscFunctionReturn(0);
815: }
817: /*@
818: NEPCISSSetExtraction - Sets the extraction technique used in the CISS solver.
820: Logically Collective on nep
822: Input Parameters:
823: + nep - the nonlinear eigensolver context
824: - extraction - the extraction technique
826: Options Database Key:
827: . -nep_ciss_extraction - Sets the extraction technique (either 'ritz', 'hankel' or 'caa')
829: Notes:
830: By default, the Rayleigh-Ritz extraction is used (NEP_CISS_EXTRACTION_RITZ).
832: If the 'hankel' or the 'caa' option is specified (NEP_CISS_EXTRACTION_HANKEL or
833: NEP_CISS_EXTRACTION_CAA), then the Block Hankel method, or the Communication-avoiding
834: Arnoldi method, respectively, is used for extracting eigenpairs.
836: Level: advanced
838: .seealso: NEPCISSGetExtraction(), NEPCISSExtraction
839: @*/
840: PetscErrorCode NEPCISSSetExtraction(NEP nep,NEPCISSExtraction extraction)
841: {
844: PetscTryMethod(nep,"NEPCISSSetExtraction_C",(NEP,NEPCISSExtraction),(nep,extraction));
845: PetscFunctionReturn(0);
846: }
848: static PetscErrorCode NEPCISSGetExtraction_CISS(NEP nep,NEPCISSExtraction *extraction)
849: {
850: NEP_CISS *ctx = (NEP_CISS*)nep->data;
852: *extraction = ctx->extraction;
853: PetscFunctionReturn(0);
854: }
856: /*@
857: NEPCISSGetExtraction - Gets the extraction technique used in the CISS solver.
859: Not Collective
861: Input Parameter:
862: . nep - the nonlinear eigensolver context
864: Output Parameters:
865: . extraction - extraction technique
867: Level: advanced
869: .seealso: NEPCISSSetExtraction() NEPCISSExtraction
870: @*/
871: PetscErrorCode NEPCISSGetExtraction(NEP nep,NEPCISSExtraction *extraction)
872: {
875: PetscUseMethod(nep,"NEPCISSGetExtraction_C",(NEP,NEPCISSExtraction*),(nep,extraction));
876: PetscFunctionReturn(0);
877: }
879: static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
880: {
881: NEP_CISS *ctx = (NEP_CISS*)nep->data;
882: SlepcContourData contour;
883: PetscInt i;
884: PC pc;
885: MPI_Comm child;
887: if (!ctx->contour) { /* initialize contour data structure first */
888: RGCanUseConjugates(nep->rg,ctx->isreal,&ctx->useconj);
889: SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)nep,&ctx->contour);
890: }
891: contour = ctx->contour;
892: if (!contour->ksp) {
893: PetscMalloc1(contour->npoints,&contour->ksp);
894: PetscSubcommGetChild(contour->subcomm,&child);
895: for (i=0;i<contour->npoints;i++) {
896: KSPCreate(child,&contour->ksp[i]);
897: PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)nep,1);
898: KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)nep)->prefix);
899: KSPAppendOptionsPrefix(contour->ksp[i],"nep_ciss_");
900: PetscLogObjectParent((PetscObject)nep,(PetscObject)contour->ksp[i]);
901: PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)nep)->options);
902: KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE);
903: KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
904: KSPGetPC(contour->ksp[i],&pc);
905: if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
906: KSPSetType(contour->ksp[i],KSPBCGS);
907: PCSetType(pc,PCBJACOBI);
908: } else {
909: KSPSetType(contour->ksp[i],KSPPREONLY);
910: PCSetType(pc,PCLU);
911: }
912: }
913: }
914: if (nsolve) *nsolve = contour->npoints;
915: if (ksp) *ksp = contour->ksp;
916: PetscFunctionReturn(0);
917: }
919: /*@C
920: NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
921: the CISS solver.
923: Not Collective
925: Input Parameter:
926: . nep - nonlinear eigenvalue solver
928: Output Parameters:
929: + nsolve - number of solver objects
930: - ksp - array of linear solver object
932: Notes:
933: The number of KSP solvers is equal to the number of integration points divided by
934: the number of partitions. This value is halved in the case of real matrices with
935: a region centered at the real axis.
937: Level: advanced
939: .seealso: NEPCISSSetSizes()
940: @*/
941: PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
942: {
944: PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
945: PetscFunctionReturn(0);
946: }
948: PetscErrorCode NEPReset_CISS(NEP nep)
949: {
950: NEP_CISS *ctx = (NEP_CISS*)nep->data;
952: BVDestroy(&ctx->S);
953: BVDestroy(&ctx->V);
954: BVDestroy(&ctx->Y);
955: SlepcContourDataReset(ctx->contour);
956: MatDestroy(&ctx->J);
957: BVDestroy(&ctx->pV);
958: if (ctx->extraction == NEP_CISS_EXTRACTION_RITZ && nep->fui==NEP_USER_INTERFACE_CALLBACK) PetscFree(ctx->dsctxf);
959: PetscFunctionReturn(0);
960: }
962: PetscErrorCode NEPSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,NEP nep)
963: {
964: NEP_CISS *ctx = (NEP_CISS*)nep->data;
965: PetscReal r1,r2;
966: PetscInt i,i1,i2,i3,i4,i5,i6,i7;
967: PetscBool b1,flg,flg2,flg3,flg4,flg5,flg6;
968: NEPCISSExtraction extraction;
970: PetscOptionsHead(PetscOptionsObject,"NEP CISS Options");
972: NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1);
973: PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,&flg);
974: PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,&flg2);
975: PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,&flg3);
976: PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,&flg4);
977: PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,&flg5);
978: PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,&flg6);
979: if (flg || flg2 || flg3 || flg4 || flg5 || flg6) NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1);
981: NEPCISSGetThreshold(nep,&r1,&r2);
982: PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,&flg);
983: PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,&flg2);
984: if (flg || flg2) NEPCISSSetThreshold(nep,r1,r2);
986: NEPCISSGetRefinement(nep,&i6,&i7);
987: PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,&flg);
988: PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,&flg2);
989: if (flg || flg2) NEPCISSSetRefinement(nep,i6,i7);
991: PetscOptionsEnum("-nep_ciss_extraction","Extraction technique","NEPCISSSetExtraction",NEPCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
992: if (flg) NEPCISSSetExtraction(nep,extraction);
994: PetscOptionsTail();
996: if (!nep->rg) NEPGetRG(nep,&nep->rg);
997: RGSetFromOptions(nep->rg); /* this is necessary here to set useconj */
998: if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
999: for (i=0;i<ctx->contour->npoints;i++) KSPSetFromOptions(ctx->contour->ksp[i]);
1000: PetscSubcommSetFromOptions(ctx->contour->subcomm);
1001: PetscFunctionReturn(0);
1002: }
1004: PetscErrorCode NEPDestroy_CISS(NEP nep)
1005: {
1006: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1008: SlepcContourDataDestroy(&ctx->contour);
1009: PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1010: PetscFree(nep->data);
1011: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL);
1012: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL);
1013: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL);
1014: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL);
1015: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL);
1016: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL);
1017: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NULL);
1018: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NULL);
1019: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL);
1020: PetscFunctionReturn(0);
1021: }
1023: PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1024: {
1025: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1026: PetscBool isascii;
1027: PetscViewer sviewer;
1029: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1030: if (isascii) {
1031: 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);
1032: if (ctx->isreal) PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n");
1033: PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1034: PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize);
1035: PetscViewerASCIIPrintf(viewer," extraction: %s\n",NEPCISSExtractions[ctx->extraction]);
1036: if (!ctx->contour || !ctx->contour->ksp) NEPCISSGetKSPs(nep,NULL,NULL);
1037: PetscViewerASCIIPushTab(viewer);
1038: if (ctx->npart>1 && ctx->contour->subcomm) {
1039: PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1040: if (!ctx->contour->subcomm->color) KSPView(ctx->contour->ksp[0],sviewer);
1041: PetscViewerFlush(sviewer);
1042: PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer);
1043: PetscViewerFlush(viewer);
1044: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1045: PetscViewerASCIIPopSynchronized(viewer);
1046: } else KSPView(ctx->contour->ksp[0],viewer);
1047: PetscViewerASCIIPopTab(viewer);
1048: }
1049: PetscFunctionReturn(0);
1050: }
1052: SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1053: {
1054: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1056: PetscNewLog(nep,&ctx);
1057: nep->data = ctx;
1058: /* set default values of parameters */
1059: ctx->N = 32;
1060: ctx->L = 16;
1061: ctx->M = ctx->N/4;
1062: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
1063: ctx->L_max = 64;
1064: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1065: ctx->isreal = PETSC_FALSE;
1066: ctx->npart = 1;
1068: nep->useds = PETSC_TRUE;
1070: nep->ops->solve = NEPSolve_CISS;
1071: nep->ops->setup = NEPSetUp_CISS;
1072: nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1073: nep->ops->reset = NEPReset_CISS;
1074: nep->ops->destroy = NEPDestroy_CISS;
1075: nep->ops->view = NEPView_CISS;
1077: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS);
1078: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS);
1079: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS);
1080: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS);
1081: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS);
1082: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS);
1083: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetExtraction_C",NEPCISSSetExtraction_CISS);
1084: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetExtraction_C",NEPCISSGetExtraction_CISS);
1085: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS);
1086: PetscFunctionReturn(0);
1087: }