Actual source code: feast.c
slepc-3.17.2 2022-08-09
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: This file implements a wrapper to the FEAST solver in MKL
12: */
14: #include <petscsys.h>
15: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
16: #define MKL_ILP64
17: #endif
18: #include <mkl.h>
19: #include <slepc/private/epsimpl.h>
21: #if defined(PETSC_USE_COMPLEX)
22: # if defined(PETSC_USE_REAL_SINGLE)
23: # define FEAST_RCI cfeast_hrci
24: # define SCALAR_CAST (MKL_Complex8*)
25: # else
26: # define FEAST_RCI zfeast_hrci
27: # define SCALAR_CAST (MKL_Complex16*)
28: # endif
29: #else
30: # if defined(PETSC_USE_REAL_SINGLE)
31: # define FEAST_RCI sfeast_srci
32: # else
33: # define FEAST_RCI dfeast_srci
34: # endif
35: # define SCALAR_CAST
36: #endif
38: typedef struct {
39: PetscInt npoints; /* number of contour points */
40: PetscScalar *work1,*Aq,*Bq; /* workspace */
41: #if defined(PETSC_USE_REAL_SINGLE)
42: MKL_Complex8 *work2;
43: #else
44: MKL_Complex16 *work2;
45: #endif
46: } EPS_FEAST;
48: PetscErrorCode EPSSetUp_FEAST(EPS eps)
49: {
50: PetscInt ncv;
51: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
52: PetscMPIInt size;
54: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
56: EPSCheckHermitianDefinite(eps);
57: EPSCheckSinvertCayley(eps);
58: if (eps->ncv!=PETSC_DEFAULT) {
60: } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
61: if (eps->mpd!=PETSC_DEFAULT) PetscInfo(eps,"Warning: parameter mpd ignored\n");
62: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 20;
63: if (!eps->which) eps->which = EPS_ALL;
65: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
66: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
68: if (!ctx->npoints) ctx->npoints = 8;
70: ncv = eps->ncv;
71: PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq);
72: PetscMalloc4(eps->nloc*ncv,&ctx->work1,eps->nloc*ncv,&ctx->work2,ncv*ncv,&ctx->Aq,ncv*ncv,&ctx->Bq);
73: PetscLogObjectMemory((PetscObject)eps,(2*eps->nloc*ncv+2*ncv*ncv)*sizeof(PetscScalar));
75: EPSAllocateSolution(eps,0);
76: EPSSetWorkVecs(eps,2);
77: PetscFunctionReturn(0);
78: }
80: PetscErrorCode EPSSolve_FEAST(EPS eps)
81: {
82: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
83: MKL_INT fpm[128],ijob,n,ncv,nconv,loop,info;
84: PetscReal *evals,epsout=0.0;
85: PetscInt i,k,nmat;
86: PetscScalar *pV,*pz;
87: Vec x,y,w=eps->work[0],z=eps->work[1];
88: Mat A,B;
89: #if defined(PETSC_USE_REAL_SINGLE)
90: MKL_Complex8 Ze;
91: #else
92: MKL_Complex16 Ze;
93: #endif
95: ncv = eps->ncv;
96: n = eps->nloc;
98: /* parameters */
99: feastinit(fpm);
100: fpm[0] = (eps->numbermonitors>0)? 1: 0; /* runtime comments */
101: fpm[1] = ctx->npoints; /* contour points */
102: #if !defined(PETSC_USE_REAL_SINGLE)
103: fpm[2] = -PetscLog10Real(eps->tol); /* tolerance for trace */
104: #endif
105: fpm[3] = eps->max_it; /* refinement loops */
106: fpm[5] = 1; /* second stopping criterion */
107: #if defined(PETSC_USE_REAL_SINGLE)
108: fpm[6] = -PetscLog10Real(eps->tol); /* tolerance for trace */
109: #endif
111: PetscMalloc1(eps->ncv,&evals);
112: BVGetArray(eps->V,&pV);
114: ijob = -1; /* first call to reverse communication interface */
115: STGetNumMatrices(eps->st,&nmat);
116: STGetMatrix(eps->st,0,&A);
117: if (nmat>1) STGetMatrix(eps->st,1,&B);
118: else B = NULL;
119: MatCreateVecsEmpty(A,&x,&y);
121: do {
123: FEAST_RCI(&ijob,&n,&Ze,SCALAR_CAST ctx->work1,ctx->work2,SCALAR_CAST ctx->Aq,SCALAR_CAST ctx->Bq,fpm,&epsout,&loop,&eps->inta,&eps->intb,&ncv,evals,SCALAR_CAST pV,&nconv,eps->errest,&info);
126: if (ijob == 10) {
127: /* set new quadrature point */
128: STSetShift(eps->st,Ze.real);
129: } else if (ijob == 20) {
130: /* use same quadrature point and factorization for transpose solve */
131: } else if (ijob == 11 || ijob == 21) {
132: /* linear solve (A-sigma*B)\work2, overwrite work2 */
133: for (k=0;k<ncv;k++) {
134: VecGetArray(z,&pz);
135: #if defined(PETSC_USE_COMPLEX)
136: for (i=0;i<eps->nloc;i++) pz[i] = PetscCMPLX(ctx->work2[eps->nloc*k+i].real,ctx->work2[eps->nloc*k+i].imag);
137: #else
138: for (i=0;i<eps->nloc;i++) pz[i] = ctx->work2[eps->nloc*k+i].real;
139: #endif
140: VecRestoreArray(z,&pz);
141: if (ijob == 11) STMatSolve(eps->st,z,w);
142: else {
143: VecConjugate(z);
144: STMatSolveTranspose(eps->st,z,w);
145: VecConjugate(w);
146: }
147: VecGetArray(w,&pz);
148: #if defined(PETSC_USE_COMPLEX)
149: for (i=0;i<eps->nloc;i++) {
150: ctx->work2[eps->nloc*k+i].real = PetscRealPart(pz[i]);
151: ctx->work2[eps->nloc*k+i].imag = PetscImaginaryPart(pz[i]);
152: }
153: #else
154: for (i=0;i<eps->nloc;i++) ctx->work2[eps->nloc*k+i].real = pz[i];
155: #endif
156: VecRestoreArray(w,&pz);
157: }
158: } else if (ijob == 30 || ijob == 40) {
159: /* multiplication A*V or B*V, result in work1 */
160: for (k=fpm[23]-1;k<fpm[23]+fpm[24]-1;k++) {
161: VecPlaceArray(x,&pV[k*eps->nloc]);
162: VecPlaceArray(y,&ctx->work1[k*eps->nloc]);
163: if (ijob == 30) MatMult(A,x,y);
164: else if (nmat>1) MatMult(B,x,y);
165: else VecCopy(x,y);
166: VecResetArray(x);
167: VecResetArray(y);
168: }
171: } while (ijob);
173: eps->reason = EPS_CONVERGED_TOL;
174: eps->its = loop;
175: eps->nconv = nconv;
176: if (info) {
177: switch (info) {
178: case 1: /* No eigenvalue has been found in the proposed search interval */
179: eps->nconv = 0;
180: break;
181: case 2: /* FEAST did not converge "yet" */
182: eps->reason = EPS_DIVERGED_ITS;
183: break;
184: default:
185: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by FEAST (%d)",info);
186: }
187: }
189: for (i=0;i<eps->nconv;i++) eps->eigr[i] = evals[i];
191: BVRestoreArray(eps->V,&pV);
192: VecDestroy(&x);
193: VecDestroy(&y);
194: PetscFree(evals);
195: PetscFunctionReturn(0);
196: }
198: PetscErrorCode EPSReset_FEAST(EPS eps)
199: {
200: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
202: PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq);
203: PetscFunctionReturn(0);
204: }
206: PetscErrorCode EPSDestroy_FEAST(EPS eps)
207: {
208: PetscFree(eps->data);
209: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",NULL);
210: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",NULL);
211: PetscFunctionReturn(0);
212: }
214: PetscErrorCode EPSSetFromOptions_FEAST(PetscOptionItems *PetscOptionsObject,EPS eps)
215: {
216: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
217: PetscInt n;
218: PetscBool flg;
220: PetscOptionsHead(PetscOptionsObject,"EPS FEAST Options");
222: n = ctx->npoints;
223: PetscOptionsInt("-eps_feast_num_points","Number of contour integration points","EPSFEASTSetNumPoints",n,&n,&flg);
224: if (flg) EPSFEASTSetNumPoints(eps,n);
226: PetscOptionsTail();
227: PetscFunctionReturn(0);
228: }
230: PetscErrorCode EPSView_FEAST(EPS eps,PetscViewer viewer)
231: {
232: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
233: PetscBool isascii;
235: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
236: if (isascii) PetscViewerASCIIPrintf(viewer," number of contour integration points=%" PetscInt_FMT "\n",ctx->npoints);
237: PetscFunctionReturn(0);
238: }
240: PetscErrorCode EPSSetDefaultST_FEAST(EPS eps)
241: {
242: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
243: PetscFunctionReturn(0);
244: }
246: static PetscErrorCode EPSFEASTSetNumPoints_FEAST(EPS eps,PetscInt npoints)
247: {
248: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
250: if (npoints == PETSC_DEFAULT) ctx->npoints = 8;
251: else ctx->npoints = npoints;
252: PetscFunctionReturn(0);
253: }
255: /*@
256: EPSFEASTSetNumPoints - Sets the number of contour integration points for
257: the FEAST package.
259: Collective on EPS
261: Input Parameters:
262: + eps - the eigenproblem solver context
263: - npoints - number of contour integration points
265: Options Database Key:
266: . -eps_feast_num_points - Sets the number of points
268: Level: advanced
270: .seealso: EPSFEASTGetNumPoints()
271: @*/
272: PetscErrorCode EPSFEASTSetNumPoints(EPS eps,PetscInt npoints)
273: {
276: PetscTryMethod(eps,"EPSFEASTSetNumPoints_C",(EPS,PetscInt),(eps,npoints));
277: PetscFunctionReturn(0);
278: }
280: static PetscErrorCode EPSFEASTGetNumPoints_FEAST(EPS eps,PetscInt *npoints)
281: {
282: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
284: *npoints = ctx->npoints;
285: PetscFunctionReturn(0);
286: }
288: /*@
289: EPSFEASTGetNumPoints - Gets the number of contour integration points for
290: the FEAST package.
292: Collective on EPS
294: Input Parameter:
295: . eps - the eigenproblem solver context
297: Output Parameter:
298: . npoints - number of contour integration points
300: Level: advanced
302: .seealso: EPSFEASTSetNumPoints()
303: @*/
304: PetscErrorCode EPSFEASTGetNumPoints(EPS eps,PetscInt *npoints)
305: {
308: PetscUseMethod(eps,"EPSFEASTGetNumPoints_C",(EPS,PetscInt*),(eps,npoints));
309: PetscFunctionReturn(0);
310: }
312: SLEPC_EXTERN PetscErrorCode EPSCreate_FEAST(EPS eps)
313: {
314: EPS_FEAST *ctx;
316: PetscNewLog(eps,&ctx);
317: eps->data = (void*)ctx;
319: eps->categ = EPS_CATEGORY_CONTOUR;
321: eps->ops->solve = EPSSolve_FEAST;
322: eps->ops->setup = EPSSetUp_FEAST;
323: eps->ops->setupsort = EPSSetUpSort_Basic;
324: eps->ops->setfromoptions = EPSSetFromOptions_FEAST;
325: eps->ops->destroy = EPSDestroy_FEAST;
326: eps->ops->reset = EPSReset_FEAST;
327: eps->ops->view = EPSView_FEAST;
328: eps->ops->setdefaultst = EPSSetDefaultST_FEAST;
330: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",EPSFEASTSetNumPoints_FEAST);
331: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",EPSFEASTGetNumPoints_FEAST);
332: PetscFunctionReturn(0);
333: }