Actual source code: feast.c
slepc-main 2024-12-17
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: static PetscErrorCode EPSSetUp_FEAST(EPS eps)
49: {
50: PetscInt ncv;
51: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
52: PetscMPIInt size;
54: PetscFunctionBegin;
55: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
56: PetscCheck(size==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The FEAST interface is supported for sequential runs only");
57: EPSCheckHermitianDefinite(eps);
58: EPSCheckSinvertCayley(eps);
59: if (eps->nev==0) eps->nev = 1;
60: if (eps->ncv!=PETSC_DETERMINE) {
61: PetscCheck(eps->ncv>=eps->nev+2,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The value of ncv must be at least nev+2");
62: } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
63: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
64: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 20;
65: if (!eps->which) eps->which = EPS_ALL;
66: PetscCheck(eps->which==EPS_ALL && eps->inta!=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver must be used with a computational interval");
67: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
68: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
70: if (!ctx->npoints) ctx->npoints = 8;
72: ncv = eps->ncv;
73: PetscCall(PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq));
74: PetscCall(PetscMalloc4(eps->nloc*ncv,&ctx->work1,eps->nloc*ncv,&ctx->work2,ncv*ncv,&ctx->Aq,ncv*ncv,&ctx->Bq));
76: PetscCall(EPSAllocateSolution(eps,0));
77: PetscCall(EPSSetWorkVecs(eps,2));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: static PetscErrorCode EPSSolve_FEAST(EPS eps)
82: {
83: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
84: MKL_INT fpm[128],ijob,n,ncv,nconv,loop,info;
85: PetscReal *evals,epsout=0.0;
86: PetscInt i,k,nmat,ld;
87: PetscScalar *pV,*pz,*X=NULL;
88: Vec x,y,w=eps->work[0],z=eps->work[1];
89: Mat A,B;
90: #if defined(PETSC_USE_REAL_SINGLE)
91: MKL_Complex8 Ze;
92: #else
93: MKL_Complex16 Ze;
94: #endif
96: PetscFunctionBegin;
97: #if !defined(PETSC_HAVE_64BIT_BLAS_INDICES)
98: PetscCall(PetscBLASIntCast(eps->ncv,&ncv));
99: PetscCall(PetscBLASIntCast(eps->nloc,&n));
100: #else
101: ncv = eps->ncv;
102: n = eps->nloc;
103: #endif
105: /* parameters */
106: feastinit(fpm);
107: fpm[0] = (eps->numbermonitors>0)? 1: 0; /* runtime comments */
108: fpm[1] = (MKL_INT)ctx->npoints; /* contour points */
109: #if !defined(PETSC_USE_REAL_SINGLE)
110: fpm[2] = -PetscLog10Real(eps->tol); /* tolerance for trace */
111: #endif
112: fpm[3] = (MKL_INT)eps->max_it; /* refinement loops */
113: fpm[5] = 1; /* second stopping criterion */
114: #if defined(PETSC_USE_REAL_SINGLE)
115: fpm[6] = -PetscLog10Real(eps->tol); /* tolerance for trace */
116: #endif
118: PetscCall(PetscMalloc1(eps->ncv,&evals));
119: PetscCall(BVGetLeadingDimension(eps->V,&ld));
120: PetscCall(BVGetArray(eps->V,&pV));
121: if (ld==n) X = pV;
122: else PetscCall(PetscMalloc1(eps->ncv*n,&X));
124: ijob = -1; /* first call to reverse communication interface */
125: PetscCall(STGetNumMatrices(eps->st,&nmat));
126: PetscCall(STGetMatrix(eps->st,0,&A));
127: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
128: else B = NULL;
129: PetscCall(MatCreateVecsEmpty(A,&x,&y));
131: do {
133: 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 X,&nconv,eps->errest,&info);
135: PetscCheck(ncv==eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"FEAST changed value of ncv to %d",(int)ncv);
136: if (ijob == 10) {
137: /* set new quadrature point */
138: PetscCall(STSetShift(eps->st,Ze.real));
139: } else if (ijob == 20) {
140: /* use same quadrature point and factorization for transpose solve */
141: } else if (ijob == 11 || ijob == 21) {
142: /* linear solve (A-sigma*B)\work2, overwrite work2 */
143: for (k=0;k<ncv;k++) {
144: PetscCall(VecGetArray(z,&pz));
145: #if defined(PETSC_USE_COMPLEX)
146: for (i=0;i<eps->nloc;i++) pz[i] = PetscCMPLX(ctx->work2[eps->nloc*k+i].real,ctx->work2[eps->nloc*k+i].imag);
147: #else
148: for (i=0;i<eps->nloc;i++) pz[i] = ctx->work2[eps->nloc*k+i].real;
149: #endif
150: PetscCall(VecRestoreArray(z,&pz));
151: if (ijob == 11) PetscCall(STMatSolve(eps->st,z,w));
152: else PetscCall(STMatSolveHermitianTranspose(eps->st,z,w));
153: PetscCall(VecGetArray(w,&pz));
154: #if defined(PETSC_USE_COMPLEX)
155: for (i=0;i<eps->nloc;i++) {
156: ctx->work2[eps->nloc*k+i].real = PetscRealPart(pz[i]);
157: ctx->work2[eps->nloc*k+i].imag = PetscImaginaryPart(pz[i]);
158: }
159: #else
160: for (i=0;i<eps->nloc;i++) ctx->work2[eps->nloc*k+i].real = pz[i];
161: #endif
162: PetscCall(VecRestoreArray(w,&pz));
163: }
164: } else if (ijob == 30 || ijob == 40) {
165: /* multiplication A*V or B*V, result in work1 */
166: for (k=fpm[23]-1;k<fpm[23]+fpm[24]-1;k++) {
167: PetscCall(VecPlaceArray(x,&X[k*eps->nloc]));
168: PetscCall(VecPlaceArray(y,&ctx->work1[k*eps->nloc]));
169: if (ijob == 30) PetscCall(MatMult(A,x,y));
170: else if (nmat>1) PetscCall(MatMult(B,x,y));
171: else PetscCall(VecCopy(x,y));
172: PetscCall(VecResetArray(x));
173: PetscCall(VecResetArray(y));
174: }
175: } else PetscCheck(ijob==0 || ijob==-2,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Internal error in FEAST reverse communication interface (ijob=%d)",(int)ijob);
177: } while (ijob);
179: eps->reason = EPS_CONVERGED_TOL;
180: eps->its = loop;
181: eps->nconv = nconv;
182: if (info) {
183: switch (info) {
184: case 1: /* No eigenvalue has been found in the proposed search interval */
185: eps->nconv = 0;
186: break;
187: case 2: /* FEAST did not converge "yet" */
188: eps->reason = EPS_DIVERGED_ITS;
189: break;
190: default:
191: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by FEAST (%d)",(int)info);
192: }
193: }
195: for (i=0;i<eps->nconv;i++) eps->eigr[i] = evals[i];
196: if (ld!=n) {
197: for (i=0;i<eps->nconv;i++) PetscCall(PetscArraycpy(pV+i*ld,X+i*n,n));
198: PetscCall(PetscFree(X));
199: }
200: PetscCall(BVRestoreArray(eps->V,&pV));
201: PetscCall(VecDestroy(&x));
202: PetscCall(VecDestroy(&y));
203: PetscCall(PetscFree(evals));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode EPSReset_FEAST(EPS eps)
208: {
209: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
211: PetscFunctionBegin;
212: PetscCall(PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: static PetscErrorCode EPSDestroy_FEAST(EPS eps)
217: {
218: PetscFunctionBegin;
219: PetscCall(PetscFree(eps->data));
220: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",NULL));
221: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",NULL));
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: static PetscErrorCode EPSSetFromOptions_FEAST(EPS eps,PetscOptionItems *PetscOptionsObject)
226: {
227: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
228: PetscInt n;
229: PetscBool flg;
231: PetscFunctionBegin;
232: PetscOptionsHeadBegin(PetscOptionsObject,"EPS FEAST Options");
234: n = ctx->npoints;
235: PetscCall(PetscOptionsInt("-eps_feast_num_points","Number of contour integration points","EPSFEASTSetNumPoints",n,&n,&flg));
236: if (flg) PetscCall(EPSFEASTSetNumPoints(eps,n));
238: PetscOptionsHeadEnd();
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: static PetscErrorCode EPSView_FEAST(EPS eps,PetscViewer viewer)
243: {
244: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
245: PetscBool isascii;
247: PetscFunctionBegin;
248: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
249: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer," number of contour integration points=%" PetscInt_FMT "\n",ctx->npoints));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: static PetscErrorCode EPSSetDefaultST_FEAST(EPS eps)
254: {
255: PetscFunctionBegin;
256: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode EPSFEASTSetNumPoints_FEAST(EPS eps,PetscInt npoints)
261: {
262: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
264: PetscFunctionBegin;
265: if (npoints == PETSC_DEFAULT || npoints == PETSC_DECIDE) ctx->npoints = 8;
266: else ctx->npoints = npoints;
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: EPSFEASTSetNumPoints - Sets the number of contour integration points for
272: the FEAST package.
274: Logically Collective
276: Input Parameters:
277: + eps - the eigenproblem solver context
278: - npoints - number of contour integration points
280: Options Database Key:
281: . -eps_feast_num_points - Sets the number of points
283: Level: advanced
285: .seealso: EPSFEASTGetNumPoints()
286: @*/
287: PetscErrorCode EPSFEASTSetNumPoints(EPS eps,PetscInt npoints)
288: {
289: PetscFunctionBegin;
292: PetscTryMethod(eps,"EPSFEASTSetNumPoints_C",(EPS,PetscInt),(eps,npoints));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode EPSFEASTGetNumPoints_FEAST(EPS eps,PetscInt *npoints)
297: {
298: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
300: PetscFunctionBegin;
301: *npoints = ctx->npoints;
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: EPSFEASTGetNumPoints - Gets the number of contour integration points for
307: the FEAST package.
309: Not Collective
311: Input Parameter:
312: . eps - the eigenproblem solver context
314: Output Parameter:
315: . npoints - number of contour integration points
317: Level: advanced
319: .seealso: EPSFEASTSetNumPoints()
320: @*/
321: PetscErrorCode EPSFEASTGetNumPoints(EPS eps,PetscInt *npoints)
322: {
323: PetscFunctionBegin;
325: PetscAssertPointer(npoints,2);
326: PetscUseMethod(eps,"EPSFEASTGetNumPoints_C",(EPS,PetscInt*),(eps,npoints));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: SLEPC_EXTERN PetscErrorCode EPSCreate_FEAST(EPS eps)
331: {
332: EPS_FEAST *ctx;
334: PetscFunctionBegin;
335: PetscCall(PetscNew(&ctx));
336: eps->data = (void*)ctx;
338: eps->categ = EPS_CATEGORY_CONTOUR;
340: eps->ops->solve = EPSSolve_FEAST;
341: eps->ops->setup = EPSSetUp_FEAST;
342: eps->ops->setupsort = EPSSetUpSort_Basic;
343: eps->ops->setfromoptions = EPSSetFromOptions_FEAST;
344: eps->ops->destroy = EPSDestroy_FEAST;
345: eps->ops->reset = EPSReset_FEAST;
346: eps->ops->view = EPSView_FEAST;
347: eps->ops->setdefaultst = EPSSetDefaultST_FEAST;
349: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",EPSFEASTSetNumPoints_FEAST));
350: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",EPSFEASTGetNumPoints_FEAST));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }