Actual source code: ciss.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/epsimpl.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; /* A and B are real */
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 pV;
58: BV Y;
59: Vec xsub;
60: Vec xdup;
61: KSP *ksp; /* ksp array for storing factorizations at integration points */
62: PetscBool useconj;
63: PetscReal est_eig;
64: VecScatter scatterin;
65: Mat pA,pB;
66: PetscSubcomm subcomm;
67: PetscBool usest;
68: PetscBool usest_set; /* whether the user set the usest flag or not */
69: EPSCISSQuadRule quad;
70: EPSCISSExtraction extraction;
71: } EPS_CISS;
73: /* destroy KSP objects when the number of solve points changes */
74: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSolvers(EPS eps)
75: {
77: PetscInt i;
78: EPS_CISS *ctx = (EPS_CISS*)eps->data;
81: if (ctx->ksp) {
82: for (i=0;i<ctx->num_solve_point;i++) {
83: KSPDestroy(&ctx->ksp[i]);
84: }
85: PetscFree(ctx->ksp);
86: }
87: return(0);
88: }
90: /* clean PetscSubcomm object when the number of partitions changes */
91: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSubcomm(EPS eps)
92: {
94: EPS_CISS *ctx = (EPS_CISS*)eps->data;
97: EPSCISSResetSolvers(eps);
98: PetscSubcommDestroy(&ctx->subcomm);
99: return(0);
100: }
102: /* determine whether half of integration points can be avoided (use its conjugates);
103: depends on isreal and the center of the region */
104: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUseConj(EPS eps,PetscBool *useconj)
105: {
107: PetscScalar center;
108: PetscReal c,d;
109: PetscBool isellipse,isinterval;
110: #if defined(PETSC_USE_COMPLEX)
111: EPS_CISS *ctx = (EPS_CISS*)eps->data;
112: #endif
115: *useconj = PETSC_FALSE;
116: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
117: if (isellipse) {
118: RGEllipseGetParameters(eps->rg,¢er,NULL,NULL);
119: #if defined(PETSC_USE_COMPLEX)
120: *useconj = (ctx->isreal && PetscImaginaryPart(center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
121: #endif
122: } else {
123: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
124: if (isinterval) {
125: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
126: #if defined(PETSC_USE_COMPLEX)
127: *useconj = (ctx->isreal && c==d)? PETSC_TRUE: PETSC_FALSE;
128: #endif
129: }
130: }
131: return(0);
132: }
134: /* create PetscSubcomm object and determine num_solve_point (depends on npart, N, useconj) */
135: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUpSubComm(EPS eps,PetscInt *num_solve_point)
136: {
138: EPS_CISS *ctx = (EPS_CISS*)eps->data;
139: PetscInt N = ctx->N;
142: PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subcomm);
143: PetscSubcommSetNumber(ctx->subcomm,ctx->npart);
144: PetscSubcommSetType(ctx->subcomm,PETSC_SUBCOMM_INTERLACED);
145: PetscLogObjectMemory((PetscObject)eps,sizeof(PetscSubcomm));
146: ctx->subcomm_id = ctx->subcomm->color;
147: EPSCISSSetUseConj(eps,&ctx->useconj);
148: if (ctx->useconj) N = N/2;
149: *num_solve_point = N / ctx->npart;
150: if (N%ctx->npart > ctx->subcomm_id) (*num_solve_point)++;
151: return(0);
152: }
154: static PetscErrorCode CISSRedundantMat(EPS eps)
155: {
157: EPS_CISS *ctx = (EPS_CISS*)eps->data;
158: Mat A,B;
159: PetscInt nmat;
162: STGetNumMatrices(eps->st,&nmat);
163: if (ctx->subcomm->n != 1) {
164: STGetMatrix(eps->st,0,&A);
165: MatDestroy(&ctx->pA);
166: MatCreateRedundantMatrix(A,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pA);
167: if (nmat>1) {
168: STGetMatrix(eps->st,1,&B);
169: MatDestroy(&ctx->pB);
170: MatCreateRedundantMatrix(B,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pB);
171: } else ctx->pB = NULL;
172: } else {
173: ctx->pA = NULL;
174: ctx->pB = NULL;
175: }
176: return(0);
177: }
179: static PetscErrorCode CISSScatterVec(EPS eps)
180: {
182: EPS_CISS *ctx = (EPS_CISS*)eps->data;
183: IS is1,is2;
184: Vec v0;
185: PetscInt i,j,k,mstart,mend,mlocal;
186: PetscInt *idx1,*idx2,mloc_sub;
189: VecDestroy(&ctx->xsub);
190: MatCreateVecs(ctx->pA,&ctx->xsub,NULL);
192: VecDestroy(&ctx->xdup);
193: MatGetLocalSize(ctx->pA,&mloc_sub,NULL);
194: VecCreateMPI(PetscSubcommContiguousParent(ctx->subcomm),mloc_sub,PETSC_DECIDE,&ctx->xdup);
196: VecScatterDestroy(&ctx->scatterin);
197: BVGetColumn(ctx->V,0,&v0);
198: VecGetOwnershipRange(v0,&mstart,&mend);
199: mlocal = mend - mstart;
200: PetscMalloc2(ctx->subcomm->n*mlocal,&idx1,ctx->subcomm->n*mlocal,&idx2);
201: j = 0;
202: for (k=0;k<ctx->subcomm->n;k++) {
203: for (i=mstart;i<mend;i++) {
204: idx1[j] = i;
205: idx2[j++] = i + eps->n*k;
206: }
207: }
208: ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
209: ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
210: VecScatterCreate(v0,is1,ctx->xdup,is2,&ctx->scatterin);
211: ISDestroy(&is1);
212: ISDestroy(&is2);
213: PetscFree2(idx1,idx2);
214: BVRestoreColumn(ctx->V,0,&v0);
215: return(0);
216: }
218: static PetscErrorCode SetPathParameter(EPS eps)
219: {
221: EPS_CISS *ctx = (EPS_CISS*)eps->data;
222: PetscInt i,j;
223: PetscScalar center=0.0,tmp,tmp2,*omegai;
224: PetscReal theta,radius=1.0,vscale,a,b,c,d,max_w=0.0,rgscale;
225: #if defined(PETSC_USE_COMPLEX)
226: PetscReal start_ang,end_ang;
227: #endif
228: PetscBool isring=PETSC_FALSE,isellipse=PETSC_FALSE,isinterval=PETSC_FALSE;
231: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
232: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
233: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
234: RGGetScale(eps->rg,&rgscale);
235: PetscMalloc1(ctx->N+1l,&omegai);
236: RGComputeContour(eps->rg,ctx->N,ctx->omega,omegai);
237: if (isellipse) {
238: RGEllipseGetParameters(eps->rg,¢er,&radius,&vscale);
239: for (i=0;i<ctx->N;i++) {
240: #if defined(PETSC_USE_COMPLEX)
241: theta = 2.0*PETSC_PI*(i+0.5)/ctx->N;
242: ctx->pp[i] = PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
243: ctx->weight[i] = rgscale*radius*(PetscCMPLX(vscale*PetscCosReal(theta),PetscSinReal(theta)))/(PetscReal)ctx->N;
244: #else
245: theta = (PETSC_PI/ctx->N)*(i+0.5);
246: ctx->pp[i] = PetscCosReal(theta);
247: ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
248: ctx->omega[i] = rgscale*(center + radius*ctx->pp[i]);
249: #endif
250: }
251: } else if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
252: for (i=0;i<ctx->N;i++) {
253: theta = (PETSC_PI/ctx->N)*(i+0.5);
254: ctx->pp[i] = PetscCosReal(theta);
255: ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
256: }
257: if (isinterval) {
258: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
259: if ((c!=d || c!=0.0) && (a!=b || a!=0.0)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Endpoints of the imaginary axis or the real axis must be both zero");
260: for (i=0;i<ctx->N;i++) {
261: if (c==d) ctx->omega[i] = ((b-a)*(ctx->pp[i]+1.0)/2.0+a)*rgscale;
262: if (a==b) {
263: #if defined(PETSC_USE_COMPLEX)
264: ctx->omega[i] = ((d-c)*(ctx->pp[i]+1.0)/2.0+c)*rgscale*PETSC_i;
265: #else
266: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
267: #endif
268: }
269: }
270: }
271: if (isring) { /* only supported in complex scalars */
272: #if defined(PETSC_USE_COMPLEX)
273: RGRingGetParameters(eps->rg,¢er,&radius,&vscale,&start_ang,&end_ang,NULL);
274: for (i=0;i<ctx->N;i++) {
275: theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(ctx->pp[i])+1.0))*PETSC_PI;
276: ctx->omega[i] = rgscale*(center + radius*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta)));
277: }
278: #endif
279: }
280: } else {
281: if (isinterval) {
282: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
283: center = rgscale*((b+a)/2.0+(d+c)/2.0*PETSC_PI);
284: radius = PetscSqrtReal(PetscPowRealInt(rgscale*(b-a)/2.0,2)+PetscPowRealInt(rgscale*(d-c)/2.0,2));
285: } else if (isring) {
286: RGRingGetParameters(eps->rg,¢er,&radius,NULL,NULL,NULL,NULL);
287: center *= rgscale;
288: radius *= rgscale;
289: }
290: for (i=0;i<ctx->N;i++) {
291: ctx->pp[i] = (ctx->omega[i]-center)/radius;
292: tmp = 1; tmp2 = 1;
293: for (j=0;j<ctx->N;j++) {
294: tmp *= ctx->omega[j];
295: if (i != j) tmp2 *= ctx->omega[j]-ctx->omega[i];
296: }
297: ctx->weight[i] = tmp/tmp2;
298: max_w = PetscMax(PetscAbsScalar(ctx->weight[i]),max_w);
299: }
300: for (i=0;i<ctx->N;i++) ctx->weight[i] /= (PetscScalar)max_w;
301: }
302: PetscFree(omegai);
303: return(0);
304: }
306: static PetscErrorCode CISSVecSetRandom(BV V,PetscInt i0,PetscInt i1)
307: {
309: PetscInt i,j,nlocal;
310: PetscScalar *vdata;
311: Vec x;
314: BVGetSizes(V,&nlocal,NULL,NULL);
315: for (i=i0;i<i1;i++) {
316: BVSetRandomColumn(V,i);
317: BVGetColumn(V,i,&x);
318: VecGetArray(x,&vdata);
319: for (j=0;j<nlocal;j++) {
320: vdata[j] = PetscRealPart(vdata[j]);
321: if (PetscRealPart(vdata[j]) < 0.5) vdata[j] = -1.0;
322: else vdata[j] = 1.0;
323: }
324: VecRestoreArray(x,&vdata);
325: BVRestoreColumn(V,i,&x);
326: }
327: return(0);
328: }
330: static PetscErrorCode VecScatterVecs(EPS eps,BV Vin,PetscInt n)
331: {
332: PetscErrorCode ierr;
333: EPS_CISS *ctx = (EPS_CISS*)eps->data;
334: PetscInt i;
335: Vec vi,pvi;
336: const PetscScalar *array;
339: for (i=0;i<n;i++) {
340: BVGetColumn(Vin,i,&vi);
341: VecScatterBegin(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
342: VecScatterEnd(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
343: BVRestoreColumn(Vin,i,&vi);
344: VecGetArrayRead(ctx->xdup,&array);
345: VecPlaceArray(ctx->xsub,array);
346: BVGetColumn(ctx->pV,i,&pvi);
347: VecCopy(ctx->xsub,pvi);
348: BVRestoreColumn(ctx->pV,i,&pvi);
349: VecResetArray(ctx->xsub);
350: VecRestoreArrayRead(ctx->xdup,&array);
351: }
352: return(0);
353: }
355: static PetscErrorCode SolveLinearSystem(EPS eps,Mat A,Mat B,BV V,PetscInt L_start,PetscInt L_end,PetscBool initksp)
356: {
358: EPS_CISS *ctx = (EPS_CISS*)eps->data;
359: PetscInt i,j,p_id;
360: Mat Fz,kspMat;
361: Vec Bvj,vj,yj;
362: KSP ksp;
365: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
366: BVCreateVec(V,&Bvj);
367: if (ctx->usest) {
368: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Fz);
369: }
370: for (i=0;i<ctx->num_solve_point;i++) {
371: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
372: if (!ctx->usest && initksp) {
373: MatDuplicate(A,MAT_COPY_VALUES,&kspMat);
374: if (B) {
375: MatAXPY(kspMat,-ctx->omega[p_id],B,DIFFERENT_NONZERO_PATTERN);
376: } else {
377: MatShift(kspMat,-ctx->omega[p_id]);
378: }
379: KSPSetOperators(ctx->ksp[i],kspMat,kspMat);
380: MatDestroy(&kspMat);
381: } else if (ctx->usest) {
382: STSetShift(eps->st,ctx->omega[p_id]);
383: STGetKSP(eps->st,&ksp);
384: }
385: for (j=L_start;j<L_end;j++) {
386: BVGetColumn(V,j,&vj);
387: BVGetColumn(ctx->Y,i*ctx->L_max+j,&yj);
388: if (B) {
389: MatMult(B,vj,Bvj);
390: if (ctx->usest) {
391: KSPSolve(ksp,Bvj,yj);
392: } else {
393: KSPSolve(ctx->ksp[i],Bvj,yj);
394: }
395: } else {
396: if (ctx->usest) {
397: KSPSolve(ksp,vj,yj);
398: } else {
399: KSPSolve(ctx->ksp[i],vj,yj);
400: }
401: }
402: BVRestoreColumn(V,j,&vj);
403: BVRestoreColumn(ctx->Y,i*ctx->L_max+j,&yj);
404: }
405: if (ctx->usest && i<ctx->num_solve_point-1) { KSPReset(ksp); }
406: }
407: if (ctx->usest) { MatDestroy(&Fz); }
408: VecDestroy(&Bvj);
409: return(0);
410: }
412: #if defined(PETSC_USE_COMPLEX)
413: static PetscErrorCode EstimateNumberEigs(EPS eps,PetscInt *L_add)
414: {
416: EPS_CISS *ctx = (EPS_CISS*)eps->data;
417: PetscInt i,j,p_id;
418: PetscScalar tmp,m = 1,sum = 0.0;
419: PetscReal eta;
420: Vec v,vtemp,vj;
423: BVCreateVec(ctx->Y,&v);
424: BVCreateVec(ctx->V,&vtemp);
425: for (j=0;j<ctx->L;j++) {
426: VecSet(v,0);
427: for (i=0;i<ctx->num_solve_point; i++) {
428: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
429: BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
430: BVMultVec(ctx->Y,ctx->weight[p_id],1,v,&m);
431: }
432: BVGetColumn(ctx->V,j,&vj);
433: if (ctx->pA) {
434: VecSet(vtemp,0);
435: VecScatterBegin(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
436: VecScatterEnd(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
437: VecDot(vj,vtemp,&tmp);
438: } else {
439: VecDot(vj,v,&tmp);
440: }
441: BVRestoreColumn(ctx->V,j,&vj);
442: if (ctx->useconj) sum += PetscRealPart(tmp)*2;
443: else sum += tmp;
444: }
445: ctx->est_eig = PetscAbsScalar(sum/(PetscReal)ctx->L);
446: eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
447: PetscInfo1(eps,"Estimation_#Eig %f\n",(double)ctx->est_eig);
448: *L_add = (PetscInt)PetscCeilReal((ctx->est_eig*eta)/ctx->M) - ctx->L;
449: if (*L_add < 0) *L_add = 0;
450: if (*L_add>ctx->L_max-ctx->L) {
451: PetscInfo(eps,"Number of eigenvalues around the contour path may be too large\n");
452: *L_add = ctx->L_max-ctx->L;
453: }
454: VecDestroy(&v);
455: VecDestroy(&vtemp);
456: return(0);
457: }
458: #endif
460: static PetscErrorCode CalcMu(EPS eps,PetscScalar *Mu)
461: {
463: PetscMPIInt sub_size,len;
464: PetscInt i,j,k,s;
465: PetscScalar *m,*temp,*temp2,*ppk,alp;
466: EPS_CISS *ctx = (EPS_CISS*)eps->data;
467: Mat M;
470: MPI_Comm_size(PetscSubcommChild(ctx->subcomm),&sub_size);
471: PetscMalloc3(ctx->num_solve_point*ctx->L*(ctx->L+1),&temp,2*ctx->M*ctx->L*ctx->L,&temp2,ctx->num_solve_point,&ppk);
472: MatCreateSeqDense(PETSC_COMM_SELF,ctx->L,ctx->L_max*ctx->num_solve_point,NULL,&M);
473: for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] = 0;
474: BVSetActiveColumns(ctx->Y,0,ctx->L_max*ctx->num_solve_point);
475: if (ctx->pA) {
476: BVSetActiveColumns(ctx->pV,0,ctx->L);
477: BVDot(ctx->Y,ctx->pV,M);
478: } else {
479: BVSetActiveColumns(ctx->V,0,ctx->L);
480: BVDot(ctx->Y,ctx->V,M);
481: }
482: MatDenseGetArray(M,&m);
483: for (i=0;i<ctx->num_solve_point;i++) {
484: for (j=0;j<ctx->L;j++) {
485: for (k=0;k<ctx->L;k++) {
486: temp[k+j*ctx->L+i*ctx->L*ctx->L]=m[k+j*ctx->L+i*ctx->L*ctx->L_max];
487: }
488: }
489: }
490: MatDenseRestoreArray(M,&m);
491: for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
492: for (k=0;k<2*ctx->M;k++) {
493: for (j=0;j<ctx->L;j++) {
494: for (i=0;i<ctx->num_solve_point;i++) {
495: alp = ppk[i]*ctx->weight[i*ctx->subcomm->n + ctx->subcomm_id];
496: for (s=0;s<ctx->L;s++) {
497: if (ctx->useconj) temp2[s+(j+k*ctx->L)*ctx->L] += PetscRealPart(alp*temp[s+(j+i*ctx->L)*ctx->L])*2;
498: else temp2[s+(j+k*ctx->L)*ctx->L] += alp*temp[s+(j+i*ctx->L)*ctx->L];
499: }
500: }
501: }
502: for (i=0;i<ctx->num_solve_point;i++)
503: ppk[i] *= ctx->pp[i*ctx->subcomm->n + ctx->subcomm_id];
504: }
505: for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] /= sub_size;
506: PetscMPIIntCast(2*ctx->M*ctx->L*ctx->L,&len);
507: MPI_Allreduce(temp2,Mu,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)eps));
508: PetscFree3(temp,temp2,ppk);
509: MatDestroy(&M);
510: return(0);
511: }
513: static PetscErrorCode BlockHankel(EPS eps,PetscScalar *Mu,PetscInt s,PetscScalar *H)
514: {
515: EPS_CISS *ctx = (EPS_CISS*)eps->data;
516: PetscInt i,j,k,L=ctx->L,M=ctx->M;
519: for (k=0;k<L*M;k++)
520: for (j=0;j<M;j++)
521: for (i=0;i<L;i++)
522: H[j*L+i+k*L*M] = Mu[i+k*L+(j+s)*L*L];
523: return(0);
524: }
526: static PetscErrorCode SVD_H0(EPS eps,PetscScalar *S,PetscInt *K)
527: {
529: EPS_CISS *ctx = (EPS_CISS*)eps->data;
530: PetscInt i,ml=ctx->L*ctx->M;
531: PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
532: PetscScalar *work;
533: #if defined(PETSC_USE_COMPLEX)
534: PetscReal *rwork;
535: #endif
538: PetscMalloc1(5*ml,&work);
539: #if defined(PETSC_USE_COMPLEX)
540: PetscMalloc1(5*ml,&rwork);
541: #endif
542: PetscBLASIntCast(ml,&m);
543: n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
544: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
545: #if defined(PETSC_USE_COMPLEX)
546: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
547: #else
548: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
549: #endif
550: SlepcCheckLapackInfo("gesvd",info);
551: PetscFPTrapPop();
552: (*K) = 0;
553: for (i=0;i<ml;i++) {
554: if (ctx->sigma[i]/PetscMax(ctx->sigma[0],1)>ctx->delta) (*K)++;
555: }
556: PetscFree(work);
557: #if defined(PETSC_USE_COMPLEX)
558: PetscFree(rwork);
559: #endif
560: return(0);
561: }
563: static PetscErrorCode ConstructS(EPS eps)
564: {
566: EPS_CISS *ctx = (EPS_CISS*)eps->data;
567: PetscInt i,j,k,vec_local_size,p_id;
568: Vec v,sj;
569: PetscScalar *ppk, *v_data, m = 1;
572: BVGetSizes(ctx->Y,&vec_local_size,NULL,NULL);
573: PetscMalloc1(ctx->num_solve_point,&ppk);
574: for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
575: BVCreateVec(ctx->Y,&v);
576: for (k=0;k<ctx->M;k++) {
577: for (j=0;j<ctx->L;j++) {
578: VecSet(v,0);
579: for (i=0;i<ctx->num_solve_point;i++) {
580: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
581: BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
582: BVMultVec(ctx->Y,ppk[i]*ctx->weight[p_id],1.0,v,&m);
583: }
584: if (ctx->useconj) {
585: VecGetArray(v,&v_data);
586: for (i=0;i<vec_local_size;i++) v_data[i] = PetscRealPart(v_data[i])*2;
587: VecRestoreArray(v,&v_data);
588: }
589: BVGetColumn(ctx->S,k*ctx->L+j,&sj);
590: if (ctx->pA) {
591: VecSet(sj,0);
592: VecScatterBegin(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
593: VecScatterEnd(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
594: } else {
595: VecCopy(v,sj);
596: }
597: BVRestoreColumn(ctx->S,k*ctx->L+j,&sj);
598: }
599: for (i=0;i<ctx->num_solve_point;i++) {
600: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
601: ppk[i] *= ctx->pp[p_id];
602: }
603: }
604: PetscFree(ppk);
605: VecDestroy(&v);
606: return(0);
607: }
609: static PetscErrorCode SVD_S(BV S,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *K)
610: {
612: PetscInt i,j,k,local_size;
613: PetscMPIInt len;
614: PetscScalar *work,*temp,*B,*tempB,*s_data,*Q1,*Q2,*temp2,alpha=1,beta=0;
615: PetscBLASInt l,m,n,lda,ldu,ldvt,lwork,info,ldb,ldc;
616: #if defined(PETSC_USE_COMPLEX)
617: PetscReal *rwork;
618: #endif
621: BVGetSizes(S,&local_size,NULL,NULL);
622: BVGetArray(S,&s_data);
623: PetscMalloc7(ml*ml,&temp,ml*ml,&temp2,local_size*ml,&Q1,local_size*ml,&Q2,ml*ml,&B,ml*ml,&tempB,5*ml,&work);
624: PetscArrayzero(B,ml*ml);
625: #if defined(PETSC_USE_COMPLEX)
626: PetscMalloc1(5*ml,&rwork);
627: #endif
628: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
630: for (i=0;i<ml;i++) B[i*ml+i]=1;
632: for (k=0;k<2;k++) {
633: PetscBLASIntCast(local_size,&m);
634: PetscBLASIntCast(ml,&l);
635: n = l; lda = m; ldb = m; ldc = l;
636: if (k == 0) {
637: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,s_data,&lda,s_data,&ldb,&beta,temp,&ldc));
638: } else if ((k%2)==1) {
639: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,temp,&ldc));
640: } else {
641: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q2,&lda,Q2,&ldb,&beta,temp,&ldc));
642: }
643: PetscArrayzero(temp2,ml*ml);
644: PetscMPIIntCast(ml*ml,&len);
645: MPI_Allreduce(temp,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S));
647: PetscBLASIntCast(ml,&m);
648: n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
649: #if defined(PETSC_USE_COMPLEX)
650: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
651: #else
652: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
653: #endif
654: SlepcCheckLapackInfo("gesvd",info);
656: PetscBLASIntCast(local_size,&l);
657: PetscBLASIntCast(ml,&n);
658: m = n; lda = l; ldb = m; ldc = l;
659: if (k==0) {
660: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,s_data,&lda,temp2,&ldb,&beta,Q1,&ldc));
661: } else if ((k%2)==1) {
662: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
663: } else {
664: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q2,&lda,temp2,&ldb,&beta,Q1,&ldc));
665: }
667: PetscBLASIntCast(ml,&l);
668: m = l; n = l; lda = l; ldb = m; ldc = l;
669: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
670: for (i=0;i<ml;i++) {
671: sigma[i] = sqrt(sigma[i]);
672: for (j=0;j<local_size;j++) {
673: if ((k%2)==1) Q2[j+i*local_size]/=sigma[i];
674: else Q1[j+i*local_size]/=sigma[i];
675: }
676: for (j=0;j<ml;j++) {
677: B[j+i*ml]=tempB[j+i*ml]*sigma[i];
678: }
679: }
680: }
682: PetscBLASIntCast(ml,&m);
683: n = m; lda = m; ldu=1; ldvt=1;
684: #if defined(PETSC_USE_COMPLEX)
685: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
686: #else
687: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
688: #endif
689: SlepcCheckLapackInfo("gesvd",info);
691: PetscBLASIntCast(local_size,&l);
692: PetscBLASIntCast(ml,&n);
693: m = n; lda = l; ldb = m; ldc = l;
694: if ((k%2)==1) {
695: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,s_data,&ldc));
696: } else {
697: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,s_data,&ldc));
698: }
700: PetscFPTrapPop();
701: BVRestoreArray(S,&s_data);
703: (*K) = 0;
704: for (i=0;i<ml;i++) {
705: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*K)++;
706: }
707: PetscFree7(temp,temp2,Q1,Q2,B,tempB,work);
708: #if defined(PETSC_USE_COMPLEX)
709: PetscFree(rwork);
710: #endif
711: return(0);
712: }
714: static PetscErrorCode isGhost(EPS eps,PetscInt ld,PetscInt nv,PetscBool *fl)
715: {
717: EPS_CISS *ctx = (EPS_CISS*)eps->data;
718: PetscInt i,j;
719: PetscScalar *pX;
720: PetscReal *tau,s1,s2,tau_max=0.0;
723: PetscMalloc1(nv,&tau);
724: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
725: DSGetArray(eps->ds,DS_MAT_X,&pX);
727: for (i=0;i<nv;i++) {
728: s1 = 0;
729: s2 = 0;
730: for (j=0;j<nv;j++) {
731: s1 += PetscAbsScalar(PetscPowScalarInt(pX[i*ld+j],2));
732: s2 += PetscPowRealInt(PetscAbsScalar(pX[i*ld+j]),2)/ctx->sigma[j];
733: }
734: tau[i] = s1/s2;
735: tau_max = PetscMax(tau_max,tau[i]);
736: }
737: DSRestoreArray(eps->ds,DS_MAT_X,&pX);
738: for (i=0;i<nv;i++) {
739: tau[i] /= tau_max;
740: }
741: for (i=0;i<nv;i++) {
742: if (tau[i]>=ctx->spurious_threshold) fl[i] = PETSC_TRUE;
743: else fl[i] = PETSC_FALSE;
744: }
745: PetscFree(tau);
746: return(0);
747: }
749: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
750: {
752: EPS_CISS *ctx = (EPS_CISS*)eps->data;
753: PetscInt i;
754: PetscScalar center;
755: PetscReal radius,a,b,c,d,rgscale;
756: #if defined(PETSC_USE_COMPLEX)
757: PetscReal start_ang,end_ang,vscale,theta;
758: #endif
759: PetscBool isring,isellipse,isinterval;
762: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
763: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
764: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
765: RGGetScale(eps->rg,&rgscale);
766: if (isinterval) {
767: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
768: if (c==d) {
769: for (i=0;i<nv;i++) {
770: #if defined(PETSC_USE_COMPLEX)
771: eps->eigr[i] = PetscRealPart(eps->eigr[i]);
772: #else
773: eps->eigi[i] = 0;
774: #endif
775: }
776: }
777: }
778: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
779: if (isellipse) {
780: RGEllipseGetParameters(eps->rg,¢er,&radius,NULL);
781: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
782: } else if (isinterval) {
783: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
784: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
785: for (i=0;i<nv;i++) {
786: if (c==d) eps->eigr[i] = ((b-a)*(eps->eigr[i]+1.0)/2.0+a)*rgscale;
787: if (a==b) {
788: #if defined(PETSC_USE_COMPLEX)
789: eps->eigr[i] = ((d-c)*(eps->eigr[i]+1.0)/2.0+c)*rgscale*PETSC_i;
790: #else
791: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
792: #endif
793: }
794: }
795: } else {
796: center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
797: radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
798: for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
799: }
800: } else if (isring) { /* only supported in complex scalars */
801: #if defined(PETSC_USE_COMPLEX)
802: RGRingGetParameters(eps->rg,¢er,&radius,&vscale,&start_ang,&end_ang,NULL);
803: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
804: for (i=0;i<nv;i++) {
805: theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
806: eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
807: }
808: } else {
809: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
810: }
811: #endif
812: }
813: }
814: return(0);
815: }
817: PetscErrorCode EPSSetUp_CISS(EPS eps)
818: {
820: EPS_CISS *ctx = (EPS_CISS*)eps->data;
821: PetscBool issinvert,istrivial,isring,isellipse,isinterval,flg,useconj;
822: PetscReal c,d;
823: Mat A;
826: if (!eps->ncv) {
827: eps->ncv = ctx->L_max*ctx->M;
828: if (eps->ncv>eps->n) {
829: eps->ncv = eps->n;
830: ctx->L_max = eps->ncv/ctx->M;
831: if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
832: }
833: } else {
834: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
835: ctx->L_max = eps->ncv/ctx->M;
836: if (!ctx->L_max) {
837: ctx->L_max = 1;
838: eps->ncv = ctx->L_max*ctx->M;
839: }
840: }
841: ctx->L = PetscMin(ctx->L,ctx->L_max);
842: if (!eps->max_it) eps->max_it = 1;
843: if (!eps->mpd) eps->mpd = eps->ncv;
844: if (!eps->which) eps->which = EPS_ALL;
845: if (!eps->extraction) { EPSSetExtraction(eps,EPS_RITZ); }
846: else if (eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
847: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
848: if (eps->stopping!=EPSStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support user-defined stopping test");
850: /* check region */
851: RGIsTrivial(eps->rg,&istrivial);
852: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
853: RGGetComplement(eps->rg,&flg);
854: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
855: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
856: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
857: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
858: if (!isellipse && !isring && !isinterval) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");
859: /* if useconj has changed, then reset subcomm data */
860: EPSCISSSetUseConj(eps,&useconj);
861: if (useconj!=ctx->useconj) { EPSCISSResetSubcomm(eps); }
863: #if !defined(PETSC_USE_COMPLEX)
864: if (isring) {
865: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
866: }
867: #endif
868: if (isinterval) {
869: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
870: #if !defined(PETSC_USE_COMPLEX)
871: if (c!=d || c!=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"In real scalars, endpoints of the imaginary axis must be both zero");
872: #endif
873: if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
874: }
875: if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;
877: /* create split comm */
878: if (!ctx->subcomm) { EPSCISSSetUpSubComm(eps,&ctx->num_solve_point); }
880: EPSAllocateSolution(eps,0);
881: if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
882: PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
883: PetscLogObjectMemory((PetscObject)eps,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));
885: /* allocate basis vectors */
886: BVDestroy(&ctx->S);
887: BVDuplicateResize(eps->V,ctx->L_max*ctx->M,&ctx->S);
888: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->S);
889: BVDestroy(&ctx->V);
890: BVDuplicateResize(eps->V,ctx->L_max,&ctx->V);
891: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->V);
893: STGetMatrix(eps->st,0,&A);
894: PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
895: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
897: if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
898: if (ctx->usest && ctx->npart>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");
900: CISSRedundantMat(eps);
901: if (ctx->pA) {
902: CISSScatterVec(eps);
903: BVDestroy(&ctx->pV);
904: BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->pV);
905: BVSetSizesFromVec(ctx->pV,ctx->xsub,eps->n);
906: BVSetFromOptions(ctx->pV);
907: BVResize(ctx->pV,ctx->L_max,PETSC_FALSE);
908: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->pV);
909: }
911: if (ctx->usest) {
912: PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&issinvert);
913: if (!issinvert) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"If the usest flag is set, you must select the STSINVERT spectral transformation");
914: }
916: BVDestroy(&ctx->Y);
917: if (ctx->pA) {
918: BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->Y);
919: BVSetSizesFromVec(ctx->Y,ctx->xsub,eps->n);
920: BVSetFromOptions(ctx->Y);
921: BVResize(ctx->Y,ctx->num_solve_point*ctx->L_max,PETSC_FALSE);
922: } else {
923: BVDuplicateResize(eps->V,ctx->num_solve_point*ctx->L_max,&ctx->Y);
924: }
925: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->Y);
927: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
928: DSSetType(eps->ds,DSGNHEP);
929: } else if (eps->isgeneralized) {
930: if (eps->ishermitian && eps->ispositive) {
931: DSSetType(eps->ds,DSGHEP);
932: } else {
933: DSSetType(eps->ds,DSGNHEP);
934: }
935: } else {
936: if (eps->ishermitian) {
937: DSSetType(eps->ds,DSHEP);
938: } else {
939: DSSetType(eps->ds,DSNHEP);
940: }
941: }
942: DSAllocate(eps->ds,eps->ncv);
943: EPSSetWorkVecs(eps,2);
945: #if !defined(PETSC_USE_COMPLEX)
946: if (!eps->ishermitian) { PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"); }
947: #endif
948: return(0);
949: }
951: PetscErrorCode EPSSolve_CISS(EPS eps)
952: {
954: EPS_CISS *ctx = (EPS_CISS*)eps->data;
955: Mat A,B,X,M,pA,pB;
956: PetscInt i,j,ld,nmat,L_add=0,nv=0,L_base=ctx->L,inner,nlocal,*inside;
957: PetscScalar *Mu,*H0,*H1=NULL,*rr,*temp;
958: PetscReal error,max_error,norm;
959: PetscBool *fl1;
960: Vec si,w[3];
961: SlepcSC sc;
962: PetscRandom rand;
963: #if defined(PETSC_USE_COMPLEX)
964: PetscBool isellipse;
965: #endif
968: w[0] = eps->work[0];
969: w[1] = NULL;
970: w[2] = eps->work[1];
971: /* override SC settings */
972: DSGetSlepcSC(eps->ds,&sc);
973: sc->comparison = SlepcCompareLargestMagnitude;
974: sc->comparisonctx = NULL;
975: sc->map = NULL;
976: sc->mapobj = NULL;
977: VecGetLocalSize(w[0],&nlocal);
978: DSGetLeadingDimension(eps->ds,&ld);
979: STGetNumMatrices(eps->st,&nmat);
980: STGetMatrix(eps->st,0,&A);
981: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
982: else B = NULL;
983: SetPathParameter(eps);
984: CISSVecSetRandom(ctx->V,0,ctx->L);
985: BVGetRandomContext(ctx->V,&rand);
987: if (ctx->pA) {
988: VecScatterVecs(eps,ctx->V,ctx->L);
989: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_TRUE);
990: } else {
991: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_TRUE);
992: }
993: #if defined(PETSC_USE_COMPLEX)
994: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
995: if (isellipse) {
996: EstimateNumberEigs(eps,&L_add);
997: } else {
998: L_add = 0;
999: }
1000: #else
1001: L_add = 0;
1002: #endif
1003: if (L_add>0) {
1004: PetscInfo2(eps,"Changing L %D -> %D by Estimate #Eig\n",ctx->L,ctx->L+L_add);
1005: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1006: if (ctx->pA) {
1007: VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1008: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1009: } else {
1010: SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1011: }
1012: ctx->L += L_add;
1013: }
1014: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
1015: for (i=0;i<ctx->refine_blocksize;i++) {
1016: CalcMu(eps,Mu);
1017: BlockHankel(eps,Mu,0,H0);
1018: SVD_H0(eps,H0,&nv);
1019: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
1020: L_add = L_base;
1021: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
1022: PetscInfo2(eps,"Changing L %D -> %D by SVD(H0)\n",ctx->L,ctx->L+L_add);
1023: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1024: if (ctx->pA) {
1025: VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1026: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1027: } else {
1028: SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1029: }
1030: ctx->L += L_add;
1031: }
1032: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1033: PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);
1034: }
1036: while (eps->reason == EPS_CONVERGED_ITERATING) {
1037: eps->its++;
1038: for (inner=0;inner<=ctx->refine_inner;inner++) {
1039: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1040: CalcMu(eps,Mu);
1041: BlockHankel(eps,Mu,0,H0);
1042: SVD_H0(eps,H0,&nv);
1043: break;
1044: } else {
1045: ConstructS(eps);
1046: BVSetActiveColumns(ctx->S,0,ctx->L);
1047: BVCopy(ctx->S,ctx->V);
1048: SVD_S(ctx->S,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
1049: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
1050: if (ctx->pA) {
1051: VecScatterVecs(eps,ctx->V,ctx->L);
1052: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1053: } else {
1054: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1055: }
1056: } else break;
1057: }
1058: }
1059: eps->nconv = 0;
1060: if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
1061: else {
1062: DSSetDimensions(eps->ds,nv,0,0,0);
1063: DSSetState(eps->ds,DS_STATE_RAW);
1065: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1066: BlockHankel(eps,Mu,0,H0);
1067: BlockHankel(eps,Mu,1,H1);
1068: DSGetArray(eps->ds,DS_MAT_A,&temp);
1069: for (j=0;j<nv;j++) {
1070: for (i=0;i<nv;i++) {
1071: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
1072: }
1073: }
1074: DSRestoreArray(eps->ds,DS_MAT_A,&temp);
1075: DSGetArray(eps->ds,DS_MAT_B,&temp);
1076: for (j=0;j<nv;j++) {
1077: for (i=0;i<nv;i++) {
1078: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
1079: }
1080: }
1081: DSRestoreArray(eps->ds,DS_MAT_B,&temp);
1082: } else {
1083: BVSetActiveColumns(ctx->S,0,nv);
1084: DSGetMat(eps->ds,DS_MAT_A,&pA);
1085: MatZeroEntries(pA);
1086: BVMatProject(ctx->S,A,ctx->S,pA);
1087: DSRestoreMat(eps->ds,DS_MAT_A,&pA);
1088: if (B) {
1089: DSGetMat(eps->ds,DS_MAT_B,&pB);
1090: MatZeroEntries(pB);
1091: BVMatProject(ctx->S,B,ctx->S,pB);
1092: DSRestoreMat(eps->ds,DS_MAT_B,&pB);
1093: }
1094: }
1096: DSSolve(eps->ds,eps->eigr,eps->eigi);
1097: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
1099: PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
1100: rescale_eig(eps,nv);
1101: isGhost(eps,ld,nv,fl1);
1102: RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside);
1103: for (i=0;i<nv;i++) {
1104: if (fl1[i] && inside[i]>=0) {
1105: rr[i] = 1.0;
1106: eps->nconv++;
1107: } else rr[i] = 0.0;
1108: }
1109: DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv);
1110: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
1111: rescale_eig(eps,nv);
1112: PetscFree3(fl1,inside,rr);
1113: BVSetActiveColumns(eps->V,0,nv);
1114: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1115: ConstructS(eps);
1116: BVSetActiveColumns(ctx->S,0,ctx->L);
1117: BVCopy(ctx->S,ctx->V);
1118: BVSetActiveColumns(ctx->S,0,nv);
1119: }
1120: BVCopy(ctx->S,eps->V);
1122: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
1123: DSGetMat(eps->ds,DS_MAT_X,&X);
1124: BVMultInPlace(ctx->S,X,0,eps->nconv);
1125: if (eps->ishermitian) {
1126: BVMultInPlace(eps->V,X,0,eps->nconv);
1127: }
1128: MatDestroy(&X);
1129: max_error = 0.0;
1130: for (i=0;i<eps->nconv;i++) {
1131: BVGetColumn(ctx->S,i,&si);
1132: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,NULL,w,&error);
1133: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) { /* vector is not normalized */
1134: VecNorm(si,NORM_2,&norm);
1135: error /= norm;
1136: }
1137: (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx);
1138: BVRestoreColumn(ctx->S,i,&si);
1139: max_error = PetscMax(max_error,error);
1140: }
1142: if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
1143: else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
1144: else {
1145: if (eps->nconv > ctx->L) {
1146: MatCreateSeqDense(PETSC_COMM_SELF,eps->nconv,ctx->L,NULL,&M);
1147: MatDenseGetArray(M,&temp);
1148: for (i=0;i<ctx->L*eps->nconv;i++) {
1149: PetscRandomGetValue(rand,&temp[i]);
1150: temp[i] = PetscRealPart(temp[i]);
1151: }
1152: MatDenseRestoreArray(M,&temp);
1153: BVSetActiveColumns(ctx->S,0,eps->nconv);
1154: BVMultInPlace(ctx->S,M,0,ctx->L);
1155: MatDestroy(&M);
1156: BVSetActiveColumns(ctx->S,0,ctx->L);
1157: BVCopy(ctx->S,ctx->V);
1158: }
1159: if (ctx->pA) {
1160: VecScatterVecs(eps,ctx->V,ctx->L);
1161: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1162: } else {
1163: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1164: }
1165: }
1166: }
1167: }
1168: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1169: PetscFree(H1);
1170: }
1171: PetscFree2(Mu,H0);
1172: return(0);
1173: }
1175: PetscErrorCode EPSComputeVectors_CISS(EPS eps)
1176: {
1178: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1179: PetscReal norm;
1180: PetscInt i;
1181: Mat B=NULL;
1184: EPSComputeVectors_Schur(eps);
1185: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) { /* normalize */
1186: if (eps->isgeneralized && eps->ishermitian && eps->ispositive) {
1187: STGetMatrix(eps->st,1,&B);
1188: BVSetMatrix(eps->V,B,PETSC_FALSE);
1189: }
1190: for (i=0;i<eps->nconv;i++) {
1191: BVNormColumn(eps->V,i,NORM_2,&norm);
1192: BVScaleColumn(eps->V,i,1.0/norm);
1193: }
1194: if (B) { BVSetMatrix(eps->V,NULL,PETSC_FALSE); }
1195: }
1196: return(0);
1197: }
1199: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1200: {
1202: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1203: PetscInt oN,onpart;
1206: oN = ctx->N;
1207: if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
1208: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
1209: } else {
1210: if (ip<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
1211: if (ip%2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
1212: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
1213: }
1214: if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
1215: ctx->L = 16;
1216: } else {
1217: if (bs<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
1218: ctx->L = bs;
1219: }
1220: if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
1221: ctx->M = ctx->N/4;
1222: } else {
1223: if (ms<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
1224: if (ms>ctx->N) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be less than or equal to the number of integration points");
1225: ctx->M = ms;
1226: }
1227: onpart = ctx->npart;
1228: if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
1229: ctx->npart = 1;
1230: } else {
1231: if (npart<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
1232: ctx->npart = npart;
1233: }
1234: if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
1235: ctx->L_max = 64;
1236: } else {
1237: if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
1238: ctx->L_max = PetscMax(bsmax,ctx->L);
1239: }
1240: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) { EPSCISSResetSubcomm(eps); }
1241: ctx->isreal = realmats;
1242: eps->state = EPS_STATE_INITIAL;
1243: return(0);
1244: }
1246: /*@
1247: EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.
1249: Logically Collective on eps
1251: Input Parameters:
1252: + eps - the eigenproblem solver context
1253: . ip - number of integration points
1254: . bs - block size
1255: . ms - moment size
1256: . npart - number of partitions when splitting the communicator
1257: . bsmax - max block size
1258: - realmats - A and B are real
1260: Options Database Keys:
1261: + -eps_ciss_integration_points - Sets the number of integration points
1262: . -eps_ciss_blocksize - Sets the block size
1263: . -eps_ciss_moments - Sets the moment size
1264: . -eps_ciss_partitions - Sets the number of partitions
1265: . -eps_ciss_maxblocksize - Sets the maximum block size
1266: - -eps_ciss_realmats - A and B are real
1268: Note:
1269: The default number of partitions is 1. This means the internal KSP object is shared
1270: among all processes of the EPS communicator. Otherwise, the communicator is split
1271: into npart communicators, so that npart KSP solves proceed simultaneously.
1273: Level: advanced
1275: .seealso: EPSCISSGetSizes()
1276: @*/
1277: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1278: {
1289: PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
1290: return(0);
1291: }
1293: static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1294: {
1295: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1298: if (ip) *ip = ctx->N;
1299: if (bs) *bs = ctx->L;
1300: if (ms) *ms = ctx->M;
1301: if (npart) *npart = ctx->npart;
1302: if (bsmax) *bsmax = ctx->L_max;
1303: if (realmats) *realmats = ctx->isreal;
1304: return(0);
1305: }
1307: /*@
1308: EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.
1310: Not Collective
1312: Input Parameter:
1313: . eps - the eigenproblem solver context
1315: Output Parameters:
1316: + ip - number of integration points
1317: . bs - block size
1318: . ms - moment size
1319: . npart - number of partitions when splitting the communicator
1320: . bsmax - max block size
1321: - realmats - A and B are real
1323: Level: advanced
1325: .seealso: EPSCISSSetSizes()
1326: @*/
1327: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1328: {
1333: PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
1334: return(0);
1335: }
1337: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
1338: {
1339: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1342: if (delta == PETSC_DEFAULT) {
1343: ctx->delta = 1e-12;
1344: } else {
1345: if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
1346: ctx->delta = delta;
1347: }
1348: if (spur == PETSC_DEFAULT) {
1349: ctx->spurious_threshold = 1e-4;
1350: } else {
1351: if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
1352: ctx->spurious_threshold = spur;
1353: }
1354: return(0);
1355: }
1357: /*@
1358: EPSCISSSetThreshold - Sets the values of various threshold parameters in
1359: the CISS solver.
1361: Logically Collective on eps
1363: Input Parameters:
1364: + eps - the eigenproblem solver context
1365: . delta - threshold for numerical rank
1366: - spur - spurious threshold (to discard spurious eigenpairs)
1368: Options Database Keys:
1369: + -eps_ciss_delta - Sets the delta
1370: - -eps_ciss_spurious_threshold - Sets the spurious threshold
1372: Level: advanced
1374: .seealso: EPSCISSGetThreshold()
1375: @*/
1376: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
1377: {
1384: PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
1385: return(0);
1386: }
1388: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
1389: {
1390: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1393: if (delta) *delta = ctx->delta;
1394: if (spur) *spur = ctx->spurious_threshold;
1395: return(0);
1396: }
1398: /*@
1399: EPSCISSGetThreshold - Gets the values of various threshold parameters
1400: in the CISS solver.
1402: Not Collective
1404: Input Parameter:
1405: . eps - the eigenproblem solver context
1407: Output Parameters:
1408: + delta - threshold for numerical rank
1409: - spur - spurious threshold (to discard spurious eigenpairs)
1411: Level: advanced
1413: .seealso: EPSCISSSetThreshold()
1414: @*/
1415: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
1416: {
1421: PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
1422: return(0);
1423: }
1425: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
1426: {
1427: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1430: if (inner == PETSC_DEFAULT) {
1431: ctx->refine_inner = 0;
1432: } else {
1433: if (inner<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
1434: ctx->refine_inner = inner;
1435: }
1436: if (blsize == PETSC_DEFAULT) {
1437: ctx->refine_blocksize = 0;
1438: } else {
1439: if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
1440: ctx->refine_blocksize = blsize;
1441: }
1442: return(0);
1443: }
1445: /*@
1446: EPSCISSSetRefinement - Sets the values of various refinement parameters
1447: in the CISS solver.
1449: Logically Collective on eps
1451: Input Parameters:
1452: + eps - the eigenproblem solver context
1453: . inner - number of iterative refinement iterations (inner loop)
1454: - blsize - number of iterative refinement iterations (blocksize loop)
1456: Options Database Keys:
1457: + -eps_ciss_refine_inner - Sets number of inner iterations
1458: - -eps_ciss_refine_blocksize - Sets number of blocksize iterations
1460: Level: advanced
1462: .seealso: EPSCISSGetRefinement()
1463: @*/
1464: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
1465: {
1472: PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
1473: return(0);
1474: }
1476: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
1477: {
1478: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1481: if (inner) *inner = ctx->refine_inner;
1482: if (blsize) *blsize = ctx->refine_blocksize;
1483: return(0);
1484: }
1486: /*@
1487: EPSCISSGetRefinement - Gets the values of various refinement parameters
1488: in the CISS solver.
1490: Not Collective
1492: Input Parameter:
1493: . eps - the eigenproblem solver context
1495: Output Parameters:
1496: + inner - number of iterative refinement iterations (inner loop)
1497: - blsize - number of iterative refinement iterations (blocksize loop)
1499: Level: advanced
1501: .seealso: EPSCISSSetRefinement()
1502: @*/
1503: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
1504: {
1509: PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
1510: return(0);
1511: }
1513: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
1514: {
1515: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1518: ctx->usest = usest;
1519: ctx->usest_set = PETSC_TRUE;
1520: eps->state = EPS_STATE_INITIAL;
1521: return(0);
1522: }
1524: /*@
1525: EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
1526: use the ST object for the linear solves.
1528: Logically Collective on eps
1530: Input Parameters:
1531: + eps - the eigenproblem solver context
1532: - usest - boolean flag to use the ST object or not
1534: Options Database Keys:
1535: . -eps_ciss_usest <bool> - whether the ST object will be used or not
1537: Level: advanced
1539: .seealso: EPSCISSGetUseST()
1540: @*/
1541: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
1542: {
1548: PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1549: return(0);
1550: }
1552: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1553: {
1554: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1557: *usest = ctx->usest;
1558: return(0);
1559: }
1561: /*@
1562: EPSCISSGetUseST - Gets the flag for using the ST object
1563: in the CISS solver.
1565: Not Collective
1567: Input Parameter:
1568: . eps - the eigenproblem solver context
1570: Output Parameters:
1571: . usest - boolean flag indicating if the ST object is being used
1573: Level: advanced
1575: .seealso: EPSCISSSetUseST()
1576: @*/
1577: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1578: {
1584: PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1585: return(0);
1586: }
1588: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1589: {
1590: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1593: ctx->quad = quad;
1594: return(0);
1595: }
1597: /*@
1598: EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.
1600: Logically Collective on eps
1602: Input Parameters:
1603: + eps - the eigenproblem solver context
1604: - quad - the quadrature rule
1606: Options Database Key:
1607: . -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1608: 'chebyshev')
1610: Notes:
1611: By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).
1613: If the 'chebyshev' option is specified (EPS_CISS_QUADRULE_CHEBYSHEV), then
1614: Chebyshev points are used as quadrature points.
1616: Level: advanced
1618: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1619: @*/
1620: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1621: {
1627: PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1628: return(0);
1629: }
1631: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1632: {
1633: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1636: *quad = ctx->quad;
1637: return(0);
1638: }
1640: /*@
1641: EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.
1643: Not Collective
1645: Input Parameter:
1646: . eps - the eigenproblem solver context
1648: Output Parameters:
1649: . quad - quadrature rule
1651: Level: advanced
1653: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1654: @*/
1655: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1656: {
1662: PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1663: return(0);
1664: }
1666: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1667: {
1668: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1671: ctx->extraction = extraction;
1672: return(0);
1673: }
1675: /*@
1676: EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.
1678: Logically Collective on eps
1680: Input Parameters:
1681: + eps - the eigenproblem solver context
1682: - extraction - the extraction technique
1684: Options Database Key:
1685: . -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1686: 'hankel')
1688: Notes:
1689: By default, the Rayleigh-Ritz extraction is used (EPS_CISS_EXTRACTION_RITZ).
1691: If the 'hankel' option is specified (EPS_CISS_EXTRACTION_HANKEL), then
1692: the Block Hankel method is used for extracting eigenpairs.
1694: Level: advanced
1696: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1697: @*/
1698: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1699: {
1705: PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1706: return(0);
1707: }
1709: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1710: {
1711: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1714: *extraction = ctx->extraction;
1715: return(0);
1716: }
1718: /*@
1719: EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.
1721: Not Collective
1723: Input Parameter:
1724: . eps - the eigenproblem solver context
1726: Output Parameters:
1727: + extraction - extraction technique
1729: Level: advanced
1731: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1732: @*/
1733: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1734: {
1740: PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1741: return(0);
1742: }
1744: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1745: {
1747: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1748: PetscInt i;
1749: PC pc;
1752: if (!ctx->ksp) {
1753: if (!ctx->subcomm) { /* initialize subcomm first */
1754: EPSCISSSetUseConj(eps,&ctx->useconj);
1755: EPSCISSSetUpSubComm(eps,&ctx->num_solve_point);
1756: }
1757: PetscMalloc1(ctx->num_solve_point,&ctx->ksp);
1758: for (i=0;i<ctx->num_solve_point;i++) {
1759: KSPCreate(PetscSubcommChild(ctx->subcomm),&ctx->ksp[i]);
1760: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)eps,1);
1761: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)eps)->prefix);
1762: KSPAppendOptionsPrefix(ctx->ksp[i],"eps_ciss_");
1763: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ksp[i]);
1764: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)eps)->options);
1765: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1766: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1767: KSPGetPC(ctx->ksp[i],&pc);
1768: KSPSetType(ctx->ksp[i],KSPPREONLY);
1769: PCSetType(pc,PCLU);
1770: }
1771: }
1772: if (nsolve) *nsolve = ctx->num_solve_point;
1773: if (ksp) *ksp = ctx->ksp;
1774: return(0);
1775: }
1777: /*@C
1778: EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1779: the CISS solver.
1781: Not Collective
1783: Input Parameter:
1784: . eps - the eigenproblem solver solver
1786: Output Parameters:
1787: + nsolve - number of solver objects
1788: - ksp - array of linear solver object
1790: Notes:
1791: The number of KSP solvers is equal to the number of integration points divided by
1792: the number of partitions. This value is halved in the case of real matrices with
1793: a region centered at the real axis.
1795: Level: advanced
1797: .seealso: EPSCISSSetSizes()
1798: @*/
1799: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1800: {
1805: PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1806: return(0);
1807: }
1809: PetscErrorCode EPSReset_CISS(EPS eps)
1810: {
1812: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1813: PetscInt i;
1816: BVDestroy(&ctx->S);
1817: BVDestroy(&ctx->V);
1818: BVDestroy(&ctx->Y);
1819: if (!ctx->usest) {
1820: for (i=0;i<ctx->num_solve_point;i++) {
1821: KSPReset(ctx->ksp[i]);
1822: }
1823: }
1824: VecScatterDestroy(&ctx->scatterin);
1825: VecDestroy(&ctx->xsub);
1826: VecDestroy(&ctx->xdup);
1827: if (ctx->pA) {
1828: MatDestroy(&ctx->pA);
1829: MatDestroy(&ctx->pB);
1830: BVDestroy(&ctx->pV);
1831: }
1832: return(0);
1833: }
1835: PetscErrorCode EPSSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,EPS eps)
1836: {
1837: PetscErrorCode ierr;
1838: PetscReal r3,r4;
1839: PetscInt i,i1,i2,i3,i4,i5,i6,i7;
1840: PetscBool b1,b2,flg;
1841: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1842: EPSCISSQuadRule quad;
1843: EPSCISSExtraction extraction;
1846: PetscOptionsHead(PetscOptionsObject,"EPS CISS Options");
1848: EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1);
1849: PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,NULL);
1850: PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,NULL);
1851: PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,NULL);
1852: PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,NULL);
1853: PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,NULL);
1854: PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,NULL);
1855: EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1);
1857: EPSCISSGetThreshold(eps,&r3,&r4);
1858: PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,NULL);
1859: PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,NULL);
1860: EPSCISSSetThreshold(eps,r3,r4);
1862: EPSCISSGetRefinement(eps,&i6,&i7);
1863: PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,NULL);
1864: PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,NULL);
1865: EPSCISSSetRefinement(eps,i6,i7);
1867: EPSCISSGetUseST(eps,&b2);
1868: PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg);
1869: if (flg) { EPSCISSSetUseST(eps,b2); }
1871: PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg);
1872: if (flg) { EPSCISSSetQuadRule(eps,quad); }
1874: PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
1875: if (flg) { EPSCISSSetExtraction(eps,extraction); }
1877: PetscOptionsTail();
1879: if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
1880: RGSetFromOptions(eps->rg); /* this is necessary here to set useconj */
1881: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1882: for (i=0;i<ctx->num_solve_point;i++) {
1883: KSPSetFromOptions(ctx->ksp[i]);
1884: }
1885: PetscSubcommSetFromOptions(ctx->subcomm);
1886: return(0);
1887: }
1889: PetscErrorCode EPSDestroy_CISS(EPS eps)
1890: {
1892: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1895: EPSCISSResetSubcomm(eps);
1896: PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1897: PetscFree(eps->data);
1898: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL);
1899: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL);
1900: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL);
1901: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL);
1902: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL);
1903: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL);
1904: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL);
1905: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL);
1906: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL);
1907: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL);
1908: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL);
1909: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL);
1910: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL);
1911: return(0);
1912: }
1914: PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1915: {
1917: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1918: PetscBool isascii;
1921: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1922: if (isascii) {
1923: 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);
1924: if (ctx->isreal) {
1925: PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n");
1926: }
1927: PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1928: PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1929: PetscViewerASCIIPrintf(viewer," extraction: %s\n",EPSCISSExtractions[ctx->extraction]);
1930: PetscViewerASCIIPrintf(viewer," quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]);
1931: if (ctx->usest) {
1932: PetscViewerASCIIPrintf(viewer," using ST for linear solves\n");
1933: } else {
1934: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1935: PetscViewerASCIIPushTab(viewer);
1936: KSPView(ctx->ksp[0],viewer);
1937: PetscViewerASCIIPopTab(viewer);
1938: }
1939: }
1940: return(0);
1941: }
1943: PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
1944: {
1946: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1947: PetscBool usest = ctx->usest;
1950: if (!((PetscObject)eps->st)->type_name) {
1951: if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
1952: if (usest) {
1953: STSetType(eps->st,STSINVERT);
1954: } else {
1955: /* we are not going to use ST, so avoid factorizing the matrix */
1956: STSetType(eps->st,STSHIFT);
1957: }
1958: }
1959: return(0);
1960: }
1962: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
1963: {
1965: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1968: PetscNewLog(eps,&ctx);
1969: eps->data = ctx;
1971: eps->useds = PETSC_TRUE;
1972: eps->categ = EPS_CATEGORY_CONTOUR;
1974: eps->ops->solve = EPSSolve_CISS;
1975: eps->ops->setup = EPSSetUp_CISS;
1976: eps->ops->setfromoptions = EPSSetFromOptions_CISS;
1977: eps->ops->destroy = EPSDestroy_CISS;
1978: eps->ops->reset = EPSReset_CISS;
1979: eps->ops->view = EPSView_CISS;
1980: eps->ops->computevectors = EPSComputeVectors_CISS;
1981: eps->ops->setdefaultst = EPSSetDefaultST_CISS;
1983: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS);
1984: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS);
1985: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS);
1986: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS);
1987: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS);
1988: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS);
1989: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS);
1990: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS);
1991: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS);
1992: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS);
1993: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS);
1994: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS);
1995: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS);
1996: /* set default values of parameters */
1997: ctx->N = 32;
1998: ctx->L = 16;
1999: ctx->M = ctx->N/4;
2000: ctx->delta = 1e-12;
2001: ctx->L_max = 64;
2002: ctx->spurious_threshold = 1e-4;
2003: ctx->usest = PETSC_TRUE;
2004: ctx->usest_set = PETSC_FALSE;
2005: ctx->isreal = PETSC_FALSE;
2006: ctx->refine_inner = 0;
2007: ctx->refine_blocksize = 0;
2008: ctx->npart = 1;
2009: ctx->quad = (EPSCISSQuadRule)0;
2010: ctx->extraction = EPS_CISS_EXTRACTION_RITZ;
2011: return(0);
2012: }