Actual source code: ciss.c
slepc-3.15.0 2021-03-31
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: PetscLogEvent EPS_CISS_SVD;
38: typedef struct {
39: /* parameters */
40: PetscInt N; /* number of integration points (32) */
41: PetscInt L; /* block size (16) */
42: PetscInt M; /* moment degree (N/4 = 4) */
43: PetscReal delta; /* threshold of singular value (1e-12) */
44: PetscInt L_max; /* maximum number of columns of the source matrix V */
45: PetscReal spurious_threshold; /* discard spurious eigenpairs */
46: PetscBool isreal; /* A and B are real */
47: PetscInt npart; /* number of partitions */
48: PetscInt refine_inner;
49: PetscInt refine_blocksize;
50: /* private data */
51: PetscReal *sigma; /* threshold for numerical rank */
52: PetscInt subcomm_id;
53: PetscInt num_solve_point;
54: PetscScalar *weight;
55: PetscScalar *omega;
56: PetscScalar *pp;
57: BV V;
58: BV S;
59: BV pV;
60: BV Y;
61: Vec xsub;
62: Vec xdup;
63: KSP *ksp; /* ksp array for storing factorizations at integration points */
64: PetscBool useconj;
65: PetscReal est_eig;
66: VecScatter scatterin;
67: Mat pA,pB;
68: PetscSubcomm subcomm;
69: PetscBool usest;
70: PetscBool usest_set; /* whether the user set the usest flag or not */
71: EPSCISSQuadRule quad;
72: EPSCISSExtraction extraction;
73: } EPS_CISS;
75: /* destroy KSP objects when the number of solve points changes */
76: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSolvers(EPS eps)
77: {
79: PetscInt i;
80: EPS_CISS *ctx = (EPS_CISS*)eps->data;
83: if (ctx->ksp) {
84: for (i=0;i<ctx->num_solve_point;i++) {
85: KSPDestroy(&ctx->ksp[i]);
86: }
87: PetscFree(ctx->ksp);
88: }
89: return(0);
90: }
92: /* clean PetscSubcomm object when the number of partitions changes */
93: PETSC_STATIC_INLINE PetscErrorCode EPSCISSResetSubcomm(EPS eps)
94: {
96: EPS_CISS *ctx = (EPS_CISS*)eps->data;
99: EPSCISSResetSolvers(eps);
100: PetscSubcommDestroy(&ctx->subcomm);
101: return(0);
102: }
104: /* determine whether half of integration points can be avoided (use its conjugates);
105: depends on isreal and the center of the region */
106: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUseConj(EPS eps,PetscBool *useconj)
107: {
109: PetscScalar center;
110: PetscReal c,d;
111: PetscBool isellipse,isinterval;
112: #if defined(PETSC_USE_COMPLEX)
113: EPS_CISS *ctx = (EPS_CISS*)eps->data;
114: #endif
117: *useconj = PETSC_FALSE;
118: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
119: if (isellipse) {
120: RGEllipseGetParameters(eps->rg,¢er,NULL,NULL);
121: #if defined(PETSC_USE_COMPLEX)
122: *useconj = (ctx->isreal && PetscImaginaryPart(center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
123: #endif
124: } else {
125: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
126: if (isinterval) {
127: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
128: #if defined(PETSC_USE_COMPLEX)
129: *useconj = (ctx->isreal && c==d)? PETSC_TRUE: PETSC_FALSE;
130: #endif
131: }
132: }
133: return(0);
134: }
136: /* create PetscSubcomm object and determine num_solve_point (depends on npart, N, useconj) */
137: PETSC_STATIC_INLINE PetscErrorCode EPSCISSSetUpSubComm(EPS eps,PetscInt *num_solve_point)
138: {
140: EPS_CISS *ctx = (EPS_CISS*)eps->data;
141: PetscInt N = ctx->N;
144: PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subcomm);
145: PetscSubcommSetNumber(ctx->subcomm,ctx->npart);
146: PetscSubcommSetType(ctx->subcomm,PETSC_SUBCOMM_INTERLACED);
147: PetscLogObjectMemory((PetscObject)eps,sizeof(PetscSubcomm));
148: ctx->subcomm_id = ctx->subcomm->color;
149: EPSCISSSetUseConj(eps,&ctx->useconj);
150: if (ctx->useconj) N = N/2;
151: *num_solve_point = N / ctx->npart;
152: if (N%ctx->npart > ctx->subcomm_id) (*num_solve_point)++;
153: return(0);
154: }
156: static PetscErrorCode CISSRedundantMat(EPS eps)
157: {
159: EPS_CISS *ctx = (EPS_CISS*)eps->data;
160: Mat A,B;
161: PetscInt nmat;
164: STGetNumMatrices(eps->st,&nmat);
165: if (ctx->subcomm->n != 1) {
166: STGetMatrix(eps->st,0,&A);
167: MatDestroy(&ctx->pA);
168: MatCreateRedundantMatrix(A,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pA);
169: if (nmat>1) {
170: STGetMatrix(eps->st,1,&B);
171: MatDestroy(&ctx->pB);
172: MatCreateRedundantMatrix(B,ctx->subcomm->n,PetscSubcommChild(ctx->subcomm),MAT_INITIAL_MATRIX,&ctx->pB);
173: } else ctx->pB = NULL;
174: } else {
175: ctx->pA = NULL;
176: ctx->pB = NULL;
177: }
178: return(0);
179: }
181: static PetscErrorCode CISSScatterVec(EPS eps)
182: {
184: EPS_CISS *ctx = (EPS_CISS*)eps->data;
185: IS is1,is2;
186: Vec v0;
187: PetscInt i,j,k,mstart,mend,mlocal;
188: PetscInt *idx1,*idx2,mloc_sub;
191: VecDestroy(&ctx->xsub);
192: MatCreateVecs(ctx->pA,&ctx->xsub,NULL);
194: VecDestroy(&ctx->xdup);
195: MatGetLocalSize(ctx->pA,&mloc_sub,NULL);
196: VecCreateMPI(PetscSubcommContiguousParent(ctx->subcomm),mloc_sub,PETSC_DECIDE,&ctx->xdup);
198: VecScatterDestroy(&ctx->scatterin);
199: BVGetColumn(ctx->V,0,&v0);
200: VecGetOwnershipRange(v0,&mstart,&mend);
201: mlocal = mend - mstart;
202: PetscMalloc2(ctx->subcomm->n*mlocal,&idx1,ctx->subcomm->n*mlocal,&idx2);
203: j = 0;
204: for (k=0;k<ctx->subcomm->n;k++) {
205: for (i=mstart;i<mend;i++) {
206: idx1[j] = i;
207: idx2[j++] = i + eps->n*k;
208: }
209: }
210: ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
211: ISCreateGeneral(PetscObjectComm((PetscObject)eps),ctx->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
212: VecScatterCreate(v0,is1,ctx->xdup,is2,&ctx->scatterin);
213: ISDestroy(&is1);
214: ISDestroy(&is2);
215: PetscFree2(idx1,idx2);
216: BVRestoreColumn(ctx->V,0,&v0);
217: return(0);
218: }
220: static PetscErrorCode SetPathParameter(EPS eps)
221: {
223: EPS_CISS *ctx = (EPS_CISS*)eps->data;
224: PetscInt i,j;
225: PetscScalar center=0.0,tmp,tmp2,*omegai;
226: PetscReal theta,radius=1.0,vscale,a,b,c,d,max_w=0.0,rgscale;
227: #if defined(PETSC_USE_COMPLEX)
228: PetscReal start_ang,end_ang;
229: #endif
230: PetscBool isring=PETSC_FALSE,isellipse=PETSC_FALSE,isinterval=PETSC_FALSE;
233: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
234: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
235: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
236: RGGetScale(eps->rg,&rgscale);
237: PetscMalloc1(ctx->N+1l,&omegai);
238: RGComputeContour(eps->rg,ctx->N,ctx->omega,omegai);
239: if (isellipse) {
240: RGEllipseGetParameters(eps->rg,¢er,&radius,&vscale);
241: for (i=0;i<ctx->N;i++) {
242: #if defined(PETSC_USE_COMPLEX)
243: theta = 2.0*PETSC_PI*(i+0.5)/ctx->N;
244: ctx->pp[i] = PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
245: ctx->weight[i] = rgscale*radius*(PetscCMPLX(vscale*PetscCosReal(theta),PetscSinReal(theta)))/(PetscReal)ctx->N;
246: #else
247: theta = (PETSC_PI/ctx->N)*(i+0.5);
248: ctx->pp[i] = PetscCosReal(theta);
249: ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
250: ctx->omega[i] = rgscale*(center + radius*ctx->pp[i]);
251: #endif
252: }
253: } else if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
254: for (i=0;i<ctx->N;i++) {
255: theta = (PETSC_PI/ctx->N)*(i+0.5);
256: ctx->pp[i] = PetscCosReal(theta);
257: ctx->weight[i] = PetscCosReal((ctx->N-1)*theta)/ctx->N;
258: }
259: if (isinterval) {
260: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
261: 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");
262: for (i=0;i<ctx->N;i++) {
263: if (c==d) ctx->omega[i] = ((b-a)*(ctx->pp[i]+1.0)/2.0+a)*rgscale;
264: if (a==b) {
265: #if defined(PETSC_USE_COMPLEX)
266: ctx->omega[i] = ((d-c)*(ctx->pp[i]+1.0)/2.0+c)*rgscale*PETSC_i;
267: #else
268: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
269: #endif
270: }
271: }
272: }
273: if (isring) { /* only supported in complex scalars */
274: #if defined(PETSC_USE_COMPLEX)
275: RGRingGetParameters(eps->rg,¢er,&radius,&vscale,&start_ang,&end_ang,NULL);
276: for (i=0;i<ctx->N;i++) {
277: theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(ctx->pp[i])+1.0))*PETSC_PI;
278: ctx->omega[i] = rgscale*(center + radius*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta)));
279: }
280: #endif
281: }
282: } else {
283: if (isinterval) {
284: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
285: center = rgscale*((b+a)/2.0+(d+c)/2.0*PETSC_PI);
286: radius = PetscSqrtReal(PetscPowRealInt(rgscale*(b-a)/2.0,2)+PetscPowRealInt(rgscale*(d-c)/2.0,2));
287: } else if (isring) {
288: RGRingGetParameters(eps->rg,¢er,&radius,NULL,NULL,NULL,NULL);
289: center *= rgscale;
290: radius *= rgscale;
291: }
292: for (i=0;i<ctx->N;i++) {
293: ctx->pp[i] = (ctx->omega[i]-center)/radius;
294: tmp = 1; tmp2 = 1;
295: for (j=0;j<ctx->N;j++) {
296: tmp *= ctx->omega[j];
297: if (i != j) tmp2 *= ctx->omega[j]-ctx->omega[i];
298: }
299: ctx->weight[i] = tmp/tmp2;
300: max_w = PetscMax(PetscAbsScalar(ctx->weight[i]),max_w);
301: }
302: for (i=0;i<ctx->N;i++) ctx->weight[i] /= (PetscScalar)max_w;
303: }
304: PetscFree(omegai);
305: return(0);
306: }
308: static PetscErrorCode CISSVecSetRandom(BV V,PetscInt i0,PetscInt i1)
309: {
311: PetscInt i,j,nlocal;
312: PetscScalar *vdata;
313: Vec x;
316: BVGetSizes(V,&nlocal,NULL,NULL);
317: for (i=i0;i<i1;i++) {
318: BVSetRandomColumn(V,i);
319: BVGetColumn(V,i,&x);
320: VecGetArray(x,&vdata);
321: for (j=0;j<nlocal;j++) {
322: vdata[j] = PetscRealPart(vdata[j]);
323: if (PetscRealPart(vdata[j]) < 0.5) vdata[j] = -1.0;
324: else vdata[j] = 1.0;
325: }
326: VecRestoreArray(x,&vdata);
327: BVRestoreColumn(V,i,&x);
328: }
329: return(0);
330: }
332: static PetscErrorCode VecScatterVecs(EPS eps,BV Vin,PetscInt n)
333: {
334: PetscErrorCode ierr;
335: EPS_CISS *ctx = (EPS_CISS*)eps->data;
336: PetscInt i;
337: Vec vi,pvi;
338: const PetscScalar *array;
341: for (i=0;i<n;i++) {
342: BVGetColumn(Vin,i,&vi);
343: VecScatterBegin(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
344: VecScatterEnd(ctx->scatterin,vi,ctx->xdup,INSERT_VALUES,SCATTER_FORWARD);
345: BVRestoreColumn(Vin,i,&vi);
346: VecGetArrayRead(ctx->xdup,&array);
347: VecPlaceArray(ctx->xsub,array);
348: BVGetColumn(ctx->pV,i,&pvi);
349: VecCopy(ctx->xsub,pvi);
350: BVRestoreColumn(ctx->pV,i,&pvi);
351: VecResetArray(ctx->xsub);
352: VecRestoreArrayRead(ctx->xdup,&array);
353: }
354: return(0);
355: }
357: static PetscErrorCode SolveLinearSystem(EPS eps,Mat A,Mat B,BV V,PetscInt L_start,PetscInt L_end,PetscBool initksp)
358: {
360: EPS_CISS *ctx = (EPS_CISS*)eps->data;
361: PetscInt i,p_id;
362: Mat Fz,kspMat,MV,BMV=NULL,MC;
363: KSP ksp;
364: const char *prefix;
367: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
368: if (ctx->usest) {
369: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Fz);
370: }
371: BVSetActiveColumns(V,L_start,L_end);
372: BVGetMat(V,&MV);
373: if (B) {
374: MatProductCreate(B,MV,NULL,&BMV);
375: MatProductSetType(BMV,MATPRODUCT_AB);
376: MatProductSetFromOptions(BMV);
377: MatProductSymbolic(BMV);
378: }
379: for (i=0;i<ctx->num_solve_point;i++) {
380: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
381: if (!ctx->usest && initksp) {
382: MatDuplicate(A,MAT_COPY_VALUES,&kspMat);
383: if (B) {
384: MatAXPY(kspMat,-ctx->omega[p_id],B,UNKNOWN_NONZERO_PATTERN);
385: } else {
386: MatShift(kspMat,-ctx->omega[p_id]);
387: }
388: KSPSetOperators(ctx->ksp[i],kspMat,kspMat);
389: /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS) */
390: KSPGetOptionsPrefix(ctx->ksp[i],&prefix);
391: MatSetOptionsPrefix(kspMat,prefix);
392: MatDestroy(&kspMat);
393: } else if (ctx->usest) {
394: STSetShift(eps->st,ctx->omega[p_id]);
395: STGetKSP(eps->st,&ksp);
396: }
397: BVSetActiveColumns(ctx->Y,i*ctx->L_max+L_start,i*ctx->L_max+L_end);
398: BVGetMat(ctx->Y,&MC);
399: if (B) {
400: MatProductNumeric(BMV);
401: if (ctx->usest) {
402: KSPMatSolve(ksp,BMV,MC);
403: } else {
404: KSPMatSolve(ctx->ksp[i],BMV,MC);
405: }
406: } else {
407: if (ctx->usest) {
408: KSPMatSolve(ksp,MV,MC);
409: } else {
410: KSPMatSolve(ctx->ksp[i],MV,MC);
411: }
412: }
413: if (ctx->usest && i<ctx->num_solve_point-1) { KSPReset(ksp); }
414: BVRestoreMat(ctx->Y,&MC);
415: }
416: MatDestroy(&BMV);
417: BVRestoreMat(V,&MV);
418: if (ctx->usest) { MatDestroy(&Fz); }
419: return(0);
420: }
422: #if defined(PETSC_USE_COMPLEX)
423: static PetscErrorCode EstimateNumberEigs(EPS eps,PetscInt *L_add)
424: {
426: EPS_CISS *ctx = (EPS_CISS*)eps->data;
427: PetscInt i,j,p_id;
428: PetscScalar tmp,m = 1,sum = 0.0;
429: PetscReal eta;
430: Vec v,vtemp,vj;
433: BVCreateVec(ctx->Y,&v);
434: BVCreateVec(ctx->V,&vtemp);
435: for (j=0;j<ctx->L;j++) {
436: VecSet(v,0);
437: for (i=0;i<ctx->num_solve_point; i++) {
438: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
439: BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
440: BVMultVec(ctx->Y,ctx->weight[p_id],1,v,&m);
441: }
442: BVGetColumn(ctx->V,j,&vj);
443: if (ctx->pA) {
444: VecSet(vtemp,0);
445: VecScatterBegin(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
446: VecScatterEnd(ctx->scatterin,v,vtemp,ADD_VALUES,SCATTER_REVERSE);
447: VecDot(vj,vtemp,&tmp);
448: } else {
449: VecDot(vj,v,&tmp);
450: }
451: BVRestoreColumn(ctx->V,j,&vj);
452: if (ctx->useconj) sum += PetscRealPart(tmp)*2;
453: else sum += tmp;
454: }
455: ctx->est_eig = PetscAbsScalar(sum/(PetscReal)ctx->L);
456: eta = PetscPowReal(10.0,-PetscLog10Real(eps->tol)/ctx->N);
457: PetscInfo1(eps,"Estimation_#Eig %f\n",(double)ctx->est_eig);
458: *L_add = (PetscInt)PetscCeilReal((ctx->est_eig*eta)/ctx->M) - ctx->L;
459: if (*L_add < 0) *L_add = 0;
460: if (*L_add>ctx->L_max-ctx->L) {
461: PetscInfo(eps,"Number of eigenvalues around the contour path may be too large\n");
462: *L_add = ctx->L_max-ctx->L;
463: }
464: VecDestroy(&v);
465: VecDestroy(&vtemp);
466: return(0);
467: }
468: #endif
470: static PetscErrorCode CalcMu(EPS eps,PetscScalar *Mu)
471: {
472: PetscErrorCode ierr;
473: PetscMPIInt sub_size,len;
474: PetscInt i,j,k,s;
475: PetscScalar *temp,*temp2,*ppk,alp;
476: const PetscScalar *pM;
477: EPS_CISS *ctx = (EPS_CISS*)eps->data;
478: Mat M;
481: MPI_Comm_size(PetscSubcommChild(ctx->subcomm),&sub_size);CHKERRMPI(ierr);
482: PetscMalloc3(ctx->num_solve_point*ctx->L*(ctx->L+1),&temp,2*ctx->M*ctx->L*ctx->L,&temp2,ctx->num_solve_point,&ppk);
483: MatCreateSeqDense(PETSC_COMM_SELF,ctx->L,ctx->L_max*ctx->num_solve_point,NULL,&M);
484: for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] = 0;
485: BVSetActiveColumns(ctx->Y,0,ctx->L_max*ctx->num_solve_point);
486: if (ctx->pA) {
487: BVSetActiveColumns(ctx->pV,0,ctx->L);
488: BVDot(ctx->Y,ctx->pV,M);
489: } else {
490: BVSetActiveColumns(ctx->V,0,ctx->L);
491: BVDot(ctx->Y,ctx->V,M);
492: }
493: MatDenseGetArrayRead(M,&pM);
494: for (i=0;i<ctx->num_solve_point;i++) {
495: for (j=0;j<ctx->L;j++) {
496: for (k=0;k<ctx->L;k++) {
497: temp[k+j*ctx->L+i*ctx->L*ctx->L]=pM[k+j*ctx->L+i*ctx->L*ctx->L_max];
498: }
499: }
500: }
501: MatDenseRestoreArrayRead(M,&pM);
502: for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
503: for (k=0;k<2*ctx->M;k++) {
504: for (j=0;j<ctx->L;j++) {
505: for (i=0;i<ctx->num_solve_point;i++) {
506: alp = ppk[i]*ctx->weight[i*ctx->subcomm->n + ctx->subcomm_id];
507: for (s=0;s<ctx->L;s++) {
508: if (ctx->useconj) temp2[s+(j+k*ctx->L)*ctx->L] += PetscRealPart(alp*temp[s+(j+i*ctx->L)*ctx->L])*2;
509: else temp2[s+(j+k*ctx->L)*ctx->L] += alp*temp[s+(j+i*ctx->L)*ctx->L];
510: }
511: }
512: }
513: for (i=0;i<ctx->num_solve_point;i++)
514: ppk[i] *= ctx->pp[i*ctx->subcomm->n + ctx->subcomm_id];
515: }
516: for (i=0;i<2*ctx->M*ctx->L*ctx->L;i++) temp2[i] /= sub_size;
517: PetscMPIIntCast(2*ctx->M*ctx->L*ctx->L,&len);
518: MPIU_Allreduce(temp2,Mu,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)eps));
519: PetscFree3(temp,temp2,ppk);
520: MatDestroy(&M);
521: return(0);
522: }
524: static PetscErrorCode BlockHankel(EPS eps,PetscScalar *Mu,PetscInt s,PetscScalar *H)
525: {
526: EPS_CISS *ctx = (EPS_CISS*)eps->data;
527: PetscInt i,j,k,L=ctx->L,M=ctx->M;
530: for (k=0;k<L*M;k++)
531: for (j=0;j<M;j++)
532: for (i=0;i<L;i++)
533: H[j*L+i+k*L*M] = Mu[i+k*L+(j+s)*L*L];
534: return(0);
535: }
537: static PetscErrorCode SVD_H0(EPS eps,PetscScalar *S,PetscInt *K)
538: {
540: EPS_CISS *ctx = (EPS_CISS*)eps->data;
541: PetscInt i,ml=ctx->L*ctx->M;
542: PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
543: PetscScalar *work;
544: #if defined(PETSC_USE_COMPLEX)
545: PetscReal *rwork;
546: #endif
549: PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0);
550: PetscMalloc1(5*ml,&work);
551: #if defined(PETSC_USE_COMPLEX)
552: PetscMalloc1(5*ml,&rwork);
553: #endif
554: PetscBLASIntCast(ml,&m);
555: n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
556: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
557: #if defined(PETSC_USE_COMPLEX)
558: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
559: #else
560: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,S,&lda,ctx->sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
561: #endif
562: SlepcCheckLapackInfo("gesvd",info);
563: PetscFPTrapPop();
564: (*K) = 0;
565: for (i=0;i<ml;i++) {
566: if (ctx->sigma[i]/PetscMax(ctx->sigma[0],1)>ctx->delta) (*K)++;
567: }
568: PetscFree(work);
569: #if defined(PETSC_USE_COMPLEX)
570: PetscFree(rwork);
571: #endif
572: PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0);
573: return(0);
574: }
576: static PetscErrorCode ConstructS(EPS eps)
577: {
579: EPS_CISS *ctx = (EPS_CISS*)eps->data;
580: PetscInt i,j,k,vec_local_size,p_id;
581: Vec v,sj;
582: PetscScalar *ppk, *v_data, m = 1;
585: BVGetSizes(ctx->Y,&vec_local_size,NULL,NULL);
586: PetscMalloc1(ctx->num_solve_point,&ppk);
587: for (i=0;i<ctx->num_solve_point;i++) ppk[i] = 1;
588: BVCreateVec(ctx->Y,&v);
589: for (k=0;k<ctx->M;k++) {
590: for (j=0;j<ctx->L;j++) {
591: VecSet(v,0);
592: for (i=0;i<ctx->num_solve_point;i++) {
593: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
594: BVSetActiveColumns(ctx->Y,i*ctx->L_max+j,i*ctx->L_max+j+1);
595: BVMultVec(ctx->Y,ppk[i]*ctx->weight[p_id],1.0,v,&m);
596: }
597: if (ctx->useconj) {
598: VecGetArray(v,&v_data);
599: for (i=0;i<vec_local_size;i++) v_data[i] = PetscRealPart(v_data[i])*2;
600: VecRestoreArray(v,&v_data);
601: }
602: BVGetColumn(ctx->S,k*ctx->L+j,&sj);
603: if (ctx->pA) {
604: VecSet(sj,0);
605: VecScatterBegin(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
606: VecScatterEnd(ctx->scatterin,v,sj,ADD_VALUES,SCATTER_REVERSE);
607: } else {
608: VecCopy(v,sj);
609: }
610: BVRestoreColumn(ctx->S,k*ctx->L+j,&sj);
611: }
612: for (i=0;i<ctx->num_solve_point;i++) {
613: p_id = i*ctx->subcomm->n + ctx->subcomm_id;
614: ppk[i] *= ctx->pp[p_id];
615: }
616: }
617: PetscFree(ppk);
618: VecDestroy(&v);
619: return(0);
620: }
622: static PetscErrorCode SVD_S(BV S,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *K)
623: {
625: PetscInt i,j,k,local_size;
626: PetscMPIInt len;
627: PetscScalar *work,*temp,*B,*tempB,*s_data,*Q1,*Q2,*temp2,alpha=1,beta=0;
628: PetscBLASInt l,m,n,lda,ldu,ldvt,lwork,info,ldb,ldc;
629: #if defined(PETSC_USE_COMPLEX)
630: PetscReal *rwork;
631: #endif
634: BVGetSizes(S,&local_size,NULL,NULL);
635: BVGetArray(S,&s_data);
636: PetscMalloc7(ml*ml,&temp,ml*ml,&temp2,local_size*ml,&Q1,local_size*ml,&Q2,ml*ml,&B,ml*ml,&tempB,5*ml,&work);
637: PetscArrayzero(B,ml*ml);
638: #if defined(PETSC_USE_COMPLEX)
639: PetscMalloc1(5*ml,&rwork);
640: #endif
641: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
643: for (i=0;i<ml;i++) B[i*ml+i]=1;
645: for (k=0;k<2;k++) {
646: PetscBLASIntCast(local_size,&m);
647: PetscBLASIntCast(ml,&l);
648: n = l; lda = m; ldb = m; ldc = l;
649: if (k == 0) {
650: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,s_data,&lda,s_data,&ldb,&beta,temp,&ldc));
651: } else if ((k%2)==1) {
652: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,temp,&ldc));
653: } else {
654: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q2,&lda,Q2,&ldb,&beta,temp,&ldc));
655: }
656: PetscArrayzero(temp2,ml*ml);
657: PetscMPIIntCast(ml*ml,&len);
658: MPIU_Allreduce(temp,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S));
660: PetscBLASIntCast(ml,&m);
661: n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
662: #if defined(PETSC_USE_COMPLEX)
663: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
664: #else
665: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
666: #endif
667: SlepcCheckLapackInfo("gesvd",info);
669: PetscBLASIntCast(local_size,&l);
670: PetscBLASIntCast(ml,&n);
671: m = n; lda = l; ldb = m; ldc = l;
672: if (k==0) {
673: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,s_data,&lda,temp2,&ldb,&beta,Q1,&ldc));
674: } else if ((k%2)==1) {
675: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
676: } else {
677: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q2,&lda,temp2,&ldb,&beta,Q1,&ldc));
678: }
680: PetscBLASIntCast(ml,&l);
681: m = l; n = l; lda = l; ldb = m; ldc = l;
682: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
683: for (i=0;i<ml;i++) {
684: sigma[i] = sqrt(sigma[i]);
685: for (j=0;j<local_size;j++) {
686: if ((k%2)==1) Q2[j+i*local_size]/=sigma[i];
687: else Q1[j+i*local_size]/=sigma[i];
688: }
689: for (j=0;j<ml;j++) {
690: B[j+i*ml]=tempB[j+i*ml]*sigma[i];
691: }
692: }
693: }
695: PetscBLASIntCast(ml,&m);
696: n = m; lda = m; ldu=1; ldvt=1;
697: #if defined(PETSC_USE_COMPLEX)
698: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
699: #else
700: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
701: #endif
702: SlepcCheckLapackInfo("gesvd",info);
704: PetscBLASIntCast(local_size,&l);
705: PetscBLASIntCast(ml,&n);
706: m = n; lda = l; ldb = m; ldc = l;
707: if ((k%2)==1) {
708: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,s_data,&ldc));
709: } else {
710: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,s_data,&ldc));
711: }
713: PetscFPTrapPop();
714: BVRestoreArray(S,&s_data);
716: (*K) = 0;
717: for (i=0;i<ml;i++) {
718: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*K)++;
719: }
720: PetscFree7(temp,temp2,Q1,Q2,B,tempB,work);
721: #if defined(PETSC_USE_COMPLEX)
722: PetscFree(rwork);
723: #endif
724: return(0);
725: }
727: static PetscErrorCode isGhost(EPS eps,PetscInt ld,PetscInt nv,PetscBool *fl)
728: {
730: EPS_CISS *ctx = (EPS_CISS*)eps->data;
731: PetscInt i,j;
732: PetscScalar *pX;
733: PetscReal *tau,s1,s2,tau_max=0.0;
736: PetscMalloc1(nv,&tau);
737: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
738: DSGetArray(eps->ds,DS_MAT_X,&pX);
740: for (i=0;i<nv;i++) {
741: s1 = 0;
742: s2 = 0;
743: for (j=0;j<nv;j++) {
744: s1 += PetscAbsScalar(PetscPowScalarInt(pX[i*ld+j],2));
745: s2 += PetscPowRealInt(PetscAbsScalar(pX[i*ld+j]),2)/ctx->sigma[j];
746: }
747: tau[i] = s1/s2;
748: tau_max = PetscMax(tau_max,tau[i]);
749: }
750: DSRestoreArray(eps->ds,DS_MAT_X,&pX);
751: for (i=0;i<nv;i++) {
752: tau[i] /= tau_max;
753: }
754: for (i=0;i<nv;i++) {
755: if (tau[i]>=ctx->spurious_threshold) fl[i] = PETSC_TRUE;
756: else fl[i] = PETSC_FALSE;
757: }
758: PetscFree(tau);
759: return(0);
760: }
762: static PetscErrorCode rescale_eig(EPS eps,PetscInt nv)
763: {
765: EPS_CISS *ctx = (EPS_CISS*)eps->data;
766: PetscInt i;
767: PetscScalar center;
768: PetscReal radius,a,b,c,d,rgscale;
769: #if defined(PETSC_USE_COMPLEX)
770: PetscReal start_ang,end_ang,vscale,theta;
771: #endif
772: PetscBool isring,isellipse,isinterval;
775: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
776: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
777: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
778: RGGetScale(eps->rg,&rgscale);
779: if (isinterval) {
780: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
781: if (c==d) {
782: for (i=0;i<nv;i++) {
783: #if defined(PETSC_USE_COMPLEX)
784: eps->eigr[i] = PetscRealPart(eps->eigr[i]);
785: #else
786: eps->eigi[i] = 0;
787: #endif
788: }
789: }
790: }
791: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
792: if (isellipse) {
793: RGEllipseGetParameters(eps->rg,¢er,&radius,NULL);
794: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
795: } else if (isinterval) {
796: RGIntervalGetEndpoints(eps->rg,&a,&b,&c,&d);
797: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
798: for (i=0;i<nv;i++) {
799: if (c==d) eps->eigr[i] = ((b-a)*(eps->eigr[i]+1.0)/2.0+a)*rgscale;
800: if (a==b) {
801: #if defined(PETSC_USE_COMPLEX)
802: eps->eigr[i] = ((d-c)*(eps->eigr[i]+1.0)/2.0+c)*rgscale*PETSC_i;
803: #else
804: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
805: #endif
806: }
807: }
808: } else {
809: center = (b+a)/2.0+(d+c)/2.0*PETSC_PI;
810: radius = PetscSqrtReal(PetscPowRealInt((b-a)/2.0,2)+PetscPowRealInt((d-c)/2.0,2));
811: for (i=0;i<nv;i++) eps->eigr[i] = center + radius*eps->eigr[i];
812: }
813: } else if (isring) { /* only supported in complex scalars */
814: #if defined(PETSC_USE_COMPLEX)
815: RGRingGetParameters(eps->rg,¢er,&radius,&vscale,&start_ang,&end_ang,NULL);
816: if (ctx->quad == EPS_CISS_QUADRULE_CHEBYSHEV) {
817: for (i=0;i<nv;i++) {
818: theta = (start_ang*2.0+(end_ang-start_ang)*(PetscRealPart(eps->eigr[i])+1.0))*PETSC_PI;
819: eps->eigr[i] = rgscale*center + (rgscale*radius+PetscImaginaryPart(eps->eigr[i]))*PetscCMPLX(PetscCosReal(theta),vscale*PetscSinReal(theta));
820: }
821: } else {
822: for (i=0;i<nv;i++) eps->eigr[i] = rgscale*(center + radius*eps->eigr[i]);
823: }
824: #endif
825: }
826: }
827: return(0);
828: }
830: PetscErrorCode EPSSetUp_CISS(EPS eps)
831: {
833: EPS_CISS *ctx = (EPS_CISS*)eps->data;
834: PetscBool istrivial,isring,isellipse,isinterval,flg,useconj;
835: PetscReal c,d;
836: PetscRandom rand;
837: Mat A;
840: if (eps->ncv==PETSC_DEFAULT) {
841: eps->ncv = ctx->L_max*ctx->M;
842: if (eps->ncv>eps->n) {
843: eps->ncv = eps->n;
844: ctx->L_max = eps->ncv/ctx->M;
845: if (!ctx->L_max) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot adjust solver parameters, try setting a smaller value of M (moment size)");
846: }
847: } else {
848: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
849: ctx->L_max = eps->ncv/ctx->M;
850: if (!ctx->L_max) {
851: ctx->L_max = 1;
852: eps->ncv = ctx->L_max*ctx->M;
853: }
854: }
855: ctx->L = PetscMin(ctx->L,ctx->L_max);
856: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
857: if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
858: if (!eps->which) eps->which = EPS_ALL;
859: if (eps->which!=EPS_ALL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only computing all eigenvalues");
860: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
862: /* check region */
863: RGIsTrivial(eps->rg,&istrivial);
864: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"CISS requires a nontrivial region, e.g. -rg_type ellipse ...");
865: RGGetComplement(eps->rg,&flg);
866: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"A region with complement flag set is not allowed");
867: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
868: PetscObjectTypeCompare((PetscObject)eps->rg,RGRING,&isring);
869: PetscObjectTypeCompare((PetscObject)eps->rg,RGINTERVAL,&isinterval);
870: if (!isellipse && !isring && !isinterval) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Currently only implemented for interval, elliptic or ring regions");
871: /* if useconj has changed, then reset subcomm data */
872: EPSCISSSetUseConj(eps,&useconj);
873: if (useconj!=ctx->useconj) { EPSCISSResetSubcomm(eps); }
875: #if !defined(PETSC_USE_COMPLEX)
876: if (isring) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Ring region only supported for complex scalars");
877: #endif
878: if (isinterval) {
879: RGIntervalGetEndpoints(eps->rg,NULL,NULL,&c,&d);
880: #if !defined(PETSC_USE_COMPLEX)
881: 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");
882: #endif
883: if (!ctx->quad && c==d) ctx->quad = EPS_CISS_QUADRULE_CHEBYSHEV;
884: }
885: if (!ctx->quad) ctx->quad = EPS_CISS_QUADRULE_TRAPEZOIDAL;
887: /* create split comm */
888: if (!ctx->subcomm) { EPSCISSSetUpSubComm(eps,&ctx->num_solve_point); }
890: EPSAllocateSolution(eps,0);
891: BVGetRandomContext(eps->V,&rand); /* make sure the random context is available when duplicating */
892: if (ctx->weight) { PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma); }
893: PetscMalloc4(ctx->N,&ctx->weight,ctx->N+1,&ctx->omega,ctx->N,&ctx->pp,ctx->L_max*ctx->M,&ctx->sigma);
894: PetscLogObjectMemory((PetscObject)eps,3*ctx->N*sizeof(PetscScalar)+ctx->L_max*ctx->N*sizeof(PetscReal));
896: /* allocate basis vectors */
897: BVDestroy(&ctx->S);
898: BVDuplicateResize(eps->V,ctx->L_max*ctx->M,&ctx->S);
899: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->S);
900: BVDestroy(&ctx->V);
901: BVDuplicateResize(eps->V,ctx->L_max,&ctx->V);
902: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->V);
904: STGetMatrix(eps->st,0,&A);
905: MatIsShell(A,&flg);
906: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Matrix type shell is not supported in this solver");
908: if (!ctx->usest_set) ctx->usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
909: if (ctx->usest && ctx->npart>1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The usest flag is not supported when partitions > 1");
911: CISSRedundantMat(eps);
912: if (ctx->pA) {
913: CISSScatterVec(eps);
914: BVDestroy(&ctx->pV);
915: BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->pV);
916: BVSetSizesFromVec(ctx->pV,ctx->xsub,eps->n);
917: BVSetFromOptions(ctx->pV);
918: BVResize(ctx->pV,ctx->L_max,PETSC_FALSE);
919: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->pV);
920: }
922: EPSCheckSinvertCondition(eps,ctx->usest," (with the usest flag set)");
924: BVDestroy(&ctx->Y);
925: if (ctx->pA) {
926: BVCreate(PetscObjectComm((PetscObject)ctx->xsub),&ctx->Y);
927: BVSetSizesFromVec(ctx->Y,ctx->xsub,eps->n);
928: BVSetFromOptions(ctx->Y);
929: BVResize(ctx->Y,ctx->num_solve_point*ctx->L_max,PETSC_FALSE);
930: } else {
931: BVDuplicateResize(eps->V,ctx->num_solve_point*ctx->L_max,&ctx->Y);
932: }
933: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->Y);
935: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
936: DSSetType(eps->ds,DSGNHEP);
937: } else if (eps->isgeneralized) {
938: if (eps->ishermitian && eps->ispositive) {
939: DSSetType(eps->ds,DSGHEP);
940: } else {
941: DSSetType(eps->ds,DSGNHEP);
942: }
943: } else {
944: if (eps->ishermitian) {
945: DSSetType(eps->ds,DSHEP);
946: } else {
947: DSSetType(eps->ds,DSNHEP);
948: }
949: }
950: DSAllocate(eps->ds,eps->ncv);
952: #if !defined(PETSC_USE_COMPLEX)
953: EPSSetWorkVecs(eps,3);
954: if (!eps->ishermitian) { PetscInfo(eps,"Warning: complex eigenvalues are not calculated exactly without --with-scalar-type=complex in PETSc\n"); }
955: #else
956: EPSSetWorkVecs(eps,2);
957: #endif
958: return(0);
959: }
961: PetscErrorCode EPSSetUpSort_CISS(EPS eps)
962: {
964: SlepcSC sc;
967: /* fill sorting criterion context */
968: eps->sc->comparison = SlepcCompareSmallestReal;
969: eps->sc->comparisonctx = NULL;
970: eps->sc->map = NULL;
971: eps->sc->mapobj = NULL;
973: /* fill sorting criterion for DS */
974: DSGetSlepcSC(eps->ds,&sc);
975: sc->comparison = SlepcCompareLargestMagnitude;
976: sc->comparisonctx = NULL;
977: sc->map = NULL;
978: sc->mapobj = NULL;
979: return(0);
980: }
982: PetscErrorCode EPSSolve_CISS(EPS eps)
983: {
985: EPS_CISS *ctx = (EPS_CISS*)eps->data;
986: Mat A,B,X,M,pA,pB;
987: PetscInt i,j,ld,nmat,L_add=0,nv=0,L_base=ctx->L,inner,nlocal,*inside;
988: PetscScalar *Mu,*H0,*H1=NULL,*rr,*temp;
989: PetscReal error,max_error,norm;
990: PetscBool *fl1;
991: Vec si,si1=NULL,w[3];
992: PetscRandom rand;
993: #if defined(PETSC_USE_COMPLEX)
994: PetscBool isellipse;
995: #else
996: PetscReal normi;
997: #endif
1000: w[0] = eps->work[0];
1001: #if defined(PETSC_USE_COMPLEX)
1002: w[1] = NULL;
1003: #else
1004: w[1] = eps->work[2];
1005: #endif
1006: w[2] = eps->work[1];
1007: VecGetLocalSize(w[0],&nlocal);
1008: DSGetLeadingDimension(eps->ds,&ld);
1009: STGetNumMatrices(eps->st,&nmat);
1010: STGetMatrix(eps->st,0,&A);
1011: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
1012: else B = NULL;
1013: SetPathParameter(eps);
1014: CISSVecSetRandom(ctx->V,0,ctx->L);
1015: BVGetRandomContext(ctx->V,&rand);
1017: if (ctx->pA) {
1018: VecScatterVecs(eps,ctx->V,ctx->L);
1019: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_TRUE);
1020: } else {
1021: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_TRUE);
1022: }
1023: #if defined(PETSC_USE_COMPLEX)
1024: PetscObjectTypeCompare((PetscObject)eps->rg,RGELLIPSE,&isellipse);
1025: if (isellipse) {
1026: EstimateNumberEigs(eps,&L_add);
1027: } else {
1028: L_add = 0;
1029: }
1030: #else
1031: L_add = 0;
1032: #endif
1033: if (L_add>0) {
1034: PetscInfo2(eps,"Changing L %D -> %D by Estimate #Eig\n",ctx->L,ctx->L+L_add);
1035: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1036: if (ctx->pA) {
1037: VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1038: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1039: } else {
1040: SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1041: }
1042: ctx->L += L_add;
1043: }
1044: PetscMalloc2(ctx->L*ctx->L*ctx->M*2,&Mu,ctx->L*ctx->M*ctx->L*ctx->M,&H0);
1045: for (i=0;i<ctx->refine_blocksize;i++) {
1046: CalcMu(eps,Mu);
1047: BlockHankel(eps,Mu,0,H0);
1048: SVD_H0(eps,H0,&nv);
1049: if (ctx->sigma[0]<=ctx->delta || nv < ctx->L*ctx->M || ctx->L == ctx->L_max) break;
1050: L_add = L_base;
1051: if (ctx->L+L_add>ctx->L_max) L_add = ctx->L_max-ctx->L;
1052: PetscInfo2(eps,"Changing L %D -> %D by SVD(H0)\n",ctx->L,ctx->L+L_add);
1053: CISSVecSetRandom(ctx->V,ctx->L,ctx->L+L_add);
1054: if (ctx->pA) {
1055: VecScatterVecs(eps,ctx->V,ctx->L+L_add);
1056: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,ctx->L,ctx->L+L_add,PETSC_FALSE);
1057: } else {
1058: SolveLinearSystem(eps,A,B,ctx->V,ctx->L,ctx->L+L_add,PETSC_FALSE);
1059: }
1060: ctx->L += L_add;
1061: }
1062: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1063: PetscMalloc1(ctx->L*ctx->M*ctx->L*ctx->M,&H1);
1064: }
1066: while (eps->reason == EPS_CONVERGED_ITERATING) {
1067: eps->its++;
1068: for (inner=0;inner<=ctx->refine_inner;inner++) {
1069: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1070: CalcMu(eps,Mu);
1071: BlockHankel(eps,Mu,0,H0);
1072: SVD_H0(eps,H0,&nv);
1073: break;
1074: } else {
1075: ConstructS(eps);
1076: BVSetActiveColumns(ctx->S,0,ctx->L);
1077: BVCopy(ctx->S,ctx->V);
1078: PetscLogEventBegin(EPS_CISS_SVD,eps,0,0,0);
1079: SVD_S(ctx->S,ctx->L*ctx->M,ctx->delta,ctx->sigma,&nv);
1080: PetscLogEventEnd(EPS_CISS_SVD,eps,0,0,0);
1081: if (ctx->sigma[0]>ctx->delta && nv==ctx->L*ctx->M && inner!=ctx->refine_inner) {
1082: if (ctx->pA) {
1083: VecScatterVecs(eps,ctx->V,ctx->L);
1084: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1085: } else {
1086: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1087: }
1088: } else break;
1089: }
1090: }
1091: eps->nconv = 0;
1092: if (nv == 0) eps->reason = EPS_CONVERGED_TOL;
1093: else {
1094: DSSetDimensions(eps->ds,nv,0,0,0);
1095: DSSetState(eps->ds,DS_STATE_RAW);
1097: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1098: BlockHankel(eps,Mu,0,H0);
1099: BlockHankel(eps,Mu,1,H1);
1100: DSGetArray(eps->ds,DS_MAT_A,&temp);
1101: for (j=0;j<nv;j++) {
1102: for (i=0;i<nv;i++) {
1103: temp[i+j*ld] = H1[i+j*ctx->L*ctx->M];
1104: }
1105: }
1106: DSRestoreArray(eps->ds,DS_MAT_A,&temp);
1107: DSGetArray(eps->ds,DS_MAT_B,&temp);
1108: for (j=0;j<nv;j++) {
1109: for (i=0;i<nv;i++) {
1110: temp[i+j*ld] = H0[i+j*ctx->L*ctx->M];
1111: }
1112: }
1113: DSRestoreArray(eps->ds,DS_MAT_B,&temp);
1114: } else {
1115: BVSetActiveColumns(ctx->S,0,nv);
1116: DSGetMat(eps->ds,DS_MAT_A,&pA);
1117: MatZeroEntries(pA);
1118: BVMatProject(ctx->S,A,ctx->S,pA);
1119: DSRestoreMat(eps->ds,DS_MAT_A,&pA);
1120: if (B) {
1121: DSGetMat(eps->ds,DS_MAT_B,&pB);
1122: MatZeroEntries(pB);
1123: BVMatProject(ctx->S,B,ctx->S,pB);
1124: DSRestoreMat(eps->ds,DS_MAT_B,&pB);
1125: }
1126: }
1128: DSSolve(eps->ds,eps->eigr,eps->eigi);
1129: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
1131: PetscMalloc3(nv,&fl1,nv,&inside,nv,&rr);
1132: rescale_eig(eps,nv);
1133: isGhost(eps,ld,nv,fl1);
1134: RGCheckInside(eps->rg,nv,eps->eigr,eps->eigi,inside);
1135: for (i=0;i<nv;i++) {
1136: if (fl1[i] && inside[i]>=0) {
1137: rr[i] = 1.0;
1138: eps->nconv++;
1139: } else rr[i] = 0.0;
1140: }
1141: DSSort(eps->ds,eps->eigr,eps->eigi,rr,NULL,&eps->nconv);
1142: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
1143: rescale_eig(eps,nv);
1144: PetscFree3(fl1,inside,rr);
1145: BVSetActiveColumns(eps->V,0,nv);
1146: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1147: ConstructS(eps);
1148: BVSetActiveColumns(ctx->S,0,ctx->L);
1149: BVCopy(ctx->S,ctx->V);
1150: BVSetActiveColumns(ctx->S,0,nv);
1151: }
1152: BVCopy(ctx->S,eps->V);
1154: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
1155: DSGetMat(eps->ds,DS_MAT_X,&X);
1156: BVMultInPlace(ctx->S,X,0,eps->nconv);
1157: if (eps->ishermitian) {
1158: BVMultInPlace(eps->V,X,0,eps->nconv);
1159: }
1160: MatDestroy(&X);
1161: max_error = 0.0;
1162: for (i=0;i<eps->nconv;i++) {
1163: BVGetColumn(ctx->S,i,&si);
1164: #if !defined(PETSC_USE_COMPLEX)
1165: if (eps->eigi[i]!=0.0) { BVGetColumn(ctx->S,i+1,&si1); }
1166: #endif
1167: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,eps->eigr[i],eps->eigi[i],si,si1,w,&error);
1168: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) { /* vector is not normalized */
1169: VecNorm(si,NORM_2,&norm);
1170: #if !defined(PETSC_USE_COMPLEX)
1171: if (eps->eigi[i]!=0.0) {
1172: VecNorm(si1,NORM_2,&normi);
1173: norm = SlepcAbsEigenvalue(norm,normi);
1174: }
1175: #endif
1176: error /= norm;
1177: }
1178: (*eps->converged)(eps,eps->eigr[i],eps->eigi[i],error,&error,eps->convergedctx);
1179: BVRestoreColumn(ctx->S,i,&si);
1180: #if !defined(PETSC_USE_COMPLEX)
1181: if (eps->eigi[i]!=0.0) {
1182: BVRestoreColumn(ctx->S,i+1,&si1);
1183: i++;
1184: }
1185: #endif
1186: max_error = PetscMax(max_error,error);
1187: }
1189: if (max_error <= eps->tol) eps->reason = EPS_CONVERGED_TOL;
1190: else if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
1191: else {
1192: if (eps->nconv > ctx->L) {
1193: MatCreateSeqDense(PETSC_COMM_SELF,eps->nconv,ctx->L,NULL,&M);
1194: MatDenseGetArrayWrite(M,&temp);
1195: for (i=0;i<ctx->L*eps->nconv;i++) {
1196: PetscRandomGetValue(rand,&temp[i]);
1197: temp[i] = PetscRealPart(temp[i]);
1198: }
1199: MatDenseRestoreArrayWrite(M,&temp);
1200: BVSetActiveColumns(ctx->S,0,eps->nconv);
1201: BVMultInPlace(ctx->S,M,0,ctx->L);
1202: MatDestroy(&M);
1203: BVSetActiveColumns(ctx->S,0,ctx->L);
1204: BVCopy(ctx->S,ctx->V);
1205: }
1206: if (ctx->pA) {
1207: VecScatterVecs(eps,ctx->V,ctx->L);
1208: SolveLinearSystem(eps,ctx->pA,ctx->pB,ctx->pV,0,ctx->L,PETSC_FALSE);
1209: } else {
1210: SolveLinearSystem(eps,A,B,ctx->V,0,ctx->L,PETSC_FALSE);
1211: }
1212: }
1213: }
1214: }
1215: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1216: PetscFree(H1);
1217: }
1218: PetscFree2(Mu,H0);
1219: return(0);
1220: }
1222: PetscErrorCode EPSComputeVectors_CISS(EPS eps)
1223: {
1225: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1226: PetscInt n;
1227: Mat Z,B=NULL;
1230: if (eps->ishermitian) {
1231: if (eps->isgeneralized && !eps->ispositive) {
1232: EPSComputeVectors_Indefinite(eps);
1233: } else {
1234: EPSComputeVectors_Hermitian(eps);
1235: }
1236: return(0);
1237: }
1238: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
1239: BVSetActiveColumns(eps->V,0,n);
1241: /* right eigenvectors */
1242: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
1244: /* V = V * Z */
1245: DSGetMat(eps->ds,DS_MAT_X,&Z);
1246: BVMultInPlace(eps->V,Z,0,n);
1247: MatDestroy(&Z);
1248: BVSetActiveColumns(eps->V,0,eps->nconv);
1250: /* normalize */
1251: if (ctx->extraction == EPS_CISS_EXTRACTION_HANKEL) {
1252: if (eps->isgeneralized && eps->ishermitian && eps->ispositive) {
1253: STGetMatrix(eps->st,1,&B);
1254: BVSetMatrix(eps->V,B,PETSC_FALSE);
1255: }
1256: BVNormalize(eps->V,NULL);
1257: if (B) { BVSetMatrix(eps->V,NULL,PETSC_FALSE); }
1258: }
1259: return(0);
1260: }
1262: static PetscErrorCode EPSCISSSetSizes_CISS(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1263: {
1265: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1266: PetscInt oN,onpart;
1269: oN = ctx->N;
1270: if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) {
1271: if (ctx->N!=32) { ctx->N =32; ctx->M = ctx->N/4; }
1272: } else {
1273: if (ip<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be > 0");
1274: if (ip%2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ip argument must be an even number");
1275: if (ctx->N!=ip) { ctx->N = ip; ctx->M = ctx->N/4; }
1276: }
1277: if (bs == PETSC_DECIDE || bs == PETSC_DEFAULT) {
1278: ctx->L = 16;
1279: } else {
1280: if (bs<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bs argument must be > 0");
1281: ctx->L = bs;
1282: }
1283: if (ms == PETSC_DECIDE || ms == PETSC_DEFAULT) {
1284: ctx->M = ctx->N/4;
1285: } else {
1286: if (ms<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The ms argument must be > 0");
1287: 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");
1288: ctx->M = ms;
1289: }
1290: onpart = ctx->npart;
1291: if (npart == PETSC_DECIDE || npart == PETSC_DEFAULT) {
1292: ctx->npart = 1;
1293: } else {
1294: if (npart<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npart argument must be > 0");
1295: ctx->npart = npart;
1296: }
1297: if (bsmax == PETSC_DECIDE || bsmax == PETSC_DEFAULT) {
1298: ctx->L_max = 64;
1299: } else {
1300: if (bsmax<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The bsmax argument must be > 0");
1301: ctx->L_max = PetscMax(bsmax,ctx->L);
1302: }
1303: if (onpart != ctx->npart || oN != ctx->N || realmats != ctx->isreal) { EPSCISSResetSubcomm(eps); }
1304: ctx->isreal = realmats;
1305: eps->state = EPS_STATE_INITIAL;
1306: return(0);
1307: }
1309: /*@
1310: EPSCISSSetSizes - Sets the values of various size parameters in the CISS solver.
1312: Logically Collective on eps
1314: Input Parameters:
1315: + eps - the eigenproblem solver context
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: Options Database Keys:
1324: + -eps_ciss_integration_points - Sets the number of integration points
1325: . -eps_ciss_blocksize - Sets the block size
1326: . -eps_ciss_moments - Sets the moment size
1327: . -eps_ciss_partitions - Sets the number of partitions
1328: . -eps_ciss_maxblocksize - Sets the maximum block size
1329: - -eps_ciss_realmats - A and B are real
1331: Note:
1332: The default number of partitions is 1. This means the internal KSP object is shared
1333: among all processes of the EPS communicator. Otherwise, the communicator is split
1334: into npart communicators, so that npart KSP solves proceed simultaneously.
1336: Level: advanced
1338: .seealso: EPSCISSGetSizes()
1339: @*/
1340: PetscErrorCode EPSCISSSetSizes(EPS eps,PetscInt ip,PetscInt bs,PetscInt ms,PetscInt npart,PetscInt bsmax,PetscBool realmats)
1341: {
1352: PetscTryMethod(eps,"EPSCISSSetSizes_C",(EPS,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool),(eps,ip,bs,ms,npart,bsmax,realmats));
1353: return(0);
1354: }
1356: static PetscErrorCode EPSCISSGetSizes_CISS(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1357: {
1358: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1361: if (ip) *ip = ctx->N;
1362: if (bs) *bs = ctx->L;
1363: if (ms) *ms = ctx->M;
1364: if (npart) *npart = ctx->npart;
1365: if (bsmax) *bsmax = ctx->L_max;
1366: if (realmats) *realmats = ctx->isreal;
1367: return(0);
1368: }
1370: /*@
1371: EPSCISSGetSizes - Gets the values of various size parameters in the CISS solver.
1373: Not Collective
1375: Input Parameter:
1376: . eps - the eigenproblem solver context
1378: Output Parameters:
1379: + ip - number of integration points
1380: . bs - block size
1381: . ms - moment size
1382: . npart - number of partitions when splitting the communicator
1383: . bsmax - max block size
1384: - realmats - A and B are real
1386: Level: advanced
1388: .seealso: EPSCISSSetSizes()
1389: @*/
1390: PetscErrorCode EPSCISSGetSizes(EPS eps,PetscInt *ip,PetscInt *bs,PetscInt *ms,PetscInt *npart,PetscInt *bsmax,PetscBool *realmats)
1391: {
1396: PetscUseMethod(eps,"EPSCISSGetSizes_C",(EPS,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool*),(eps,ip,bs,ms,npart,bsmax,realmats));
1397: return(0);
1398: }
1400: static PetscErrorCode EPSCISSSetThreshold_CISS(EPS eps,PetscReal delta,PetscReal spur)
1401: {
1402: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1405: if (delta == PETSC_DEFAULT) {
1406: ctx->delta = 1e-12;
1407: } else {
1408: if (delta<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The delta argument must be > 0.0");
1409: ctx->delta = delta;
1410: }
1411: if (spur == PETSC_DEFAULT) {
1412: ctx->spurious_threshold = 1e-4;
1413: } else {
1414: if (spur<=0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The spurious threshold argument must be > 0.0");
1415: ctx->spurious_threshold = spur;
1416: }
1417: return(0);
1418: }
1420: /*@
1421: EPSCISSSetThreshold - Sets the values of various threshold parameters in
1422: the CISS solver.
1424: Logically Collective on eps
1426: Input Parameters:
1427: + eps - the eigenproblem solver context
1428: . delta - threshold for numerical rank
1429: - spur - spurious threshold (to discard spurious eigenpairs)
1431: Options Database Keys:
1432: + -eps_ciss_delta - Sets the delta
1433: - -eps_ciss_spurious_threshold - Sets the spurious threshold
1435: Level: advanced
1437: .seealso: EPSCISSGetThreshold()
1438: @*/
1439: PetscErrorCode EPSCISSSetThreshold(EPS eps,PetscReal delta,PetscReal spur)
1440: {
1447: PetscTryMethod(eps,"EPSCISSSetThreshold_C",(EPS,PetscReal,PetscReal),(eps,delta,spur));
1448: return(0);
1449: }
1451: static PetscErrorCode EPSCISSGetThreshold_CISS(EPS eps,PetscReal *delta,PetscReal *spur)
1452: {
1453: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1456: if (delta) *delta = ctx->delta;
1457: if (spur) *spur = ctx->spurious_threshold;
1458: return(0);
1459: }
1461: /*@
1462: EPSCISSGetThreshold - Gets the values of various threshold parameters
1463: in the CISS solver.
1465: Not Collective
1467: Input Parameter:
1468: . eps - the eigenproblem solver context
1470: Output Parameters:
1471: + delta - threshold for numerical rank
1472: - spur - spurious threshold (to discard spurious eigenpairs)
1474: Level: advanced
1476: .seealso: EPSCISSSetThreshold()
1477: @*/
1478: PetscErrorCode EPSCISSGetThreshold(EPS eps,PetscReal *delta,PetscReal *spur)
1479: {
1484: PetscUseMethod(eps,"EPSCISSGetThreshold_C",(EPS,PetscReal*,PetscReal*),(eps,delta,spur));
1485: return(0);
1486: }
1488: static PetscErrorCode EPSCISSSetRefinement_CISS(EPS eps,PetscInt inner,PetscInt blsize)
1489: {
1490: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1493: if (inner == PETSC_DEFAULT) {
1494: ctx->refine_inner = 0;
1495: } else {
1496: if (inner<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine inner argument must be >= 0");
1497: ctx->refine_inner = inner;
1498: }
1499: if (blsize == PETSC_DEFAULT) {
1500: ctx->refine_blocksize = 0;
1501: } else {
1502: if (blsize<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The refine blocksize argument must be >= 0");
1503: ctx->refine_blocksize = blsize;
1504: }
1505: return(0);
1506: }
1508: /*@
1509: EPSCISSSetRefinement - Sets the values of various refinement parameters
1510: in the CISS solver.
1512: Logically Collective on eps
1514: Input Parameters:
1515: + eps - the eigenproblem solver context
1516: . inner - number of iterative refinement iterations (inner loop)
1517: - blsize - number of iterative refinement iterations (blocksize loop)
1519: Options Database Keys:
1520: + -eps_ciss_refine_inner - Sets number of inner iterations
1521: - -eps_ciss_refine_blocksize - Sets number of blocksize iterations
1523: Level: advanced
1525: .seealso: EPSCISSGetRefinement()
1526: @*/
1527: PetscErrorCode EPSCISSSetRefinement(EPS eps,PetscInt inner,PetscInt blsize)
1528: {
1535: PetscTryMethod(eps,"EPSCISSSetRefinement_C",(EPS,PetscInt,PetscInt),(eps,inner,blsize));
1536: return(0);
1537: }
1539: static PetscErrorCode EPSCISSGetRefinement_CISS(EPS eps,PetscInt *inner,PetscInt *blsize)
1540: {
1541: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1544: if (inner) *inner = ctx->refine_inner;
1545: if (blsize) *blsize = ctx->refine_blocksize;
1546: return(0);
1547: }
1549: /*@
1550: EPSCISSGetRefinement - Gets the values of various refinement parameters
1551: in the CISS solver.
1553: Not Collective
1555: Input Parameter:
1556: . eps - the eigenproblem solver context
1558: Output Parameters:
1559: + inner - number of iterative refinement iterations (inner loop)
1560: - blsize - number of iterative refinement iterations (blocksize loop)
1562: Level: advanced
1564: .seealso: EPSCISSSetRefinement()
1565: @*/
1566: PetscErrorCode EPSCISSGetRefinement(EPS eps, PetscInt *inner, PetscInt *blsize)
1567: {
1572: PetscUseMethod(eps,"EPSCISSGetRefinement_C",(EPS,PetscInt*,PetscInt*),(eps,inner,blsize));
1573: return(0);
1574: }
1576: static PetscErrorCode EPSCISSSetUseST_CISS(EPS eps,PetscBool usest)
1577: {
1578: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1581: ctx->usest = usest;
1582: ctx->usest_set = PETSC_TRUE;
1583: eps->state = EPS_STATE_INITIAL;
1584: return(0);
1585: }
1587: /*@
1588: EPSCISSSetUseST - Sets a flag indicating that the CISS solver will
1589: use the ST object for the linear solves.
1591: Logically Collective on eps
1593: Input Parameters:
1594: + eps - the eigenproblem solver context
1595: - usest - boolean flag to use the ST object or not
1597: Options Database Keys:
1598: . -eps_ciss_usest <bool> - whether the ST object will be used or not
1600: Level: advanced
1602: .seealso: EPSCISSGetUseST()
1603: @*/
1604: PetscErrorCode EPSCISSSetUseST(EPS eps,PetscBool usest)
1605: {
1611: PetscTryMethod(eps,"EPSCISSSetUseST_C",(EPS,PetscBool),(eps,usest));
1612: return(0);
1613: }
1615: static PetscErrorCode EPSCISSGetUseST_CISS(EPS eps,PetscBool *usest)
1616: {
1617: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1620: *usest = ctx->usest;
1621: return(0);
1622: }
1624: /*@
1625: EPSCISSGetUseST - Gets the flag for using the ST object
1626: in the CISS solver.
1628: Not Collective
1630: Input Parameter:
1631: . eps - the eigenproblem solver context
1633: Output Parameters:
1634: . usest - boolean flag indicating if the ST object is being used
1636: Level: advanced
1638: .seealso: EPSCISSSetUseST()
1639: @*/
1640: PetscErrorCode EPSCISSGetUseST(EPS eps,PetscBool *usest)
1641: {
1647: PetscUseMethod(eps,"EPSCISSGetUseST_C",(EPS,PetscBool*),(eps,usest));
1648: return(0);
1649: }
1651: static PetscErrorCode EPSCISSSetQuadRule_CISS(EPS eps,EPSCISSQuadRule quad)
1652: {
1653: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1656: ctx->quad = quad;
1657: return(0);
1658: }
1660: /*@
1661: EPSCISSSetQuadRule - Sets the quadrature rule used in the CISS solver.
1663: Logically Collective on eps
1665: Input Parameters:
1666: + eps - the eigenproblem solver context
1667: - quad - the quadrature rule
1669: Options Database Key:
1670: . -eps_ciss_quadrule - Sets the quadrature rule (either 'trapezoidal' or
1671: 'chebyshev')
1673: Notes:
1674: By default, the trapezoidal rule is used (EPS_CISS_QUADRULE_TRAPEZOIDAL).
1676: If the 'chebyshev' option is specified (EPS_CISS_QUADRULE_CHEBYSHEV), then
1677: Chebyshev points are used as quadrature points.
1679: Level: advanced
1681: .seealso: EPSCISSGetQuadRule(), EPSCISSQuadRule
1682: @*/
1683: PetscErrorCode EPSCISSSetQuadRule(EPS eps,EPSCISSQuadRule quad)
1684: {
1690: PetscTryMethod(eps,"EPSCISSSetQuadRule_C",(EPS,EPSCISSQuadRule),(eps,quad));
1691: return(0);
1692: }
1694: static PetscErrorCode EPSCISSGetQuadRule_CISS(EPS eps,EPSCISSQuadRule *quad)
1695: {
1696: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1699: *quad = ctx->quad;
1700: return(0);
1701: }
1703: /*@
1704: EPSCISSGetQuadRule - Gets the quadrature rule used in the CISS solver.
1706: Not Collective
1708: Input Parameter:
1709: . eps - the eigenproblem solver context
1711: Output Parameters:
1712: . quad - quadrature rule
1714: Level: advanced
1716: .seealso: EPSCISSSetQuadRule() EPSCISSQuadRule
1717: @*/
1718: PetscErrorCode EPSCISSGetQuadRule(EPS eps,EPSCISSQuadRule *quad)
1719: {
1725: PetscUseMethod(eps,"EPSCISSGetQuadRule_C",(EPS,EPSCISSQuadRule*),(eps,quad));
1726: return(0);
1727: }
1729: static PetscErrorCode EPSCISSSetExtraction_CISS(EPS eps,EPSCISSExtraction extraction)
1730: {
1731: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1734: ctx->extraction = extraction;
1735: return(0);
1736: }
1738: /*@
1739: EPSCISSSetExtraction - Sets the extraction technique used in the CISS solver.
1741: Logically Collective on eps
1743: Input Parameters:
1744: + eps - the eigenproblem solver context
1745: - extraction - the extraction technique
1747: Options Database Key:
1748: . -eps_ciss_extraction - Sets the extraction technique (either 'ritz' or
1749: 'hankel')
1751: Notes:
1752: By default, the Rayleigh-Ritz extraction is used (EPS_CISS_EXTRACTION_RITZ).
1754: If the 'hankel' option is specified (EPS_CISS_EXTRACTION_HANKEL), then
1755: the Block Hankel method is used for extracting eigenpairs.
1757: Level: advanced
1759: .seealso: EPSCISSGetExtraction(), EPSCISSExtraction
1760: @*/
1761: PetscErrorCode EPSCISSSetExtraction(EPS eps,EPSCISSExtraction extraction)
1762: {
1768: PetscTryMethod(eps,"EPSCISSSetExtraction_C",(EPS,EPSCISSExtraction),(eps,extraction));
1769: return(0);
1770: }
1772: static PetscErrorCode EPSCISSGetExtraction_CISS(EPS eps,EPSCISSExtraction *extraction)
1773: {
1774: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1777: *extraction = ctx->extraction;
1778: return(0);
1779: }
1781: /*@
1782: EPSCISSGetExtraction - Gets the extraction technique used in the CISS solver.
1784: Not Collective
1786: Input Parameter:
1787: . eps - the eigenproblem solver context
1789: Output Parameters:
1790: . extraction - extraction technique
1792: Level: advanced
1794: .seealso: EPSCISSSetExtraction() EPSCISSExtraction
1795: @*/
1796: PetscErrorCode EPSCISSGetExtraction(EPS eps,EPSCISSExtraction *extraction)
1797: {
1803: PetscUseMethod(eps,"EPSCISSGetExtraction_C",(EPS,EPSCISSExtraction*),(eps,extraction));
1804: return(0);
1805: }
1807: static PetscErrorCode EPSCISSGetKSPs_CISS(EPS eps,PetscInt *nsolve,KSP **ksp)
1808: {
1810: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1811: PetscInt i;
1812: PC pc;
1815: if (!ctx->ksp) {
1816: if (!ctx->subcomm) { /* initialize subcomm first */
1817: EPSCISSSetUseConj(eps,&ctx->useconj);
1818: EPSCISSSetUpSubComm(eps,&ctx->num_solve_point);
1819: }
1820: PetscMalloc1(ctx->num_solve_point,&ctx->ksp);
1821: for (i=0;i<ctx->num_solve_point;i++) {
1822: KSPCreate(PetscSubcommChild(ctx->subcomm),&ctx->ksp[i]);
1823: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)eps,1);
1824: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)eps)->prefix);
1825: KSPAppendOptionsPrefix(ctx->ksp[i],"eps_ciss_");
1826: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ksp[i]);
1827: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)eps)->options);
1828: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1829: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1830: KSPGetPC(ctx->ksp[i],&pc);
1831: KSPSetType(ctx->ksp[i],KSPPREONLY);
1832: PCSetType(pc,PCLU);
1833: }
1834: }
1835: if (nsolve) *nsolve = ctx->num_solve_point;
1836: if (ksp) *ksp = ctx->ksp;
1837: return(0);
1838: }
1840: /*@C
1841: EPSCISSGetKSPs - Retrieve the array of linear solver objects associated with
1842: the CISS solver.
1844: Not Collective
1846: Input Parameter:
1847: . eps - the eigenproblem solver solver
1849: Output Parameters:
1850: + nsolve - number of solver objects
1851: - ksp - array of linear solver object
1853: Notes:
1854: The number of KSP solvers is equal to the number of integration points divided by
1855: the number of partitions. This value is halved in the case of real matrices with
1856: a region centered at the real axis.
1858: Level: advanced
1860: .seealso: EPSCISSSetSizes()
1861: @*/
1862: PetscErrorCode EPSCISSGetKSPs(EPS eps,PetscInt *nsolve,KSP **ksp)
1863: {
1868: PetscUseMethod(eps,"EPSCISSGetKSPs_C",(EPS,PetscInt*,KSP**),(eps,nsolve,ksp));
1869: return(0);
1870: }
1872: PetscErrorCode EPSReset_CISS(EPS eps)
1873: {
1875: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1876: PetscInt i;
1879: BVDestroy(&ctx->S);
1880: BVDestroy(&ctx->V);
1881: BVDestroy(&ctx->Y);
1882: if (!ctx->usest) {
1883: for (i=0;i<ctx->num_solve_point;i++) {
1884: KSPReset(ctx->ksp[i]);
1885: }
1886: }
1887: VecScatterDestroy(&ctx->scatterin);
1888: VecDestroy(&ctx->xsub);
1889: VecDestroy(&ctx->xdup);
1890: if (ctx->pA) {
1891: MatDestroy(&ctx->pA);
1892: MatDestroy(&ctx->pB);
1893: BVDestroy(&ctx->pV);
1894: }
1895: return(0);
1896: }
1898: PetscErrorCode EPSSetFromOptions_CISS(PetscOptionItems *PetscOptionsObject,EPS eps)
1899: {
1900: PetscErrorCode ierr;
1901: PetscReal r3,r4;
1902: PetscInt i,i1,i2,i3,i4,i5,i6,i7;
1903: PetscBool b1,b2,flg;
1904: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1905: EPSCISSQuadRule quad;
1906: EPSCISSExtraction extraction;
1909: PetscOptionsHead(PetscOptionsObject,"EPS CISS Options");
1911: EPSCISSGetSizes(eps,&i1,&i2,&i3,&i4,&i5,&b1);
1912: PetscOptionsInt("-eps_ciss_integration_points","Number of integration points","EPSCISSSetSizes",i1,&i1,NULL);
1913: PetscOptionsInt("-eps_ciss_blocksize","Block size","EPSCISSSetSizes",i2,&i2,NULL);
1914: PetscOptionsInt("-eps_ciss_moments","Moment size","EPSCISSSetSizes",i3,&i3,NULL);
1915: PetscOptionsInt("-eps_ciss_partitions","Number of partitions","EPSCISSSetSizes",i4,&i4,NULL);
1916: PetscOptionsInt("-eps_ciss_maxblocksize","Maximum block size","EPSCISSSetSizes",i5,&i5,NULL);
1917: PetscOptionsBool("-eps_ciss_realmats","True if A and B are real","EPSCISSSetSizes",b1,&b1,NULL);
1918: EPSCISSSetSizes(eps,i1,i2,i3,i4,i5,b1);
1920: EPSCISSGetThreshold(eps,&r3,&r4);
1921: PetscOptionsReal("-eps_ciss_delta","Threshold for numerical rank","EPSCISSSetThreshold",r3,&r3,NULL);
1922: PetscOptionsReal("-eps_ciss_spurious_threshold","Threshold for the spurious eigenpairs","EPSCISSSetThreshold",r4,&r4,NULL);
1923: EPSCISSSetThreshold(eps,r3,r4);
1925: EPSCISSGetRefinement(eps,&i6,&i7);
1926: PetscOptionsInt("-eps_ciss_refine_inner","Number of inner iterative refinement iterations","EPSCISSSetRefinement",i6,&i6,NULL);
1927: PetscOptionsInt("-eps_ciss_refine_blocksize","Number of blocksize iterative refinement iterations","EPSCISSSetRefinement",i7,&i7,NULL);
1928: EPSCISSSetRefinement(eps,i6,i7);
1930: EPSCISSGetUseST(eps,&b2);
1931: PetscOptionsBool("-eps_ciss_usest","Use ST for linear solves","EPSCISSSetUseST",b2,&b2,&flg);
1932: if (flg) { EPSCISSSetUseST(eps,b2); }
1934: PetscOptionsEnum("-eps_ciss_quadrule","Quadrature rule","EPSCISSSetQuadRule",EPSCISSQuadRules,(PetscEnum)ctx->quad,(PetscEnum*)&quad,&flg);
1935: if (flg) { EPSCISSSetQuadRule(eps,quad); }
1937: PetscOptionsEnum("-eps_ciss_extraction","Extraction technique","EPSCISSSetExtraction",EPSCISSExtractions,(PetscEnum)ctx->extraction,(PetscEnum*)&extraction,&flg);
1938: if (flg) { EPSCISSSetExtraction(eps,extraction); }
1940: PetscOptionsTail();
1942: if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
1943: RGSetFromOptions(eps->rg); /* this is necessary here to set useconj */
1944: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1945: for (i=0;i<ctx->num_solve_point;i++) {
1946: KSPSetFromOptions(ctx->ksp[i]);
1947: }
1948: PetscSubcommSetFromOptions(ctx->subcomm);
1949: return(0);
1950: }
1952: PetscErrorCode EPSDestroy_CISS(EPS eps)
1953: {
1955: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1958: EPSCISSResetSubcomm(eps);
1959: PetscFree4(ctx->weight,ctx->omega,ctx->pp,ctx->sigma);
1960: PetscFree(eps->data);
1961: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",NULL);
1962: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",NULL);
1963: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",NULL);
1964: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",NULL);
1965: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",NULL);
1966: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",NULL);
1967: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",NULL);
1968: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",NULL);
1969: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",NULL);
1970: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",NULL);
1971: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",NULL);
1972: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",NULL);
1973: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",NULL);
1974: return(0);
1975: }
1977: PetscErrorCode EPSView_CISS(EPS eps,PetscViewer viewer)
1978: {
1980: EPS_CISS *ctx = (EPS_CISS*)eps->data;
1981: PetscBool isascii;
1984: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1985: if (isascii) {
1986: 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);
1987: if (ctx->isreal) {
1988: PetscViewerASCIIPrintf(viewer," exploiting symmetry of integration points\n");
1989: }
1990: PetscViewerASCIIPrintf(viewer," threshold { delta: %g, spurious threshold: %g }\n",(double)ctx->delta,(double)ctx->spurious_threshold);
1991: PetscViewerASCIIPrintf(viewer," iterative refinement { inner: %D, blocksize: %D }\n",ctx->refine_inner, ctx->refine_blocksize);
1992: PetscViewerASCIIPrintf(viewer," extraction: %s\n",EPSCISSExtractions[ctx->extraction]);
1993: PetscViewerASCIIPrintf(viewer," quadrature rule: %s\n",EPSCISSQuadRules[ctx->quad]);
1994: if (ctx->usest) {
1995: PetscViewerASCIIPrintf(viewer," using ST for linear solves\n");
1996: } else {
1997: if (!ctx->ksp) { EPSCISSGetKSPs(eps,&ctx->num_solve_point,&ctx->ksp); }
1998: PetscViewerASCIIPushTab(viewer);
1999: KSPView(ctx->ksp[0],viewer);
2000: PetscViewerASCIIPopTab(viewer);
2001: }
2002: }
2003: return(0);
2004: }
2006: PetscErrorCode EPSSetDefaultST_CISS(EPS eps)
2007: {
2009: EPS_CISS *ctx = (EPS_CISS*)eps->data;
2010: PetscBool usest = ctx->usest;
2013: if (!((PetscObject)eps->st)->type_name) {
2014: if (!ctx->usest_set) usest = (ctx->npart>1)? PETSC_FALSE: PETSC_TRUE;
2015: if (usest) {
2016: STSetType(eps->st,STSINVERT);
2017: } else {
2018: /* we are not going to use ST, so avoid factorizing the matrix */
2019: STSetType(eps->st,STSHIFT);
2020: }
2021: }
2022: return(0);
2023: }
2025: SLEPC_EXTERN PetscErrorCode EPSCreate_CISS(EPS eps)
2026: {
2028: EPS_CISS *ctx = (EPS_CISS*)eps->data;
2031: PetscNewLog(eps,&ctx);
2032: eps->data = ctx;
2034: eps->useds = PETSC_TRUE;
2035: eps->categ = EPS_CATEGORY_CONTOUR;
2037: eps->ops->solve = EPSSolve_CISS;
2038: eps->ops->setup = EPSSetUp_CISS;
2039: eps->ops->setupsort = EPSSetUpSort_CISS;
2040: eps->ops->setfromoptions = EPSSetFromOptions_CISS;
2041: eps->ops->destroy = EPSDestroy_CISS;
2042: eps->ops->reset = EPSReset_CISS;
2043: eps->ops->view = EPSView_CISS;
2044: eps->ops->computevectors = EPSComputeVectors_CISS;
2045: eps->ops->setdefaultst = EPSSetDefaultST_CISS;
2047: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetSizes_C",EPSCISSSetSizes_CISS);
2048: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetSizes_C",EPSCISSGetSizes_CISS);
2049: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetThreshold_C",EPSCISSSetThreshold_CISS);
2050: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetThreshold_C",EPSCISSGetThreshold_CISS);
2051: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetRefinement_C",EPSCISSSetRefinement_CISS);
2052: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetRefinement_C",EPSCISSGetRefinement_CISS);
2053: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetUseST_C",EPSCISSSetUseST_CISS);
2054: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetUseST_C",EPSCISSGetUseST_CISS);
2055: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetQuadRule_C",EPSCISSSetQuadRule_CISS);
2056: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetQuadRule_C",EPSCISSGetQuadRule_CISS);
2057: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSSetExtraction_C",EPSCISSSetExtraction_CISS);
2058: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetExtraction_C",EPSCISSGetExtraction_CISS);
2059: PetscObjectComposeFunction((PetscObject)eps,"EPSCISSGetKSPs_C",EPSCISSGetKSPs_CISS);
2061: /* log events */
2062: PetscLogEventRegister("EPSCISS_SVD",EPS_CLASSID,&EPS_CISS_SVD);
2064: /* set default values of parameters */
2065: ctx->N = 32;
2066: ctx->L = 16;
2067: ctx->M = ctx->N/4;
2068: ctx->delta = 1e-12;
2069: ctx->L_max = 64;
2070: ctx->spurious_threshold = 1e-4;
2071: ctx->usest = PETSC_TRUE;
2072: ctx->usest_set = PETSC_FALSE;
2073: ctx->isreal = PETSC_FALSE;
2074: ctx->refine_inner = 0;
2075: ctx->refine_blocksize = 0;
2076: ctx->npart = 1;
2077: ctx->quad = (EPSCISSQuadRule)0;
2078: ctx->extraction = EPS_CISS_EXTRACTION_RITZ;
2079: return(0);
2080: }