Actual source code: nciss.c
slepc-3.13.2 2020-05-12
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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>
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] = PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
150: ctx->weight[i] = radius*PetscCMPLX(vscale*PetscCosReal(theta),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: {
318: NEP_CISS *ctx = (NEP_CISS*)nep->data;
319: PetscInt i,ml=ctx->L*ctx->M;
320: PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
321: PetscScalar *work;
322: #if defined(PETSC_USE_COMPLEX)
323: PetscReal *rwork;
324: #endif
327: PetscMalloc1(5*ml,&work);
328: #if defined(PETSC_USE_COMPLEX)
329: PetscMalloc1(5*ml,&rwork);
330: #endif
331: PetscBLASIntCast(ml,&m);
332: n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
333: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
334: #if defined(PETSC_USE_COMPLEX)
335: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
336: #else
337: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
338: #endif
339: SlepcCheckLapackInfo("gesvd",info);
340: PetscFPTrapPop();
341: (*K) = 0;
342: for (i=0;i<ml;i++) {
343: if (ctx->sigma[i]/PetscMax(ctx->sigma[0],1)>ctx->delta) (*K)++;
344: }
345: PetscFree(work);
346: #if defined(PETSC_USE_COMPLEX)
347: PetscFree(rwork);
348: #endif
349: return(0);
350: }
352: static PetscErrorCode ConstructS(NEP nep)
353: {
355: NEP_CISS *ctx = (NEP_CISS*)nep->data;
356: PetscInt i,j,k,vec_local_size,p_id;
357: Vec v,sj,yj;
358: PetscScalar *ppk, *v_data, m = 1;
361: BVGetSizes(ctx->Y,&vec_local_size,NULL,NULL);
362: PetscMalloc1(ctx->num_solve_point,&ppk);
363: for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
364: BVGetColumn(ctx->Y,0,&yj);
365: VecDuplicate(yj,&v);
366: BVRestoreColumn(ctx->Y,0,&yj);
367: for (k=0;k<ctx->M;k++) {
368: for (j=0;j<ctx->L;j++) {
369: VecSet(v,0);
370: for (i=0;i<ctx->num_solve_point;i++) {
371: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
372: BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
373: BVMultVec(ctx->Y,ppk[i]*ctx->weight[p_id],1,v,&m);
374: }
375: if (ctx->useconj) {
376: VecGetArray(v,&v_data);
377: for (i=0;i<vec_local_size;i++) v_data[i] = PetscRealPart(v_data[i])*2;
378: VecRestoreArray(v,&v_data);
379: }
380: BVGetColumn(ctx->S,k*ctx->L+j,&sj);
381: VecCopy(v,sj);
382: BVRestoreColumn(ctx->S,k*ctx->L+j,&sj);
383: }
384: for (i=0;i<ctx->num_solve_point;i++) {
385: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
386: ppk[i] *= ctx->pp[p_id];
387: }
388: }
389: PetscFree(ppk);
390: VecDestroy(&v);
391: return(0);
392: }
394: static PetscErrorCode isGhost(NEP nep,PetscInt ld,PetscInt nv,PetscBool *fl)
395: {
397: NEP_CISS *ctx = (NEP_CISS*)nep->data;
398: PetscInt i,j;
399: PetscScalar *pX;
400: PetscReal *tau,s1,s2,tau_max=0.0;
403: PetscMalloc1(nv,&tau);
404: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
405: DSGetArray(nep->ds,DS_MAT_X,&pX);
407: for (i=0;i<nv;i++) {
408: s1 = 0;
409: s2 = 0;
410: for (j=0;j<nv;j++) {
411: s1 += PetscAbsScalar(PetscPowScalar(pX[i*ld+j],2));
412: s2 += PetscPowReal(PetscAbsScalar(pX[i*ld+j]),2)/ctx->sigma[j];
413: }
414: tau[i] = s1/s2;
415: tau_max = PetscMax(tau_max,tau[i]);
416: }
417: DSRestoreArray(nep->ds,DS_MAT_X,&pX);
418: for (i=0;i<nv;i++) {
419: tau[i] /= tau_max;
420: }
421: for (i=0;i<nv;i++) {
422: if (tau[i]>=ctx->spurious_threshold) fl[i] = PETSC_TRUE;
423: else fl[i] = PETSC_FALSE;
424: }
425: PetscFree(tau);
426: return(0);
427: }
429: PetscErrorCode NEPSetUp_CISS(NEP nep)
430: {
432: NEP_CISS *ctx = (NEP_CISS*)nep->data;
433: PetscInt nwork;
434: PetscBool istrivial,isellipse,flg,useconj;
437: if (!nep->ncv) {
438: nep->ncv = ctx->L_max*ctx->M;
439: if (nep->ncv>nep->n) {
440: nep->ncv = nep->n;
441: ctx->L_max = nep->ncv/ctx->M;
442: if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
443: }
444: } else {
445: ctx->L_max = nep->ncv/ctx->M;
446: if (!ctx->L_max) {
447: ctx->L_max = 1;
448: nep->ncv = ctx->L_max*ctx->M;
449: }
450: }
451: ctx->L = PetscMin(ctx->L,ctx->L_max);
452: if (!nep->max_it) nep->max_it = 1;
453: if (!nep->mpd) nep->mpd = nep->ncv;
454: if (!nep->which) nep->which = NEP_ALL;
455: if (nep->stopping!=NEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support user-defined stopping test");
457: /* check region */
458: RGIsTrivial(nep->rg,&istrivial);
459: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
460: RGGetComplement(nep->rg,&flg);
461: if (flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
462: PetscObjectTypeCompare((PetscObject)nep->rg,RGELLIPSE,&isellipse);
463: if (!isellipse) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Currently only implemented for elliptic regions");
464: /* if useconj has changed, then reset subcomm data */
465: NEPCISSSetUseConj(nep,&useconj);
466: if (useconj!=ctx->useconj) { NEPCISSResetSubcomm(nep); }
468: /* create split comm */
469: if (!ctx->subcomm) { NEPCISSSetUpSubComm(nep,&ctx->num_solve_point); }
471: NEPAllocateSolution(nep,0);
472: if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
473: PetscMalloc4(ctx->N,&ctx->weight,ctx->N,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
474: PetscLogObjectMemory((PetscObject)nep,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));
476: /* allocate basis vectors */
477: BVDestroy(&ctx->S);
478: BVDuplicateResize(nep->V,ctx->L_max*ctx->M,&ctx->S);
479: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->S);
480: BVDestroy(&ctx->V);
481: BVDuplicateResize(nep->V,ctx->L_max,&ctx->V);
482: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->V);
484: BVDestroy(&ctx->Y);
485: BVDuplicateResize(nep->V,ctx->num_solve_point*ctx->L_max,&ctx->Y);
487: DSSetType(nep->ds,DSGNHEP);
488: DSAllocate(nep->ds,nep->ncv);
489: nwork = (nep->fui==NEP_USER_INTERFACE_SPLIT)? 2: 1;
490: NEPSetWorkVecs(nep,nwork);
491: return(0);
492: }
494: PetscErrorCode NEPSolve_CISS(NEP nep)
495: {
497: NEP_CISS *ctx = (NEP_CISS*)nep->data;
498: Mat X,M;
499: PetscInt i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
500: PetscScalar *Mu,*H0,*H1,*rr,*temp,center;
501: PetscReal error,max_error,radius,rgscale;
502: PetscBool *fl1;
503: Vec si;
504: SlepcSC sc;
505: PetscRandom rand;
508: DSGetSlepcSC(nep->ds,&sc);
509: sc->comparison = SlepcCompareLargestMagnitude;
510: sc->comparisonctx = NULL;
511: sc->map = NULL;
512: sc->mapobj = NULL;
513: DSGetLeadingDimension(nep->ds,&ld);
514: SetPathParameter(nep);
515: CISSVecSetRandom(ctx->V,0,ctx->L);
516: BVGetRandomContext(ctx->V,&rand);
518: SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_TRUE);
519: EstimateNumberEigs(nep,&L_add);
520: if (L_add>0) {
521: PetscInfo2(nep,"Changing L %D -> %D by Estimate #Eig\n",ctx->L,ctx->L+L_add);
522: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
523: SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
524: ctx->L += L_add;
525: }
526: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
527: for (i=0;i<ctx->refine_blocksize;i++) {
528: CalcMu(nep,Mu);
529: BlockHankel(nep,Mu,0,H0);
530: SVD_H0(nep,H0,&nv);
531: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
532: L_add = L_base;
533: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
534: PetscInfo2(nep,"Changing L %D -> %D by SVD(H0)\n",ctx->L,ctx->L+L_add);
535: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
536: SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
537: ctx->L += L_add;
538: }
539: PetscFree2(Mu,H0);
541: RGGetScale(nep->rg,&rgscale);
542: RGEllipseGetParameters(nep->rg,¢er,&radius,NULL);
544: 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);
545: while (nep->reason == NEP_CONVERGED_ITERATING) {
546: nep->its++;
547: for (inner=0;inner<=ctx->refine_inner;inner++) {
548: CalcMu(nep,Mu);
549: BlockHankel(nep,Mu,0,H0);
550: SVD_H0(nep,H0,&nv);
551: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
552: ConstructS(nep);
553: BVSetActiveColumns(ctx->S,0,ctx->L);
554: BVCopy(ctx->S,ctx->V);
555: SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_FALSE);
556: } else break;
557: }
559: nep->nconv = 0;
560: if (nv == 0) { nep->reason = NEP_CONVERGED_TOL; break; }
561: BlockHankel(nep,Mu,0,H0);
562: BlockHankel(nep,Mu,1,H1);
563: DSSetDimensions(nep->ds,nv,0,0,0);
564: DSSetState(nep->ds,DS_STATE_RAW);
565: DSGetArray(nep->ds,DS_MAT_A,&temp);
566: for (j=0;j<nv;j++)
567: for (i=0;i<nv;i++)
568: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
569: DSRestoreArray(nep->ds,DS_MAT_A,&temp);
570: DSGetArray(nep->ds,DS_MAT_B,&temp);
571: for (j=0;j<nv;j++)
572: for (i=0;i<nv;i++)
573: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
574: DSRestoreArray(nep->ds,DS_MAT_B,&temp);
575: DSSolve(nep->ds,nep->eigr,nep->eigi);
576: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
577: for (i=0;i<nv;i++){
578: nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
579: #if !defined(PETSC_USE_COMPLEX)
580: nep->eigi[i] = nep->eigi[i]*radius*rgscale;
581: #endif
582: }
583: PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
584: isGhost(nep,ld,nv,fl1);
585: RGCheckInside(nep->rg,nv,nep->eigr,nep->eigi,inside);
586: for (i=0;i<nv;i++) {
587: if (fl1[i] && inside[i]>0) {
588: rr[i] = 1.0;
589: nep->nconv++;
590: } else rr[i] = 0.0;
591: }
592: DSSort(nep->ds,nep->eigr,nep->eigi,rr,NULL,&nep->nconv);
593: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
594: for (i=0;i<nv;i++){
595: nep->eigr[i] = (nep->eigr[i]*radius+center)*rgscale;
596: #if !defined(PETSC_USE_COMPLEX)
597: nep->eigi[i] = nep->eigi[i]*radius*rgscale;
598: #endif
599: }
600: PetscFree3(fl1,inside,rr);
601: BVSetActiveColumns(nep->V,0,nv);
602: ConstructS(nep);
603: BVSetActiveColumns(ctx->S,0,nv);
604: BVCopy(ctx->S,nep->V);
606: DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
607: DSGetMat(nep->ds,DS_MAT_X,&X);
608: BVMultInPlace(ctx->S,X,0,nep->nconv);
609: BVMultInPlace(nep->V,X,0,nep->nconv);
610: MatDestroy(&X);
611: max_error = 0.0;
612: for (i=0;i<nep->nconv;i++) {
613: BVGetColumn(nep->V,i,&si);
614: VecNormalize(si,NULL);
615: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,nep->eigr[i],si,nep->work,&error);
616: (*nep->converged)(nep,nep->eigr[i],0,error,&error,nep->convergedctx);
617: BVRestoreColumn(nep->V,i,&si);
618: max_error = PetscMax(max_error,error);
619: }
620: if (max_error <= nep->tol) nep->reason = NEP_CONVERGED_TOL;
621: else if (nep->its > nep->max_it) nep->reason = NEP_DIVERGED_ITS;
622: else {
623: if (nep->nconv > ctx->L) nv = nep->nconv;
624: else if (ctx->L > nv) nv = ctx->L;
625: MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M);
626: MatDenseGetArray(M,&temp);
627: for (i=0;i<ctx->L*nv;i++) {
628: PetscRandomGetValue(rand,&temp[i]);
629: temp[i] = PetscRealPart(temp[i]);
630: }
631: MatDenseRestoreArray(M,&temp);
632: BVSetActiveColumns(ctx->S,0,nv);
633: BVMultInPlace(ctx->S,M,0,ctx->L);
634: MatDestroy(&M);
635: BVSetActiveColumns(ctx->S,0,ctx->L);
636: BVCopy(ctx->S,ctx->V);
637: SolveLinearSystem(nep,nep->function,nep->jacobian,ctx->V,0,ctx->L,PETSC_FALSE);
638: }
639: }
640: PetscFree3(Mu,H0,H1);
641: return(0);
642: }
644: static PetscErrorCode NEPCISSSetSizes_CISS(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
645: {
647: NEP_CISS *ctx = (NEP_CISS*)nep->data;
648: PetscInt oN,onpart;
651: oN = ctx->N;
652: if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
653: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
654: } else {
655: if (ip<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
656: if (ip%2) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
657: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
658: }
659: if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
660: ctx->L = 16;
661: } else {
662: if (bs<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
663: ctx->L = bs;
664: }
665: if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
666: ctx->M = ctx->N/4;
667: } else {
668: if (ms<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
669: 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");
670: ctx->M = ms;
671: }
672: onpart = ctx->npart;
673: if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
674: ctx->npart = 1;
675: } else {
676: if (npart<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
677: ctx->npart = npart;
678: if (npart>1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Current implementation does not allow partitions");
679: }
680: if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
681: ctx->L_max = 64;
682: } else {
683: if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
684: ctx->L_max = PetscMax(bsmax,ctx->L);
685: }
686: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) { NEPCISSResetSubcomm(nep); }
687: ctx->isreal = realmats;
688: nep->state = NEP_STATE_INITIAL;
689: return(0);
690: }
692: /*@
693: NEPCISSSetSizes - Sets the values of various size parameters in the CISS solver.
695: Logically Collective on nep
697: Input Parameters:
698: + nep - the eigenproblem solver context
699: . ip - number of integration points
700: . bs - block size
701: . ms - moment size
702: . npart - number of partitions when splitting the communicator
703: . bsmax - max block size
704: - realmats - T(z) is real for real z
706: Options Database Keys:
707: + -nep_ciss_integration_points - Sets the number of integration points
708: . -nep_ciss_blocksize - Sets the block size
709: . -nep_ciss_moments - Sets the moment size
710: . -nep_ciss_partitions - Sets the number of partitions
711: . -nep_ciss_maxblocksize - Sets the maximum block size
712: - -nep_ciss_realmats - T(z) is real for real z
714: Note:
715: The default number of partitions is 1. This means the internal KSP object is shared
716: among all processes of the NEP communicator. Otherwise, the communicator is split
717: into npart communicators, so that npart KSP solves proceed simultaneously.
719: The realmats flag can be set to true when T(.) is guaranteed to be real
720: when the argument is a real value, for example, when all matrices in
721: the split form are real. When set to true, the solver avoids some computations.
723: Level: advanced
725: .seealso: NEPCISSGetSizes()
726: @*/
727: PetscErrorCode NEPCISSSetSizes(NEP nep,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
728: {
739: PetscTryMethod(nep,"NEPCISSSetSizes_C",(NEP,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(nep,ip,bs,ms,npart,bsmax,realmats));
740: return(0);
741: }
743: static PetscErrorCode NEPCISSGetSizes_CISS(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
744: {
745: NEP_CISS *ctx = (NEP_CISS*)nep->data;
748: if (ip) *ip = ctx->N;
749: if (bs) *bs = ctx->L;
750: if (ms) *ms = ctx->M;
751: if (npart) *npart = ctx->npart;
752: if (bsmax) *bsmax = ctx->L_max;
753: if (realmats) *realmats = ctx->isreal;
754: return(0);
755: }
757: /*@
758: NEPCISSGetSizes - Gets the values of various size parameters in the CISS solver.
760: Not Collective
762: Input Parameter:
763: . nep - the eigenproblem solver context
765: Output Parameters:
766: + ip - number of integration points
767: . bs - block size
768: . ms - moment size
769: . npart - number of partitions when splitting the communicator
770: . bsmax - max block size
771: - realmats - T(z) is real for real z
773: Level: advanced
775: .seealso: NEPCISSSetSizes()
776: @*/
777: PetscErrorCode NEPCISSGetSizes(NEP nep,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
778: {
783: PetscUseMethod(nep,"NEPCISSGetSizes_C",(NEP,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(nep,ip,bs,ms,npart,bsmax,realmats));
784: return(0);
785: }
787: static PetscErrorCode NEPCISSSetThreshold_CISS(NEP nep,PetscReal delta,PetscReal spur)
788: {
789: NEP_CISS *ctx = (NEP_CISS*)nep->data;
792: if (delta == PETSC_DEFAULT) {
793: ctx->delta = 1e-12;
794: } else {
795: if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
796: ctx->delta = delta;
797: }
798: if (spur == PETSC_DEFAULT) {
799: ctx->spurious_threshold = 1e-4;
800: } else {
801: if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
802: ctx->spurious_threshold = spur;
803: }
804: return(0);
805: }
807: /*@
808: NEPCISSSetThreshold - Sets the values of various threshold parameters in
809: the CISS solver.
811: Logically Collective on nep
813: Input Parameters:
814: + nep - the eigenproblem solver context
815: . delta - threshold for numerical rank
816: - spur - spurious threshold (to discard spurious eigenpairs)
818: Options Database Keys:
819: + -nep_ciss_delta - Sets the delta
820: - -nep_ciss_spurious_threshold - Sets the spurious threshold
822: Level: advanced
824: .seealso: NEPCISSGetThreshold()
825: @*/
826: PetscErrorCode NEPCISSSetThreshold(NEP nep,PetscReal delta,PetscReal spur)
827: {
834: PetscTryMethod(nep,"NEPCISSSetThreshold_C",(NEP,PetscReal,PetscReal),(nep,delta,spur));
835: return(0);
836: }
838: static PetscErrorCode NEPCISSGetThreshold_CISS(NEP nep,PetscReal *delta,PetscReal *spur)
839: {
840: NEP_CISS *ctx = (NEP_CISS*)nep->data;
843: if (delta) *delta = ctx->delta;
844: if (spur) *spur = ctx->spurious_threshold;
845: return(0);
846: }
848: /*@
849: NEPCISSGetThreshold - Gets the values of various threshold parameters in
850: the CISS solver.
852: Not Collective
854: Input Parameter:
855: . nep - the eigenproblem solver context
857: Output Parameters:
858: + delta - threshold for numerical rank
859: - spur - spurious threshold (to discard spurious eigenpairs)
861: Level: advanced
863: .seealso: NEPCISSSetThreshold()
864: @*/
865: PetscErrorCode NEPCISSGetThreshold(NEP nep,PetscReal *delta,PetscReal *spur)
866: {
871: PetscUseMethod(nep,"NEPCISSGetThreshold_C",(NEP,PetscReal*,PetscReal*),(nep,delta,spur));
872: return(0);
873: }
875: static PetscErrorCode NEPCISSSetRefinement_CISS(NEP nep,PetscInt inner,PetscInt blsize)
876: {
877: NEP_CISS *ctx = (NEP_CISS*)nep->data;
880: if (inner == PETSC_DEFAULT) {
881: ctx->refine_inner = 0;
882: } else {
883: if (inner<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
884: ctx->refine_inner = inner;
885: }
886: if (blsize == PETSC_DEFAULT) {
887: ctx->refine_blocksize = 0;
888: } else {
889: if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
890: ctx->refine_blocksize = blsize;
891: }
892: return(0);
893: }
895: /*@
896: NEPCISSSetRefinement - Sets the values of various refinement parameters
897: in the CISS solver.
899: Logically Collective on nep
901: Input Parameters:
902: + nep - the eigenproblem solver context
903: . inner - number of iterative refinement iterations (inner loop)
904: - blsize - number of iterative refinement iterations (blocksize loop)
906: Options Database Keys:
907: + -nep_ciss_refine_inner - Sets number of inner iterations
908: - -nep_ciss_refine_blocksize - Sets number of blocksize iterations
910: Level: advanced
912: .seealso: NEPCISSGetRefinement()
913: @*/
914: PetscErrorCode NEPCISSSetRefinement(NEP nep,PetscInt inner,PetscInt blsize)
915: {
922: PetscTryMethod(nep,"NEPCISSSetRefinement_C",(NEP,PetscInt,PetscInt),(nep,inner,blsize));
923: return(0);
924: }
926: static PetscErrorCode NEPCISSGetRefinement_CISS(NEP nep,PetscInt *inner,PetscInt *blsize)
927: {
928: NEP_CISS *ctx = (NEP_CISS*)nep->data;
931: if (inner) *inner = ctx->refine_inner;
932: if (blsize) *blsize = ctx->refine_blocksize;
933: return(0);
934: }
936: /*@
937: NEPCISSGetRefinement - Gets the values of various refinement parameters
938: in the CISS solver.
940: Not Collective
942: Input Parameter:
943: . nep - the eigenproblem solver context
945: Output Parameters:
946: + inner - number of iterative refinement iterations (inner loop)
947: - blsize - number of iterative refinement iterations (blocksize loop)
949: Level: advanced
951: .seealso: NEPCISSSetRefinement()
952: @*/
953: PetscErrorCode NEPCISSGetRefinement(NEP nep, PetscInt *inner, PetscInt *blsize)
954: {
959: PetscUseMethod(nep,"NEPCISSGetRefinement_C",(NEP,PetscInt*,PetscInt*),(nep,inner,blsize));
960: return(0);
961: }
963: static PetscErrorCode NEPCISSGetKSPs_CISS(NEP nep,PetscInt *nsolve,KSP **ksp)
964: {
966: NEP_CISS *ctx = (NEP_CISS*)nep->data;
967: PetscInt i;
968: PC pc;
971: if (!ctx->ksp) {
972: if (!ctx->subcomm) { /* initialize subcomm first */
973: NEPCISSSetUseConj(nep,&ctx->useconj);
974: NEPCISSSetUpSubComm(nep,&ctx->num_solve_point);
975: }
976: PetscMalloc1(ctx->num_solve_point,&ctx->ksp);
977: for (i=0;i<ctx->num_solve_point;i++) {
978: KSPCreate(PetscSubcommChild(ctx->subcomm),&ctx->ksp[i]);
979: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
980: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
981: KSPAppendOptionsPrefix(ctx->ksp[i],"nep_ciss_");
982: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
983: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
984: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
985: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
986: KSPGetPC(ctx->ksp[i],&pc);
987: KSPSetType(ctx->ksp[i],KSPPREONLY);
988: PCSetType(pc,PCLU);
989: }
990: }
991: if (nsolve) *nsolve = ctx->num_solve_point;
992: if (ksp) *ksp = ctx->ksp;
993: return(0);
994: }
996: /*@C
997: NEPCISSGetKSPs - Retrieve the array of linear solver objects associated with
998: the CISS solver.
1000: Not Collective
1002: Input Parameter:
1003: . nep - nonlinear eigenvalue solver
1005: Output Parameters:
1006: + nsolve - number of solver objects
1007: - ksp - array of linear solver object
1009: Notes:
1010: The number of KSP solvers is equal to the number of integration points divided by
1011: the number of partitions. This value is halved in the case of real matrices with
1012: a region centered at the real axis.
1014: Level: advanced
1016: .seealso: NEPCISSSetSizes()
1017: @*/
1018: PetscErrorCode NEPCISSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1019: {
1024: PetscUseMethod(nep,"NEPCISSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1025: return(0);
1026: }
1028: PetscErrorCode NEPReset_CISS(NEP nep)
1029: {
1031: PetscInt i;
1032: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1035: BVDestroy(&ctx->S);
1036: BVDestroy(&ctx->V);
1037: BVDestroy(&ctx->Y);
1038: for (i=0;i<ctx->num_solve_point;i++) {
1039: KSPReset(ctx->ksp[i]);
1040: }
1041: return(0);
1042: }
1044: PetscErrorCode NEPSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,NEP nep)
1045: {
1047: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1048: PetscReal r1,r2;
1049: PetscInt i,i1,i2,i3,i4,i5,i6,i7;
1050: PetscBool b1;
1053: PetscOptionsHead(PetscOptionsObject,"NEP CISS Options");
1055: NEPCISSGetSizes(nep,&i1,&i2,&i3,&i4,&i5,&b1);
1056: PetscOptionsInt("-nep_ciss_integration_points","Number of integration points","NEPCISSSetSizes",i1,&i1,NULL);
1057: PetscOptionsInt("-nep_ciss_blocksize","Block size","NEPCISSSetSizes",i2,&i2,NULL);
1058: PetscOptionsInt("-nep_ciss_moments","Moment size","NEPCISSSetSizes",i3,&i3,NULL);
1059: PetscOptionsInt("-nep_ciss_partitions","Number of partitions","NEPCISSSetSizes",i4,&i4,NULL);
1060: PetscOptionsInt("-nep_ciss_maxblocksize","Maximum block size","NEPCISSSetSizes",i5,&i5,NULL);
1061: PetscOptionsBool("-nep_ciss_realmats","True if T(z) is real for real z","NEPCISSSetSizes",b1,&b1,NULL);
1062: NEPCISSSetSizes(nep,i1,i2,i3,i4,i5,b1);
1064: NEPCISSGetThreshold(nep,&r1,&r2);
1065: PetscOptionsReal("-nep_ciss_delta","Threshold for numerical rank","NEPCISSSetThreshold",r1,&r1,NULL);
1066: PetscOptionsReal("-nep_ciss_spurious_threshold","Threshold for the spurious eigenpairs","NEPCISSSetThreshold",r2,&r2,NULL);
1067: NEPCISSSetThreshold(nep,r1,r2);
1069: NEPCISSGetRefinement(nep,&i6,&i7);
1070: PetscOptionsInt("-nep_ciss_refine_inner","Number of inner iterative refinement iterations","NEPCISSSetRefinement",i6,&i6,NULL);
1071: PetscOptionsInt("-nep_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","NEPCISSSetRefinement",i7,&i7,NULL);
1072: NEPCISSSetRefinement(nep,i6,i7);
1074: PetscOptionsTail();
1076: if (!nep->rg) { NEPGetRG(nep,&nep->rg); }
1077: RGSetFromOptions(nep->rg); /* this is necessary here to set useconj */
1078: if (!ctx->ksp) { NEPCISSGetKSPs(nep,&ctx->num_solve_point,&ctx->ksp); }
1079: for (i=0;i<ctx->num_solve_point;i++) {
1080: KSPSetFromOptions(ctx->ksp[i]);
1081: }
1082: PetscSubcommSetFromOptions(ctx->subcomm);
1083: return(0);
1084: }
1086: PetscErrorCode NEPDestroy_CISS(NEP nep)
1087: {
1089: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1092: NEPCISSResetSubcomm(nep);
1093: PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1094: PetscFree(nep->data);
1095: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NULL);
1096: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NULL);
1097: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NULL);
1098: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NULL);
1099: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NULL);
1100: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NULL);
1101: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NULL);
1102: return(0);
1103: }
1105: PetscErrorCode NEPView_CISS(NEP nep,PetscViewer viewer)
1106: {
1108: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1109: PetscBool isascii;
1112: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1113: if (isascii) {
1114: 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);
1115: if (ctx->isreal) {
1116: PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n");
1117: }
1118: PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1119: PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1120: if (!ctx->ksp) { NEPCISSGetKSPs(nep,&ctx->num_solve_point,&ctx->ksp); }
1121: PetscViewerASCIIPushTab(viewer);
1122: KSPView(ctx->ksp[0],viewer);
1123: PetscViewerASCIIPopTab(viewer);
1124: }
1125: return(0);
1126: }
1128: SLEPC_EXTERN PetscErrorCode NEPCreate_CISS(NEP nep)
1129: {
1131: NEP_CISS *ctx = (NEP_CISS*)nep->data;
1134: PetscNewLog(nep,&ctx);
1135: nep->data = ctx;
1136: /* set default values of parameters */
1137: ctx->N = 32;
1138: ctx->L = 16;
1139: ctx->M = ctx->N/4;
1140: ctx->delta = 1e-12;
1141: ctx->L_max = 64;
1142: ctx->spurious_threshold = 1e-4;
1143: ctx->isreal = PETSC_FALSE;
1144: ctx->npart = 1;
1146: nep->useds = PETSC_TRUE;
1148: nep->ops->solve = NEPSolve_CISS;
1149: nep->ops->setup = NEPSetUp_CISS;
1150: nep->ops->setfromoptions = NEPSetFromOptions_CISS;
1151: nep->ops->reset = NEPReset_CISS;
1152: nep->ops->destroy = NEPDestroy_CISS;
1153: nep->ops->view = NEPView_CISS;
1155: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetSizes_C",NEPCISSSetSizes_CISS);
1156: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetSizes_C",NEPCISSGetSizes_CISS);
1157: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetThreshold_C",NEPCISSSetThreshold_CISS);
1158: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetThreshold_C",NEPCISSGetThreshold_CISS);
1159: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSSetRefinement_C",NEPCISSSetRefinement_CISS);
1160: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetRefinement_C",NEPCISSGetRefinement_CISS);
1161: PetscObjectComposeFunction((PetscObject)nep,"NEPCISSGetKSPs_C",NEPCISSGetKSPs_CISS);
1162: return(0);
1163: }