Actual source code: ciss.c
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] 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/epsimpl.h>
34: #include <slepc/private/slepccontour.h>
35: #include <slepcblaslapack.h>
37: typedef struct {
38: /* user parameters */
39: PetscInt N; /* number of integration points (32) */
40: PetscInt L; /* block size (16) */
41: PetscInt M; /* moment degree (N/4 = 4) */
42: PetscReal delta; /* threshold of singular value (1e-12) */
43: PetscInt L_max; /* maximum number of columns of the source matrix V */
44: PetscReal spurious_threshold; /* discard spurious eigenpairs */
45: PetscBool isreal; /* A and B are real */
46: PetscInt npart; /* number of partitions */
47: PetscInt refine_inner;
48: PetscInt refine_blocksize;
49: EPSCISSQuadRule quad;
50: EPSCISSExtraction extraction;
51: EPSCISSStrategy strategy;
52: /* private data */
53: SlepcContourData contour;
54: PetscReal *sigma; /* threshold for numerical rank */
55: PetscScalar *weight;
56: PetscScalar *omega;
57: PetscScalar *pp;
58: BV V;
59: BV S;
60: BV pV;
61: BV Y;
62: PetscBool useconj;
63: PetscObjectId rgid;
64: PetscObjectState rgstate;
65: } EPS_CISS;
67: /* initialize contour data structure */
68: static PetscErrorCode EPSCISSGetContour_Private(EPS eps,SlepcContourData *contour)
69: {
70: EPS_CISS *ctx = (EPS_CISS*)eps->data;
72: PetscFunctionBegin;
73: if (!ctx->contour) {
74: PetscCall(RGCanUseConjugates(eps->rg,ctx->isreal,&ctx->useconj));
75: PetscCall(SlepcContourDataCreate(ctx->useconj?ctx->N/2:ctx->N,ctx->npart,(PetscObject)eps,&ctx->contour));
76: }
77: if (contour) *contour = ctx->contour;
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
80: /*
81: Set up KSP solvers for SPLIT strategy
82: */
83: static PetscErrorCode EPSCISSSetUp_SPLIT(EPS eps,Mat A,Mat B,Mat Pa,Mat Pb)
84: {
85: EPS_CISS *ctx = (EPS_CISS*)eps->data;
86: SlepcContourData contour;
87: PetscInt i,p_id,nsplit;
88: Mat Amat,Pmat;
89: MatStructure str,strp;
91: PetscFunctionBegin;
92: PetscCall(EPSCISSGetContour_Private(eps,&contour));
93: PetscCall(STGetMatStructure(eps->st,&str));
94: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,&strp));
95: for (i=0;i<contour->npoints;i++) {
96: p_id = i*contour->subcomm->n + contour->subcomm->color;
97: PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Amat));
98: if (B) PetscCall(MatAXPY(Amat,-ctx->omega[p_id],B,str));
99: else PetscCall(MatShift(Amat,-ctx->omega[p_id]));
100: if (nsplit) {
101: PetscCall(MatDuplicate(Pa,MAT_COPY_VALUES,&Pmat));
102: if (Pb) PetscCall(MatAXPY(Pmat,-ctx->omega[p_id],Pb,strp));
103: else PetscCall(MatShift(Pmat,-ctx->omega[p_id]));
104: } else Pmat = Amat;
105: PetscCall(EPS_KSPSetOperators(contour->ksp[i],Amat,Pmat));
106: if (eps->setfromoptionscalled) PetscCall(KSPSetFromOptions(contour->ksp[i]));
107: PetscCall(KSPSetUp(contour->ksp[i]));
108: PetscCall(MatDestroy(&Amat));
109: if (nsplit) PetscCall(MatDestroy(&Pmat));
110: }
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: /*
115: Set up KSP solvers for MULTISHIFT strategy
116: */
117: static PetscErrorCode EPSCISSSetUp_MULTISHIFT(EPS eps,Mat A,Mat B,Mat Pa,Mat Pb)
118: {
119: EPS_CISS *ctx = (EPS_CISS*)eps->data;
120: SlepcContourData contour;
121: PetscInt i,p_id,nsplit,nshift;
122: Mat Amat,Pmat;
123: MatStructure str,strp;
124: PetscScalar *sigma,*sigma_imaginary=NULL;
125: PetscBool explicitmat = PETSC_FALSE; // TODO: add user option
127: PetscFunctionBegin;
128: PetscCall(EPSCISSGetContour_Private(eps,&contour));
129: nshift = contour->npoints;
130: PetscCall(PetscCalloc2(nshift,&sigma,nshift,&sigma_imaginary));
131: for (i=0;i<nshift;i++) {
132: p_id = i*contour->subcomm->n + contour->subcomm->color;
133: sigma[i] = -ctx->omega[p_id];
134: }
135: PetscCall(STGetMatStructure(eps->st,&str));
136: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,&strp));
137: PetscCall(MatCreateNestFromMultipleShifts(A,nshift,sigma,sigma_imaginary,B,explicitmat,str,&Amat));
138: if (nsplit) PetscCall(MatCreateNestFromMultipleShifts(Pa,nshift,sigma,sigma_imaginary,Pb,explicitmat,strp,&Pmat));
139: else Pmat = Amat;
140: PetscCall(EPS_KSPSetOperators(contour->ksp[0],Amat,Pmat));
141: if (eps->setfromoptionscalled) PetscCall(KSPSetFromOptions(contour->ksp[0]));
142: PetscCall(KSPSetUp(contour->ksp[0]));
143: PetscCall(MatDestroy(&Amat));
144: if (nsplit) PetscCall(MatDestroy(&Pmat));
145: PetscCall(PetscFree2(sigma,sigma_imaginary));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*
150: Linear solves for the USEST strategy
151: */
152: static PetscErrorCode EPSCISSSolve_USEST(EPS eps,Mat V,PetscInt L_start,PetscInt L_end)
153: {
154: EPS_CISS *ctx = (EPS_CISS*)eps->data;
155: SlepcContourData contour;
156: PetscInt i,p_id;
157: Mat MC;
158: KSP ksp;
160: PetscFunctionBegin;
161: PetscCall(EPSCISSGetContour_Private(eps,&contour));
162: for (i=0;i<contour->npoints;i++) {
163: p_id = i*contour->subcomm->n + contour->subcomm->color;
164: PetscCall(STSetShift(eps->st,ctx->omega[p_id]));
165: PetscCall(STGetKSP(eps->st,&ksp));
166: PetscCall(BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end));
167: PetscCall(BVGetMat(ctx->Y,&MC));
168: PetscCall(KSPMatSolve(ksp,V,MC));
169: PetscCall(BVRestoreMat(ctx->Y,&MC));
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*
175: Linear solves for the SPLIT strategy
176: */
177: static PetscErrorCode EPSCISSSolve_SPLIT(EPS eps,Mat V,PetscInt L_start,PetscInt L_end)
178: {
179: EPS_CISS *ctx = (EPS_CISS*)eps->data;
180: SlepcContourData contour;
181: PetscInt i;
182: Mat MC;
183: KSP ksp;
185: PetscFunctionBegin;
186: PetscCall(EPSCISSGetContour_Private(eps,&contour));
187: for (i=0;i<contour->npoints;i++) {
188: PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
189: ksp = contour->ksp[i];
190: PetscCall(BVSetActiveColumns(ctx->Y,i*ctx->L+L_start,i*ctx->L+L_end));
191: PetscCall(BVGetMat(ctx->Y,&MC));
192: PetscCall(KSPMatSolve(ksp,V,MC));
193: PetscCall(BVRestoreMat(ctx->Y,&MC));
194: }
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*
199: Linear solves for the MULTISHIFT strategy
200: */
201: static PetscErrorCode EPSCISSSolve_MULTISHIFT(EPS eps,Mat V,PetscInt L_start,PetscInt L_end)
202: {
203: EPS_CISS *ctx = (EPS_CISS*)eps->data;
204: SlepcContourData contour;
205: PetscInt i,j,k=L_end-L_start;
206: Mat A_nest;
207: Vec b,b_nest,x,x_nest;
208: KSP ksp;
209: PC pc;
211: PetscFunctionBegin;
212: PetscCall(EPSCISSGetContour_Private(eps,&contour));
213: ksp = contour->ksp[0];
214: PetscCall(KSPGetPC(ksp,&pc));
215: PetscCall(PCGetOperators(pc,&A_nest,NULL));
216: PetscCall(MatCreateVecNestFromMultipleShifts(A_nest,NULL,&x_nest));
217: for (j=0;j<k;j++) { // TODO: solve all columns at once
218: PetscCall(MatDenseGetColumnVecRead(V,j,&b));
219: PetscCall(MatCreateVecNestFromMultipleShifts(A_nest,b,&b_nest));
220: PetscCall(KSPSolve(ksp,b_nest,x_nest));
221: for (i=0;i<contour->npoints;i++) {
222: PetscCall(VecNestGetSubVec(x_nest,i,&x));
223: PetscCall(BVInsertVec(ctx->Y,i*ctx->L+j,x));
224: }
225: PetscCall(VecDestroy(&b_nest));
226: PetscCall(MatDenseRestoreColumnVecRead(V,j,&b));
227: }
228: PetscCall(VecDestroy(&x_nest));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*
233: Y_i = (A-z_i B)^{-1}BV for every integration point
234: */
235: static PetscErrorCode EPSCISSSolve(EPS eps,Mat B,BV V,PetscInt L_start,PetscInt L_end)
236: {
237: EPS_CISS *ctx = (EPS_CISS*)eps->data;
238: Mat MV,BMV=NULL;
240: PetscFunctionBegin;
241: PetscCall(BVSetActiveColumns(V,L_start,L_end));
242: PetscCall(BVGetMat(V,&MV));
243: if (B) {
244: PetscCall(MatProductCreate(B,MV,NULL,&BMV));
245: PetscCall(MatProductSetType(BMV,MATPRODUCT_AB));
246: PetscCall(MatProductSetFromOptions(BMV));
247: PetscCall(MatProductSymbolic(BMV));
248: PetscCall(MatProductNumeric(BMV));
249: }
250: switch (ctx->strategy) {
251: case EPS_CISS_STRATEGY_USEST:
252: PetscCall(EPSCISSSolve_USEST(eps,B?BMV:MV,L_start,L_end));
253: break;
254: case EPS_CISS_STRATEGY_SPLIT:
255: PetscCall(EPSCISSSolve_SPLIT(eps,B?BMV:MV,L_start,L_end));
256: break;
257: case EPS_CISS_STRATEGY_MULTISHIFT:
258: PetscCall(EPSCISSSolve_MULTISHIFT(eps,B?BMV:MV,L_start,L_end));
259: break;
260: }
261: PetscCall(MatDestroy(&BMV));
262: PetscCall(BVRestoreMat(V,&MV));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
267: {
268: EPS_CISS *ctx = (EPS_CISS*)eps->data;
269: PetscInt i;
270: PetscScalar center;
271: PetscReal radius,a,b,c,d,rgscale;
272: #if PetscDefined(USE_COMPLEX)
273: PetscReal start_ang,end_ang,vscale,theta;
274: #endif
275: PetscBool isring,isellipse,isinterval;
277: PetscFunctionBegin;
278: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
279: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
280: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
281: PetscCall(RGGetScale(eps->rg,&rgscale));
282: if (isinterval) {
283: PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
284: if (c==d) {
285: for (i=0;i<nv;i++) {
286: #if PetscDefined(USE_COMPLEX)
287: eps->eigr[i] = PetscRealPart(eps->eigr[i]);
288: #else
289: eps->eigi[i] = 0;
290: #endif
291: }
292: }
293: }
294: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
295: if (isellipse) {
296: PetscCall(RGEllipseGetParameters(eps->rg,¢er,&radius,NULL));
297: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
298: } else if (isinterval) {
299: PetscCall(RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d));
300: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
301: for (i=0;i<nv;i++) {
302: if (c==d) eps->eigr[i] = ((eps->eigr[i]+1.0)*(b-a)/2.0+a)*rgscale;
303: if (a==b) {
304: #if PetscDefined(USE_COMPLEX)
305: eps->eigr[i] = ((eps->eigr[i]+1.0)*(d-c)/2.0+c)*rgscale*PETSC_i;
306: #else
307: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
308: #endif
309: }
310: }
311: } else {
312: center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
313: radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
314: for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
315: }
316: } else if (isring) { /* only supported in complex scalars */
317: #if PetscDefined(USE_COMPLEX)
318: PetscCall(RGRingGetParameters(eps->rg,¢er,&radius,&vscale,&start_ang,&end_ang,NULL));
319: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
320: for (i=0;i<nv;i++) {
321: theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
322: eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
323: }
324: } else {
325: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
326: }
327: #endif
328: }
329: }
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode EPSSetUp_CISS(EPS eps)
334: {
335: EPS_CISS *ctx = (EPS_CISS*)eps->data;
336: SlepcContourData contour;
337: PetscBool istrivial,isring,isellipse,isinterval,flg;
338: PetscReal c,d;
339: PetscInt nsplit;
340: PetscRandom rand;
341: PetscObjectId id;
342: PetscObjectState state;
343: Mat A[2],Psplit[2],T,J,Pa=NULL,Pb=NULL;
344: Vec v0;
346: PetscFunctionBegin;
347: EPSCheckNotStructured(eps);
348: if (eps->ncv==PETSC_DETERMINE) {
349: eps->ncv = ctx->L_max*ctx->M;
350: if (eps->ncv>eps->n) {
351: eps->ncv = eps->n;
352: ctx->L_max = eps->ncv/ctx->M;
353: PetscCheck(ctx->L_max,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
354: }
355: } else {
356: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
357: ctx->L_max = eps->ncv/ctx->M;
358: if (!ctx->L_max) {
359: ctx->L_max = 1;
360: eps->ncv = ctx->L_max*ctx->M;
361: }
362: }
363: ctx->L = PetscMin(ctx->L,ctx->L_max);
364: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 5;
365: if (eps->mpd==PETSC_DETERMINE) eps->mpd = eps->ncv;
366: if (!eps->which) eps->which = EPS_ALL;
367: PetscCheck(eps->which==EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
368: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
370: /* check region */
371: PetscCall(RGIsTrivial(eps->rg,&istrivial));
372: PetscCheck(!istrivial,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
373: PetscCall(RGGetComplement(eps->rg,&flg));
374: PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
375: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
376: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring));
377: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval));
378: PetscCheck(isellipse || isring || isinterval,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");
380: /* if the region has changed, then reset contour data */
381: PetscCall(PetscObjectGetId((PetscObject)eps->rg,&id));
382: PetscCall(PetscObjectStateGet((PetscObject)eps->rg,&state));
383: if (ctx->rgid && (id != ctx->rgid || state != ctx->rgstate)) {
384: PetscCall(SlepcContourDataDestroy(&ctx->contour));
385: PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of region\n"));
386: ctx->rgid = id; ctx->rgstate = state;
387: }
389: #if !PetscDefined(USE_COMPLEX)
390: PetscCheck(!isring,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
391: #endif
392: if (isinterval) {
393: PetscCall(RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d));
394: #if !PetscDefined(USE_COMPLEX)
395: PetscCheck(c==d && c==0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"In real scalars, endpoints of the imaginary axis must be both zero");
396: #endif
397: if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
398: }
399: if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;
401: if (!ctx->strategy) ctx->strategy = (ctx->npart>1)? EPS_CISS_STRATEGY_SPLIT: EPS_CISS_STRATEGY_USEST;
402: PetscCheck(ctx->strategy != EPS_CISS_STRATEGY_USEST || ctx->npart==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The EPS_CISS_STRATEGY_USEST strategy is not supported when partitions > 1");
404: PetscCall(EPSCISSGetContour_Private(eps,&contour));
406: if (eps->setfromoptionscalled && ctx->strategy != EPS_CISS_STRATEGY_USEST) PetscCall(PetscSubcommSetFromOptions(contour->subcomm));
408: PetscCall(EPSAllocateSolution(eps,0));
409: PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
410: if (ctx->weight) PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
411: PetscCall(PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma));
413: /* allocate basis vectors */
414: PetscCall(BVDestroy(&ctx->S));
415: PetscCall(BVDuplicateResize(eps->V,ctx->L*ctx->M,&ctx->S));
416: PetscCall(BVDestroy(&ctx->V));
417: PetscCall(BVDuplicateResize(eps->V,ctx->L,&ctx->V));
419: PetscCall(STGetMatrix(eps->st,0,&A[0]));
420: PetscCall(MatIsShell(A[0],&flg));
421: PetscCheck(!flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
422: if (eps->isgeneralized) PetscCall(STGetMatrix(eps->st,1,&A[1]));
423: else A[1] = NULL;
425: /* check if a user-defined split preconditioner has been set */
426: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
427: if (nsplit) {
428: PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Psplit[0]));
429: if (eps->isgeneralized) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Psplit[1]));
430: }
432: PetscCall(SlepcContourRedundantMat(contour,eps->isgeneralized?2:1,A,nsplit?Psplit:NULL));
433: if (contour->pA) {
434: PetscCall(BVGetColumn(ctx->V,0,&v0));
435: PetscCall(SlepcContourScatterCreate(contour,v0));
436: PetscCall(BVRestoreColumn(ctx->V,0,&v0));
437: PetscCall(BVDestroy(&ctx->pV));
438: PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->pV));
439: PetscCall(BVSetSizesFromVec(ctx->pV,contour->xsub,eps->n));
440: PetscCall(BVSetFromOptions(ctx->pV));
441: PetscCall(BVResize(ctx->pV,ctx->L,PETSC_FALSE));
442: }
444: EPSCheckDefinite(eps);
445: EPSCheckSinvertCondition(eps,ctx->strategy == EPS_CISS_STRATEGY_USEST," (with EPS_CISS_STRATEGY_USEST)");
447: PetscCall(BVDestroy(&ctx->Y));
448: if (contour->pA) {
449: PetscCall(BVCreate(PetscObjectComm((PetscObject)contour->xsub),&ctx->Y));
450: PetscCall(BVSetSizesFromVec(ctx->Y,contour->xsub,eps->n));
451: PetscCall(BVSetFromOptions(ctx->Y));
452: PetscCall(BVResize(ctx->Y,contour->npoints*ctx->L,PETSC_FALSE));
453: } else PetscCall(BVDuplicateResize(eps->V,contour->npoints*ctx->L,&ctx->Y));
455: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(DSSetType(eps->ds,DSGNHEP));
456: else if (eps->isgeneralized) {
457: if (eps->ishermitian && eps->ispositive) PetscCall(DSSetType(eps->ds,DSGHEP));
458: else PetscCall(DSSetType(eps->ds,DSGNHEP));
459: } else {
460: if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
461: else PetscCall(DSSetType(eps->ds,DSNHEP));
462: }
463: PetscCall(DSAllocate(eps->ds,eps->ncv));
465: #if !PetscDefined(USE_COMPLEX)
466: PetscCall(EPSSetWorkVecs(eps,3));
467: if (!eps->ishermitian) PetscCall(PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"));
468: #else
469: PetscCall(EPSSetWorkVecs(eps,2));
470: #endif
472: PetscCall(RGComputeQuadrature(eps->rg,ctx->quad==EPS_CISS_QUADRULE_CHEBYSHEV?RG_QUADRULE_CHEBYSHEV:RG_QUADRULE_TRAPEZOIDAL,ctx->N,ctx->omega,ctx->pp,ctx->weight));
473: J = (contour->pA && A[1])? contour->pA[1]: A[1];
474: if (ctx->strategy == EPS_CISS_STRATEGY_SPLIT || ctx->strategy == EPS_CISS_STRATEGY_MULTISHIFT) {
475: T = contour->pA? contour->pA[0]: A[0];
476: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
477: if (nsplit) {
478: if (contour->pA) {
479: Pa = contour->pP[0];
480: if (nsplit>1) Pb = contour->pP[1];
481: } else {
482: PetscCall(STGetSplitPreconditionerTerm(eps->st,0,&Pa));
483: if (nsplit>1) PetscCall(STGetSplitPreconditionerTerm(eps->st,1,&Pb));
484: }
485: }
486: PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
487: if (ctx->strategy == EPS_CISS_STRATEGY_SPLIT) PetscCall(EPSCISSSetUp_SPLIT(eps,T,J,Pa,Pb));
488: else PetscCall(EPSCISSSetUp_MULTISHIFT(eps,T,J,Pa,Pb));
489: }
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: static PetscErrorCode EPSSetUpSort_CISS(EPS eps)
494: {
495: SlepcSC sc;
497: PetscFunctionBegin;
498: /* fill sorting criterion context */
499: eps->sc->comparison = SlepcCompareSmallestReal;
500: eps->sc->comparisonctx = NULL;
501: eps->sc->map = NULL;
502: eps->sc->mapobj = NULL;
504: /* fill sorting criterion for DS */
505: PetscCall(DSGetSlepcSC(eps->ds,&sc));
506: sc->comparison = SlepcCompareLargestMagnitude;
507: sc->comparisonctx = NULL;
508: sc->map = NULL;
509: sc->mapobj = NULL;
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: static PetscErrorCode EPSSolve_CISS(EPS eps)
514: {
515: EPS_CISS *ctx = (EPS_CISS*)eps->data;
516: SlepcContourData contour;
517: Mat A,B=NULL,X,M,pA,pB,J;
518: BV V;
519: PetscInt i,j,ld,L_add=0,nv=0,L_base=ctx->L,inner,*inside;
520: PetscScalar *Mu,*H0,*H1=NULL,*rr,*temp;
521: PetscReal error,max_error,norm;
522: PetscBool *fl1;
523: Vec si,si1=NULL,w[3];
524: PetscRandom rand;
525: #if PetscDefined(USE_COMPLEX)
526: PetscBool isellipse;
527: PetscReal est_eig,eta;
528: #else
529: PetscReal normi;
530: #endif
532: PetscFunctionBegin;
533: w[0] = eps->work[0];
534: #if PetscDefined(USE_COMPLEX)
535: w[1] = NULL;
536: #else
537: w[1] = eps->work[2];
538: #endif
539: w[2] = eps->work[1];
540: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
542: PetscCall(STGetMatrix(eps->st,0,&A));
543: if (eps->isgeneralized) PetscCall(STGetMatrix(eps->st,1,&B));
544: PetscCall(EPSCISSGetContour_Private(eps,&contour));
545: J = (contour->pA && eps->isgeneralized)? contour->pA[1]: B;
546: V = contour->pA? ctx->pV: ctx->V;
548: PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
549: PetscCall(BVSetRandomSign(ctx->V));
550: PetscCall(BVGetRandomContext(ctx->V,&rand));
552: if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
553: PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
554: #if PetscDefined(USE_COMPLEX)
555: PetscCall(PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse));
556: if (isellipse) {
557: PetscCall(BVTraceQuadrature(ctx->Y,ctx->V,ctx->L,ctx->L,ctx->weight,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj,&est_eig));
558: PetscCall(PetscInfo(eps,"Estimated eigenvalue count: %f\n",(double)est_eig));
559: eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
560: L_add = PetscMax(0,(PetscInt)PetscCeilReal((est_eig*eta)/ctx->M)-ctx->L);
561: if (L_add>ctx->L_max-ctx->L) {
562: PetscCall(PetscInfo(eps,"Number of eigenvalues inside the contour path may be too large\n"));
563: L_add = ctx->L_max-ctx->L;
564: }
565: }
566: #endif
567: if (L_add>0) {
568: PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by Estimate #Eig\n",ctx->L,ctx->L+L_add));
569: PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
570: PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
571: PetscCall(BVSetRandomSign(ctx->V));
572: if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
573: ctx->L += L_add;
574: PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
575: }
576: PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
577: for (i=0;i<ctx->refine_blocksize;i++) {
578: PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
579: PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
580: PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
581: PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
582: PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
583: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
584: L_add = L_base;
585: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
586: PetscCall(PetscInfo(eps,"Changing L %" PetscInt_FMT " -> %" PetscInt_FMT " by SVD(H0)\n",ctx->L,ctx->L+L_add));
587: PetscCall(BVCISSResizeBases(ctx->S,contour->pA?ctx->pV:ctx->V,ctx->Y,ctx->L,ctx->L+L_add,ctx->M,contour->npoints));
588: PetscCall(BVSetActiveColumns(ctx->V,ctx->L,ctx->L+L_add));
589: PetscCall(BVSetRandomSign(ctx->V));
590: if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
591: ctx->L += L_add;
592: PetscCall(EPSCISSSolve(eps,J,V,ctx->L-L_add,ctx->L));
593: if (L_add) {
594: PetscCall(PetscFree2(Mu,H0));
595: PetscCall(PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0));
596: }
597: }
598: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1));
600: while (eps->reason == EPS_CONVERGED_ITERATING) {
601: eps->its++;
602: for (inner=0;inner<=ctx->refine_inner;inner++) {
603: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
604: PetscCall(BVDotQuadrature(ctx->Y,(contour->pA)?ctx->pV:ctx->V,Mu,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->subcomm,contour->npoints,ctx->useconj));
605: PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
606: PetscCall(PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0));
607: PetscCall(SlepcCISS_BH_SVD(H0,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv));
608: PetscCall(PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0));
609: break;
610: } else {
611: PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
612: PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
613: PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
614: PetscCall(BVCopy(ctx->S,ctx->V));
615: PetscCall(BVSVDAndRank(ctx->S,ctx->M,ctx->L,ctx->delta,BV_SVD_METHOD_REFINE,H0,ctx->sigma,&nv));
616: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
617: if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
618: PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
619: } else break;
620: }
621: }
622: eps->nconv = 0;
623: if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
624: else {
625: PetscCall(DSSetDimensions(eps->ds,nv,0,0));
626: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
628: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
629: PetscCall(CISS_BlockHankel(Mu,0,ctx->L,ctx->M,H0));
630: PetscCall(CISS_BlockHankel(Mu,1,ctx->L,ctx->M,H1));
631: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&temp));
632: for (j=0;j<nv;j++) {
633: for (i=0;i<nv;i++) {
634: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
635: }
636: }
637: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&temp));
638: PetscCall(DSGetArray(eps->ds,DS_MAT_B,&temp));
639: for (j=0;j<nv;j++) {
640: for (i=0;i<nv;i++) {
641: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
642: }
643: }
644: PetscCall(DSRestoreArray(eps->ds,DS_MAT_B,&temp));
645: } else {
646: PetscCall(BVSetActiveColumns(ctx->S,0,nv));
647: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&pA));
648: PetscCall(MatZeroEntries(pA));
649: PetscCall(BVMatProject(ctx->S,A,ctx->S,pA));
650: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&pA));
651: if (B) {
652: PetscCall(DSGetMat(eps->ds,DS_MAT_B,&pB));
653: PetscCall(MatZeroEntries(pB));
654: PetscCall(BVMatProject(ctx->S,B,ctx->S,pB));
655: PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&pB));
656: }
657: }
659: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
660: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
662: PetscCall(PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr));
663: PetscCall(rescale_eig(eps,nv));
664: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
665: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
666: PetscCall(SlepcCISS_isGhost(X,nv,ctx->sigma,ctx->spurious_threshold,fl1));
667: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
668: PetscCall(RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside));
669: for (i=0;i<nv;i++) {
670: if (fl1[i] && inside[i]>=0) {
671: rr[i] = 1.0;
672: eps->nconv++;
673: } else rr[i] = 0.0;
674: }
675: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv));
676: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
677: PetscCall(rescale_eig(eps,nv));
678: PetscCall(PetscFree3(fl1,inside,rr));
679: PetscCall(BVSetActiveColumns(eps->V,0,nv));
680: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
681: PetscCall(BVSumQuadrature(ctx->S,ctx->Y,ctx->M,ctx->L,ctx->L,ctx->weight,ctx->pp,contour->scatterin,contour->subcomm,contour->npoints,ctx->useconj));
682: PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
683: PetscCall(BVCopy(ctx->S,ctx->V));
684: PetscCall(BVSetActiveColumns(ctx->S,0,nv));
685: }
686: PetscCall(BVCopy(ctx->S,eps->V));
688: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
689: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
690: PetscCall(BVMultInPlace(ctx->S,X,0,eps->nconv));
691: if (eps->ishermitian) PetscCall(BVMultInPlace(eps->V,X,0,eps->nconv));
692: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
693: max_error = 0.0;
694: for (i=0;i<eps->nconv;i++) {
695: PetscCall(BVGetColumn(ctx->S,i,&si));
696: #if !PetscDefined(USE_COMPLEX)
697: if (eps->eigi[i]!=0.0) PetscCall(BVGetColumn(ctx->S,i+1,&si1));
698: #endif
699: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,si1,w,&error));
700: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) { /* vector is not normalized */
701: PetscCall(VecNorm(si,NORM_2,&norm));
702: #if !PetscDefined(USE_COMPLEX)
703: if (eps->eigi[i]!=0.0) {
704: PetscCall(VecNorm(si1,NORM_2,&normi));
705: norm = SlepcAbsEigenvalue(norm,normi);
706: }
707: #endif
708: error /= norm;
709: }
710: PetscCall((*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx));
711: PetscCall(BVRestoreColumn(ctx->S,i,&si));
712: #if !PetscDefined(USE_COMPLEX)
713: if (eps->eigi[i]!=0.0) {
714: PetscCall(BVRestoreColumn(ctx->S,i+1,&si1));
715: i++;
716: }
717: #endif
718: max_error = PetscMax(max_error,error);
719: }
721: if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
722: else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
723: else {
724: if (eps->nconv > ctx->L) nv = eps->nconv;
725: else if (ctx->L > nv) nv = ctx->L;
726: nv = PetscMin(nv,ctx->L*ctx->M);
727: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nv,ctx->L,NULL,&M));
728: PetscCall(MatSetRandom(M,rand));
729: PetscCall(BVSetActiveColumns(ctx->S,0,nv));
730: PetscCall(BVMultInPlace(ctx->S,M,0,ctx->L));
731: PetscCall(MatDestroy(&M));
732: PetscCall(BVSetActiveColumns(ctx->S,0,ctx->L));
733: PetscCall(BVSetActiveColumns(ctx->V,0,ctx->L));
734: PetscCall(BVCopy(ctx->S,ctx->V));
735: if (contour->pA) PetscCall(BVScatter(ctx->V,ctx->pV,contour->scatterin,contour->xdup));
736: PetscCall(EPSCISSSolve(eps,J,V,0,ctx->L));
737: }
738: }
739: }
740: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(PetscFree(H1));
741: PetscCall(PetscFree2(Mu,H0));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode EPSComputeVectors_CISS(EPS eps)
746: {
747: EPS_CISS *ctx = (EPS_CISS*)eps->data;
748: PetscInt n;
749: Mat Z,B=NULL;
751: PetscFunctionBegin;
752: if (eps->ishermitian) {
753: if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
754: else PetscCall(EPSComputeVectors_Hermitian(eps));
755: if (eps->isgeneralized && eps->ispositive && ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
756: /* normalize to have unit B-norm */
757: PetscCall(STGetMatrix(eps->st,1,&B));
758: PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
759: PetscCall(BVNormalize(eps->V,NULL));
760: PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
761: }
762: PetscFunctionReturn(PETSC_SUCCESS);
763: }
764: PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
765: PetscCall(BVSetActiveColumns(eps->V,0,n));
767: /* right eigenvectors */
768: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
770: /* V = V * Z */
771: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
772: PetscCall(BVMultInPlace(eps->V,Z,0,n));
773: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
774: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
776: /* normalize */
777: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) PetscCall(BVNormalize(eps->V,NULL));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
782: {
783: EPS_CISS *ctx = (EPS_CISS*)eps->data;
784: PetscInt oN,oL,oM,oLmax,onpart;
785: PetscMPIInt size;
787: PetscFunctionBegin;
788: oN = ctx->N;
789: if (ip == PETSC_DETERMINE) {
790: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
791: } else if (ip != PETSC_CURRENT) {
792: PetscCheck(ip>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
793: PetscCheck(ip%2==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
794: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
795: }
796: oL = ctx->L;
797: if (bs == PETSC_DETERMINE) {
798: ctx->L = 16;
799: } else if (bs != PETSC_CURRENT) {
800: PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
801: ctx->L = bs;
802: }
803: oM = ctx->M;
804: if (ms == PETSC_DETERMINE) {
805: ctx->M = ctx->N/4;
806: } else if (ms != PETSC_CURRENT) {
807: PetscCheck(ms>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
808: PetscCheck(ms<=ctx->N,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
809: ctx->M = ms;
810: }
811: onpart = ctx->npart;
812: if (npart == PETSC_DETERMINE) {
813: ctx->npart = 1;
814: } else if (npart != PETSC_CURRENT) {
815: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
816: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
817: ctx->npart = npart;
818: }
819: oLmax = ctx->L_max;
820: if (bsmax == PETSC_DETERMINE) {
821: ctx->L_max = 64;
822: } else if (bsmax != PETSC_CURRENT) {
823: PetscCheck(bsmax>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
824: ctx->L_max = PetscMax(bsmax,ctx->L);
825: }
826: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) {
827: PetscCall(SlepcContourDataDestroy(&ctx->contour));
828: PetscCall(PetscInfo(eps,"Resetting the contour data structure due to a change of parameters\n"));
829: eps->state = EPS_STATE_INITIAL;
830: }
831: ctx->isreal = realmats;
832: if (oL != ctx->L || oM != ctx->M || oLmax != ctx->L_max) eps->state = EPS_STATE_INITIAL;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /*@
837: EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.
839: Logically Collective
841: Input Parameters:
842: + eps - the linear eigensolver context
843: . ip - number of integration points
844: . bs - block size
845: . ms - moment size
846: . npart - number of partitions when splitting the communicator
847: . bsmax - max block size
848: - realmats - $A$ and $B$ are real
850: Options Database Keys:
851: + -eps_ciss_integration_points ip - sets the number of integration points
852: . -eps_ciss_blocksize bs - sets the block size
853: . -eps_ciss_moments ms - sets the moment size
854: . -eps_ciss_partitions npart - sets the number of partitions
855: . -eps_ciss_maxblocksize bsmax - sets the maximum block size
856: - -eps_ciss_realmats (true|false) - $A$ and $B$ are real
858: Notes:
859: For all integer arguments, you can use `PETSC_CURRENT` to keep the current value, and
860: `PETSC_DETERMINE` to set them to a default value.
862: The default number of partitions is 1. This means the internal `KSP` object is shared
863: among all processes of the `EPS` communicator. Otherwise, the communicator is split
864: into `npart` communicators, so that `npart` `KSP` solves proceed simultaneously.
866: For a detailed description of the parameters see {cite:p}`Mae16`.
868: Level: advanced
870: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetSizes()`, `EPSCISSSetThreshold()`, `EPSCISSSetRefinement()`
871: @*/
872: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
873: {
874: PetscFunctionBegin;
882: PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
886: static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
887: {
888: EPS_CISS *ctx = (EPS_CISS*)eps->data;
890: PetscFunctionBegin;
891: if (ip) *ip = ctx->N;
892: if (bs) *bs = ctx->L;
893: if (ms) *ms = ctx->M;
894: if (npart) *npart = ctx->npart;
895: if (bsmax) *bsmax = ctx->L_max;
896: if (realmats) *realmats = ctx->isreal;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.
903: Not Collective
905: Input Parameter:
906: . eps - the linear eigensolver context
908: Output Parameters:
909: + ip - number of integration points
910: . bs - block size
911: . ms - moment size
912: . npart - number of partitions when splitting the communicator
913: . bsmax - max block size
914: - realmats - $A$ and $B$ are real
916: Level: advanced
918: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetSizes()`
919: @*/
920: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
921: {
922: PetscFunctionBegin;
924: PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
925: PetscFunctionReturn(PETSC_SUCCESS);
926: }
928: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
929: {
930: EPS_CISS *ctx = (EPS_CISS*)eps->data;
932: PetscFunctionBegin;
933: if (delta == (PetscReal)PETSC_DETERMINE) {
934: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
935: } else if (delta != (PetscReal)PETSC_CURRENT) {
936: PetscCheck(delta>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
937: ctx->delta = delta;
938: }
939: if (spur == (PetscReal)PETSC_DETERMINE) {
940: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
941: } else if (spur != (PetscReal)PETSC_CURRENT) {
942: PetscCheck(spur>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
943: ctx->spurious_threshold = spur;
944: }
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: /*@
949: EPSCISSSetThreshold - Sets the values of various threshold parameters in
950: the CISS solver.
952: Logically Collective
954: Input Parameters:
955: + eps - the linear eigensolver context
956: . delta - threshold for numerical rank
957: - spur - spurious threshold (to discard spurious eigenpairs)
959: Options Database Keys:
960: + -eps_ciss_delta delta - sets the delta
961: - -eps_ciss_spurious_threshold spur - sets the spurious threshold
963: Notes:
964: `PETSC_CURRENT` can be used to preserve the current value of any of the
965: arguments, and `PETSC_DETERMINE` to set them to a default value.
967: For a detailed description of the parameters see {cite:p}`Mae16`.
969: Level: advanced
971: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetThreshold()`, `EPSCISSSetSizes()`, `EPSCISSSetRefinement()`
972: @*/
973: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
974: {
975: PetscFunctionBegin;
979: PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
984: {
985: EPS_CISS *ctx = (EPS_CISS*)eps->data;
987: PetscFunctionBegin;
988: if (delta) *delta = ctx->delta;
989: if (spur) *spur = ctx->spurious_threshold;
990: PetscFunctionReturn(PETSC_SUCCESS);
991: }
993: /*@
994: EPSCISSGetThreshold - Gets the values of various threshold parameters
995: in the CISS solver.
997: Not Collective
999: Input Parameter:
1000: . eps - the linear eigensolver context
1002: Output Parameters:
1003: + delta - threshold for numerical rank
1004: - spur - spurious threshold (to discard spurious eigenpairs)
1006: Level: advanced
1008: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetThreshold()`
1009: @*/
1010: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
1011: {
1012: PetscFunctionBegin;
1014: PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
1019: {
1020: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1022: PetscFunctionBegin;
1023: if (inner == PETSC_DETERMINE) {
1024: ctx->refine_inner = 0;
1025: } else if (inner != PETSC_CURRENT) {
1026: PetscCheck(inner>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
1027: ctx->refine_inner = inner;
1028: }
1029: if (blsize == PETSC_DETERMINE) {
1030: ctx->refine_blocksize = 0;
1031: } else if (blsize != PETSC_CURRENT) {
1032: PetscCheck(blsize>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
1033: ctx->refine_blocksize = blsize;
1034: }
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: EPSCISSSetRefinement - Sets the values of various refinement parameters
1040: in the CISS solver.
1042: Logically Collective
1044: Input Parameters:
1045: + eps - the linear eigensolver context
1046: . inner - number of iterative refinement iterations (inner loop)
1047: - blsize - number of iterative refinement iterations (blocksize loop)
1049: Options Database Keys:
1050: + -eps_ciss_refine_inner inner - sets number of inner iterations
1051: - -eps_ciss_refine_blocksize blsize - sets number of blocksize iterations
1053: Notes:
1054: `PETSC_CURRENT` can be used to preserve the current value of any of the
1055: arguments, and `PETSC_DETERMINE` to set them to a default of 0 (no refinement).
1057: For a detailed description of the parameters see {cite:p}`Mae16`.
1059: Level: advanced
1061: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetRefinement()`, `EPSCISSSetSizes()`, `EPSCISSSetThreshold()`
1062: @*/
1063: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
1064: {
1065: PetscFunctionBegin;
1069: PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
1074: {
1075: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1077: PetscFunctionBegin;
1078: if (inner) *inner = ctx->refine_inner;
1079: if (blsize) *blsize = ctx->refine_blocksize;
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: EPSCISSGetRefinement - Gets the values of various refinement parameters
1085: in the CISS solver.
1087: Not Collective
1089: Input Parameter:
1090: . eps - the linear eigensolver context
1092: Output Parameters:
1093: + inner - number of iterative refinement iterations (inner loop)
1094: - blsize - number of iterative refinement iterations (blocksize loop)
1096: Level: advanced
1098: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetRefinement()`
1099: @*/
1100: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
1101: {
1102: PetscFunctionBegin;
1104: PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: static PetscErrorCode EPSCISSSetStrategy_CISS(EPS eps,EPSCISSStrategy strategy)
1109: {
1110: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1112: PetscFunctionBegin;
1113: if (ctx->strategy != strategy) {
1114: ctx->strategy = strategy;
1115: eps->state = EPS_STATE_INITIAL;
1116: }
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*@
1121: EPSCISSSetStrategy - Sets the strategy to be used when performing linear solves
1122: associated with integration points in the CISS solver.
1124: Logically Collective
1126: Input Parameters:
1127: + eps - the linear eigensolver context
1128: - strategy - the strategy, see `EPSCISSStrategy` for possible values
1130: Options Database Key:
1131: . -eps_ciss_strategy (usest|split|multishift) - sets the strategy
1133: Notes:
1134: When the `usest` strategy is selected the linear solves can be configured by
1135: setting options for the `KSP` object obtained with `STGetKSP()`.
1136: Otherwise, several `KSP` objects are created, which can be accessed
1137: with `EPSCISSGetKSPs()`.
1139: The default is to use the `ST`, unless several partitions have been
1140: specified, see `EPSCISSSetSizes()`.
1142: Level: advanced
1144: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetStrategy()`, `EPSCISSSetSizes()`, `EPSCISSGetKSPs()`, `STGetKSP()`
1145: @*/
1146: PetscErrorCode EPSCISSSetStrategy(EPS eps,EPSCISSStrategy strategy)
1147: {
1148: PetscFunctionBegin;
1151: PetscTryMethod(eps,"EPSCISSSetStrategy_C",(EPS,EPSCISSStrategy),(eps,strategy));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: static PetscErrorCode EPSCISSGetStrategy_CISS(EPS eps,EPSCISSStrategy *strategy)
1156: {
1157: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1159: PetscFunctionBegin;
1160: *strategy = ctx->strategy;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: /*@
1165: EPSCISSGetStrategy - Gets the strategy to be used when performing linear solves
1166: associated with integration points in the CISS solver.
1168: Not Collective
1170: Input Parameter:
1171: . eps - the linear eigensolver context
1173: Output Parameter:
1174: . strategy - the strategy
1176: Level: advanced
1178: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetStrategy()`, `EPSCISSStrategy`
1179: @*/
1180: PetscErrorCode EPSCISSGetStrategy(EPS eps,EPSCISSStrategy *strategy)
1181: {
1182: PetscFunctionBegin;
1184: PetscAssertPointer(strategy,2);
1185: PetscUseMethod(eps,"EPSCISSGetStrategy_C",(EPS,EPSCISSStrategy*),(eps,strategy));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1190: {
1191: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1193: PetscFunctionBegin;
1194: if (ctx->quad != quad) {
1195: ctx->quad = quad;
1196: eps->state = EPS_STATE_INITIAL;
1197: }
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: /*@
1202: EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.
1204: Logically Collective
1206: Input Parameters:
1207: + eps - the linear eigensolver context
1208: - quad - the quadrature rule, see `EPSCISSQuadRule` for possible values
1210: Options Database Key:
1211: . -eps_ciss_quadrule (trapezoidal|chebyshev) - sets the quadrature rule
1213: Notes:
1214: By default, the trapezoidal rule is used (`EPS_CISS_QUADRULE_TRAPEZOIDAL`).
1216: If the `chebyshev` option is specified (`EPS_CISS_QUADRULE_CHEBYSHEV`), then
1217: Chebyshev points are used as quadrature points.
1219: Level: advanced
1221: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetQuadRule()`, `EPSCISSQuadRule`
1222: @*/
1223: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1224: {
1225: PetscFunctionBegin;
1228: PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1233: {
1234: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1236: PetscFunctionBegin;
1237: *quad = ctx->quad;
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*@
1242: EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.
1244: Not Collective
1246: Input Parameter:
1247: . eps - the linear eigensolver context
1249: Output Parameter:
1250: . quad - quadrature rule
1252: Level: advanced
1254: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetQuadRule()`, `EPSCISSQuadRule`
1255: @*/
1256: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1257: {
1258: PetscFunctionBegin;
1260: PetscAssertPointer(quad,2);
1261: PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }
1265: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1266: {
1267: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1269: PetscFunctionBegin;
1270: if (ctx->extraction != extraction) {
1271: ctx->extraction = extraction;
1272: eps->state = EPS_STATE_INITIAL;
1273: }
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1277: /*@
1278: EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.
1280: Logically Collective
1282: Input Parameters:
1283: + eps - the linear eigensolver context
1284: - extraction - the extraction technique, see `EPSCISSExtraction` for possible values
1286: Options Database Key:
1287: . -eps_ciss_extraction (ritz|hankel) - sets the extraction technique
1289: Notes:
1290: By default, the Rayleigh-Ritz extraction is used (`EPS_CISS_EXTRACTION_RITZ`),
1291: see {cite:p}`Sak07`.
1293: If the `hankel` option is specified (`EPS_CISS_EXTRACTION_HANKEL`), then
1294: the block Hankel method is used for extracting eigenpairs {cite:p}`Sak03`.
1296: Level: advanced
1298: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSGetExtraction()`, `EPSCISSExtraction`
1299: @*/
1300: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1301: {
1302: PetscFunctionBegin;
1305: PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1306: PetscFunctionReturn(PETSC_SUCCESS);
1307: }
1309: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1310: {
1311: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1313: PetscFunctionBegin;
1314: *extraction = ctx->extraction;
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.
1321: Not Collective
1323: Input Parameter:
1324: . eps - the linear eigensolver context
1326: Output Parameter:
1327: . extraction - extraction technique
1329: Level: advanced
1331: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetExtraction()`, `EPSCISSExtraction`
1332: @*/
1333: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1334: {
1335: PetscFunctionBegin;
1337: PetscAssertPointer(extraction,2);
1338: PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1343: {
1344: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1345: SlepcContourData contour;
1346: PetscInt i,nsplit;
1347: KSP ksps,kspm;
1348: PC pc;
1349: MPI_Comm child;
1351: PetscFunctionBegin;
1352: PetscCall(EPSCISSGetContour_Private(eps,&contour));
1353: if (!contour->ksp) {
1354: switch (ctx->strategy) {
1355: case EPS_CISS_STRATEGY_USEST:
1356: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Should not call EPSCISSGetKSPs() with EPS_CISS_STRATEGY_USEST");
1357: break;
1358: case EPS_CISS_STRATEGY_SPLIT:
1359: contour->nksp = contour->npoints;
1360: PetscCall(PetscMalloc1(contour->nksp,&contour->ksp));
1361: PetscCall(EPSGetST(eps,&eps->st));
1362: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
1363: PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
1364: for (i=0;i<contour->nksp;i++) {
1365: PetscCall(KSPCreate(child,&contour->ksp[i]));
1366: PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[i],(PetscObject)eps,1));
1367: PetscCall(KSPSetOptionsPrefix(contour->ksp[i],((PetscObject)eps)->prefix));
1368: PetscCall(KSPAppendOptionsPrefix(contour->ksp[i],"eps_ciss_"));
1369: PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[i],((PetscObject)eps)->options));
1370: PetscCall(KSPSetErrorIfNotConverged(contour->ksp[i],PETSC_TRUE));
1371: PetscCall(KSPSetTolerances(contour->ksp[i],SlepcDefaultTol(eps->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1372: PetscCall(KSPGetPC(contour->ksp[i],&pc));
1373: if (nsplit) {
1374: PetscCall(KSPSetType(contour->ksp[i],KSPBCGS));
1375: PetscCall(PCSetType(pc,PCBJACOBI));
1376: } else {
1377: PetscCall(KSPSetType(contour->ksp[i],KSPPREONLY));
1378: PetscCall(PCSetType(pc,PCLU));
1379: }
1380: }
1381: break;
1382: case EPS_CISS_STRATEGY_MULTISHIFT:
1383: contour->nksp = 1; /* one solver per subcomm */
1384: PetscCall(PetscMalloc1(contour->nksp,&contour->ksp));
1385: PetscCall(EPSGetST(eps,&eps->st));
1386: PetscCall(STGetSplitPreconditionerInfo(eps->st,&nsplit,NULL));
1387: PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
1388: PetscCall(KSPCreate(child,&contour->ksp[0]));
1389: PetscCall(PetscObjectIncrementTabLevel((PetscObject)contour->ksp[0],(PetscObject)eps,1));
1390: PetscCall(KSPSetOptionsPrefix(contour->ksp[0],((PetscObject)eps)->prefix));
1391: PetscCall(KSPAppendOptionsPrefix(contour->ksp[0],"eps_ciss_"));
1392: PetscCall(PetscObjectSetOptions((PetscObject)contour->ksp[0],((PetscObject)eps)->options));
1393: PetscCall(KSPSetErrorIfNotConverged(contour->ksp[0],PETSC_TRUE));
1394: PetscCall(KSPSetTolerances(contour->ksp[0],SlepcDefaultTol(eps->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1395: PetscCall(KSPGetPC(contour->ksp[0],&pc));
1396: PetscCall(KSPSetType(contour->ksp[0],KSPEKSM));
1397: PetscCall(PCSetType(pc,PCNONE));
1398: PetscCall(KSPEKSMGetKSP(contour->ksp[0],&ksps,eps->isgeneralized?&kspm:NULL));
1399: PetscCall(KSPGetPC(ksps,&pc));
1400: PetscCall(KSPSetType(ksps,KSPPREONLY));
1401: PetscCall(PCSetType(pc,PCLU));
1402: if (eps->isgeneralized) {
1403: PetscCall(KSPGetPC(kspm,&pc));
1404: PetscCall(KSPSetType(kspm,KSPPREONLY));
1405: PetscCall(PCSetType(pc,PCLU));
1406: }
1407: break;
1408: }
1409: }
1410: if (nsolve) *nsolve = contour->nksp;
1411: if (ksp) *ksp = contour->ksp;
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: /*@
1416: EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1417: the CISS solver.
1419: Not Collective
1421: Input Parameter:
1422: . eps - the linear eigensolver context
1424: Output Parameters:
1425: + nsolve - number of solver objects
1426: - ksp - array of linear solver object
1428: Notes:
1429: In `EPS_CISS_STRATEGY_SPLIT` the number of `KSP` solvers is equal to the number of
1430: integration points divided by the number of partitions, see `EPSCISSSetSizes()`. This
1431: value is halved in the case of real matrices with a region centered at the real axis.
1433: In `EPS_CISS_STRATEGY_MULTISHIFT` only one `KSP` solver is returned.
1435: When the number of partitions is larger than one, MPI processes belonging to different
1436: subcommunicators will obtain different `KSP` objects.
1438: Level: advanced
1440: .seealso: [](ch:eps), `EPSCISS`, `EPSCISSSetSizes()`, `EPSCISSSetStrategy()`
1441: @*/
1442: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1443: {
1444: PetscFunctionBegin;
1446: PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: static PetscErrorCode EPSReset_CISS(EPS eps)
1451: {
1452: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1454: PetscFunctionBegin;
1455: PetscCall(BVDestroy(&ctx->S));
1456: PetscCall(BVDestroy(&ctx->V));
1457: PetscCall(BVDestroy(&ctx->Y));
1458: if (ctx->strategy != EPS_CISS_STRATEGY_USEST) PetscCall(SlepcContourDataReset(ctx->contour));
1459: PetscCall(BVDestroy(&ctx->pV));
1460: PetscFunctionReturn(PETSC_SUCCESS);
1461: }
1463: static PetscErrorCode EPSSetFromOptions_CISS(EPS eps,PetscOptionItems PetscOptionsObject)
1464: {
1465: PetscReal r3,r4;
1466: PetscInt i1,i2,i3,i4,i5,i6,i7;
1467: PetscBool b1,flg,flg2,flg3,flg4,flg5,flg6;
1468: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1469: EPSCISSQuadRule quad;
1470: EPSCISSExtraction extraction;
1471: EPSCISSStrategy strategy;
1473: PetscFunctionBegin;
1474: PetscOptionsHeadBegin(PetscOptionsObject,"EPS CISS Options");
1476: PetscCall(EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1));
1477: PetscCall(PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,&flg));
1478: PetscCall(PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,&flg2));
1479: PetscCall(PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,&flg3));
1480: PetscCall(PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,&flg4));
1481: PetscCall(PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,&flg5));
1482: PetscCall(PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,&flg6));
1483: if (flg || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1));
1485: PetscCall(EPSCISSGetThreshold(eps,&r3,&r4));
1486: PetscCall(PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,&flg));
1487: PetscCall(PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,&flg2));
1488: if (flg || flg2) PetscCall(EPSCISSSetThreshold(eps,r3,r4));
1490: PetscCall(EPSCISSGetRefinement(eps,&i6,&i7));
1491: PetscCall(PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,&flg));
1492: PetscCall(PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,&flg2));
1493: if (flg || flg2) PetscCall(EPSCISSSetRefinement(eps,i6,i7));
1495: PetscCall(PetscOptionsDeprecated("-eps_ciss_usest", NULL, "3.26", "Use -eps_ciss_strategy usest"));
1496: PetscCall(PetscOptionsEnum("-eps_ciss_strategy","Strategy for linear solves","EPSCISSSetStrategy",EPSCISSStrategies,(PetscEnum)ctx->strategy,(PetscEnum*)&strategy,&flg));
1497: if (flg) PetscCall(EPSCISSSetStrategy(eps,strategy));
1499: PetscCall(PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg));
1500: if (flg) PetscCall(EPSCISSSetQuadRule(eps,quad));
1502: PetscCall(PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg));
1503: if (flg) PetscCall(EPSCISSSetExtraction(eps,extraction));
1505: PetscOptionsHeadEnd();
1506: PetscFunctionReturn(PETSC_SUCCESS);
1507: }
1509: static PetscErrorCode EPSDestroy_CISS(EPS eps)
1510: {
1511: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1513: PetscFunctionBegin;
1514: PetscCall(SlepcContourDataDestroy(&ctx->contour));
1515: PetscCall(PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma));
1516: PetscCall(PetscFree(eps->data));
1517: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL));
1518: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL));
1519: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL));
1520: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL));
1521: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL));
1522: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL));
1523: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetStrategy_C",NULL));
1524: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetStrategy_C",NULL));
1525: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL));
1526: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL));
1527: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL));
1528: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL));
1529: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL));
1530: PetscFunctionReturn(PETSC_SUCCESS);
1531: }
1533: static PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1534: {
1535: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1536: PetscBool isascii;
1537: PetscViewer sviewer;
1539: PetscFunctionBegin;
1540: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1541: if (isascii) {
1542: PetscCall(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));
1543: if (ctx->isreal) PetscCall(PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n"));
1544: PetscCall(PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold));
1545: PetscCall(PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %" PetscInt_FMT ", blocksize: %" PetscInt_FMT " }\n",ctx->refine_inner, ctx->refine_blocksize));
1546: PetscCall(PetscViewerASCIIPrintf(viewer," extraction: %s\n",EPSCISSExtractions[ctx->extraction]));
1547: PetscCall(PetscViewerASCIIPrintf(viewer," quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]));
1548: switch (ctx->strategy) {
1549: case EPS_CISS_STRATEGY_USEST:
1550: PetscCall(PetscViewerASCIIPrintf(viewer," using ST for linear solves\n"));
1551: break;
1552: case EPS_CISS_STRATEGY_SPLIT:
1553: case EPS_CISS_STRATEGY_MULTISHIFT:
1554: PetscCall(EPSCISSGetKSPs(eps,NULL,NULL));
1555: PetscCall(PetscViewerASCIIPushTab(viewer));
1556: if (ctx->npart>1 && ctx->contour->subcomm) {
1557: PetscCall(PetscViewerGetSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1558: if (!ctx->contour->subcomm->color) PetscCall(KSPView(ctx->contour->ksp[0],sviewer));
1559: PetscCall(PetscViewerFlush(sviewer));
1560: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->contour->subcomm->child,&sviewer));
1561: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1562: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1563: } else PetscCall(KSPView(ctx->contour->ksp[0],viewer));
1564: PetscCall(PetscViewerASCIIPopTab(viewer));
1565: break;
1566: }
1567: }
1568: PetscFunctionReturn(PETSC_SUCCESS);
1569: }
1571: static PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1572: {
1573: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1574: PetscBool usest;
1575: KSP ksp;
1576: PC pc;
1578: PetscFunctionBegin;
1579: if (!((PetscObject)eps->st)->type_name) {
1580: if (!ctx->strategy) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1581: else usest = (ctx->strategy==EPS_CISS_STRATEGY_USEST)? PETSC_TRUE: PETSC_FALSE;
1582: if (usest) PetscCall(STSetType(eps->st,STSINVERT));
1583: else {
1584: /* we are not going to use ST, so avoid factorizing the matrix */
1585: PetscCall(STSetType(eps->st,STSHIFT));
1586: if (eps->isgeneralized) {
1587: PetscCall(STGetKSP(eps->st,&ksp));
1588: PetscCall(KSPGetPC(ksp,&pc));
1589: PetscCall(PCSetType(pc,PCNONE));
1590: }
1591: }
1592: }
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: /*MC
1597: EPSCISS - EPSCISS = "ciss" - A contour integral eigensolver based on the
1598: Sakurai-Sugiura scheme.
1600: Notes:
1601: This solver is based on the numerical contour integration idea
1602: proposed initially by {cite:t}`Sak03` and improved later by adding
1603: a Rayleigh-Ritz projection step {cite:p}`Sak07`.
1605: Contour integral methods are able to compute all eigenvalues
1606: lying inside a region of the complex plane. Use `EPSGetRG()` to
1607: specify the region. However, the computational cost is usually high
1608: because multiple linear systems must be solved. For this, we can
1609: use the `KSP` object inside `ST`, or several independent `KSP`s,
1610: see `EPSCISSSetStrategy()`.
1612: Details of the implementation in SLEPc can be found in {cite:p}`Mae16`.
1614: Level: beginner
1616: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSGetRG()`
1617: M*/
1618: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1619: {
1620: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1622: PetscFunctionBegin;
1623: PetscCall(PetscNew(&ctx));
1624: eps->data = ctx;
1626: eps->useds = PETSC_TRUE;
1627: eps->categ = EPS_CATEGORY_CONTOUR;
1629: eps->ops->solve = EPSSolve_CISS;
1630: eps->ops->setup = EPSSetUp_CISS;
1631: eps->ops->setupsort = EPSSetUpSort_CISS;
1632: eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1633: eps->ops->destroy = EPSDestroy_CISS;
1634: eps->ops->reset = EPSReset_CISS;
1635: eps->ops->view = EPSView_CISS;
1636: eps->ops->computevectors = EPSComputeVectors_CISS;
1637: eps->ops->setdefaultst = EPSSetDefaultST_CISS;
1639: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS));
1640: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS));
1641: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS));
1642: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS));
1643: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS));
1644: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS));
1645: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetStrategy_C",EPSCISSSetStrategy_CISS));
1646: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetStrategy_C",EPSCISSGetStrategy_CISS));
1647: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS));
1648: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS));
1649: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS));
1650: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS));
1651: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS));
1653: /* set default values of parameters */
1654: ctx->N = 32;
1655: ctx->L = 16;
1656: ctx->M = ctx->N/4;
1657: ctx->delta = SLEPC_DEFAULT_TOL*1e-4;
1658: ctx->L_max = 64;
1659: ctx->spurious_threshold = PetscSqrtReal(SLEPC_DEFAULT_TOL);
1660: ctx->isreal = PETSC_FALSE;
1661: ctx->refine_inner = 0;
1662: ctx->refine_blocksize = 0;
1663: ctx->npart = 1;
1664: ctx->quad = (EPSCISSQuadRule)0;
1665: ctx->extraction = EPS_CISS_EXTRACTION_RITZ;
1666: PetscFunctionReturn(PETSC_SUCCESS);
1667: }