Actual source code: feast.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 package
12: */
14: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15: #include <../src/eps/impls/external/feast/feastp.h>
17: PetscErrorCode EPSSetUp_FEAST(EPS eps)
18: {
20: PetscInt ncv;
21: PetscBool issinv,flg;
22: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
23: PetscMPIInt size;
26: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
27: if (size!=1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The FEAST interface is supported for sequential runs only");
28: if (eps->ncv) {
29: if (eps->ncv<eps->nev+2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The value of ncv must be at least nev+2");
30: } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
31: if (eps->mpd) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
32: if (!eps->max_it) eps->max_it = PetscMax(300,(PetscInt)(2*eps->n/eps->ncv));
33: if (!eps->which) eps->which = EPS_ALL;
35: ncv = eps->ncv;
36: PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq);
37: PetscMalloc4(eps->nloc*ncv,&ctx->work1,eps->nloc*ncv,&ctx->work2,ncv*ncv,&ctx->Aq,ncv*ncv,&ctx->Bq);
38: PetscLogObjectMemory((PetscObject)eps,(2*eps->nloc*ncv+2*ncv*ncv)*sizeof(PetscScalar));
40: PetscObjectTypeCompareAny((PetscObject)eps->st,&issinv,STSINVERT,STCAYLEY,"");
41: if (!issinv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Shift-and-invert or Cayley ST is needed for FEAST");
43: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
45: if (eps->which!=EPS_ALL || (eps->inta==0.0 && eps->intb==0.0)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"FEAST must be used with a computational interval");
46: if (!eps->ishermitian) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"FEAST only available for symmetric/Hermitian eigenproblems");
47: if (eps->balance!=EPS_BALANCE_NONE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Balancing not supported in the FEAST interface");
48: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
49: if (eps->stopping!=EPSStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"External packages do not support user-defined stopping test");
51: if (!ctx->npoints) ctx->npoints = 8;
53: EPSAllocateSolution(eps,0);
54: PetscObjectTypeCompare((PetscObject)eps->V,BVVECS,&flg);
55: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires a BV with contiguous storage");
56: EPSSetWorkVecs(eps,1);
57: return(0);
58: }
60: PetscErrorCode EPSSolve_FEAST(EPS eps)
61: {
63: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
64: PetscBLASInt n,fpm[64],ijob,info,nev,ncv,loop;
65: PetscReal *evals,epsout;
66: PetscInt i,k,nmat;
67: PetscScalar *pV,Ze;
68: Vec v0,x,y,w = eps->work[0];
69: Mat A,B;
72: PetscBLASIntCast(eps->nev,&nev);
73: PetscBLASIntCast(eps->ncv,&ncv);
74: PetscBLASIntCast(eps->nloc,&n);
76: /* parameters */
77: FEASTinit_(fpm);
78: fpm[0] = (eps->numbermonitors>0)? 1: 0; /* runtime comments */
79: fpm[1] = ctx->npoints; /* contour points */
80: PetscBLASIntCast(eps->max_it,&fpm[3]); /* refinement loops */
81: #if !defined(PETSC_HAVE_MPIUNI)
82: PetscBLASIntCast(MPI_Comm_c2f(PetscObjectComm((PetscObject)eps)),&fpm[8]);
83: #endif
85: PetscMalloc1(eps->ncv,&evals);
86: BVGetColumn(eps->V,0,&v0);
87: VecGetArray(v0,&pV);
89: ijob = -1; /* first call to reverse communication interface */
90: STGetNumMatrices(eps->st,&nmat);
91: STGetMatrix(eps->st,0,&A);
92: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
93: else B = NULL;
94: MatCreateVecsEmpty(A,&x,&y);
96: do {
98: PetscStackCall("FEASTrci",FEASTrci_(&ijob,&n,&Ze,ctx->work1,ctx->work2,ctx->Aq,ctx->Bq,fpm,&epsout,&loop,&eps->inta,&eps->intb,&eps->ncv,evals,pV,&eps->nconv,eps->errest,&info));
100: if (ncv!=eps->ncv) SETERRQ1(PetscObjectComm((PetscObject)eps),1,"FEAST changed value of ncv to %d",ncv);
101: if (ijob == 10 || ijob == 20) {
102: /* set new quadrature point */
103: STSetShift(eps->st,-Ze);
104: } else if (ijob == 11 || ijob == 21) {
105: /* linear solve (A-sigma*B)\work2, overwrite work2 */
106: for (k=0;k<ncv;k++) {
107: VecPlaceArray(x,ctx->work2+eps->nloc*k);
108: if (ijob == 11) {
109: STMatSolve(eps->st,x,w);
110: } else {
111: STMatSolveTranspose(eps->st,x,w);
112: }
113: VecCopy(w,x);
114: VecScale(x,-1.0);
115: VecResetArray(x);
116: }
117: } else if (ijob == 30 || ijob == 40) {
118: /* multiplication A*V or B*V, result in work1 */
119: for (k=0;k<fpm[24];k++) {
120: VecPlaceArray(x,&pV[(fpm[23]+k-1)*eps->nloc]);
121: VecPlaceArray(y,&ctx->work1[(fpm[23]+k-1)*eps->nloc]);
122: if (ijob == 30) {
123: MatMult(A,x,y);
124: } else if (nmat>1) {
125: MatMult(B,x,y);
126: } else {
127: VecCopy(x,y);
128: }
129: VecResetArray(x);
130: VecResetArray(y);
131: }
132: } else if (ijob != 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Internal error in FEAST reverse comunication interface (ijob=%d)",ijob);
134: } while (ijob != 0);
136: eps->reason = EPS_CONVERGED_TOL;
137: eps->its = loop;
138: if (info!=0) {
139: if (info==1) { /* No eigenvalue has been found in the proposed search interval */
140: eps->nconv = 0;
141: } else if (info==2) { /* FEAST did not converge "yet" */
142: eps->reason = EPS_DIVERGED_ITS;
143: } else SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by FEAST (%d)",info);
144: }
146: for (i=0;i<eps->nconv;i++) eps->eigr[i] = evals[i];
148: VecRestoreArray(v0,&pV);
149: BVRestoreColumn(eps->V,0,&v0);
150: VecDestroy(&x);
151: VecDestroy(&y);
152: PetscFree(evals);
153: return(0);
154: }
156: PetscErrorCode EPSReset_FEAST(EPS eps)
157: {
159: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
162: PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq);
163: return(0);
164: }
166: PetscErrorCode EPSDestroy_FEAST(EPS eps)
167: {
171: PetscFree(eps->data);
172: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",NULL);
173: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",NULL);
174: return(0);
175: }
177: PetscErrorCode EPSSetFromOptions_FEAST(PetscOptionItems *PetscOptionsObject,EPS eps)
178: {
180: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
181: PetscInt n;
182: PetscBool flg;
185: PetscOptionsHead(PetscOptionsObject,"EPS FEAST Options");
187: n = ctx->npoints;
188: PetscOptionsInt("-eps_feast_num_points","Number of contour integration points","EPSFEASTSetNumPoints",n,&n,&flg);
189: if (flg) { EPSFEASTSetNumPoints(eps,n); }
191: PetscOptionsTail();
192: return(0);
193: }
195: PetscErrorCode EPSView_FEAST(EPS eps,PetscViewer viewer)
196: {
198: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
199: PetscBool isascii;
202: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
203: if (isascii) {
204: PetscViewerASCIIPrintf(viewer," number of contour integration points=%D\n",ctx->npoints);
205: }
206: return(0);
207: }
209: PetscErrorCode EPSSetDefaultST_FEAST(EPS eps)
210: {
214: if (!((PetscObject)eps->st)->type_name) {
215: STSetType(eps->st,STSINVERT);
216: }
217: return(0);
218: }
220: static PetscErrorCode EPSFEASTSetNumPoints_FEAST(EPS eps,PetscInt npoints)
221: {
223: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
226: if (npoints == PETSC_DEFAULT) ctx->npoints = 8;
227: else {
228: PetscBLASIntCast(npoints,&ctx->npoints);
229: }
230: return(0);
231: }
233: /*@
234: EPSFEASTSetNumPoints - Sets the number of contour integration points for
235: the FEAST package.
237: Collective on EPS
239: Input Parameters:
240: + eps - the eigenproblem solver context
241: - npoints - number of contour integration points
243: Options Database Key:
244: . -eps_feast_num_points - Sets the number of points
246: Level: advanced
248: .seealso: EPSFEASTGetNumPoints()
249: @*/
250: PetscErrorCode EPSFEASTSetNumPoints(EPS eps,PetscInt npoints)
251: {
257: PetscTryMethod(eps,"EPSFEASTSetNumPoints_C",(EPS,PetscInt),(eps,npoints));
258: return(0);
259: }
261: static PetscErrorCode EPSFEASTGetNumPoints_FEAST(EPS eps,PetscInt *npoints)
262: {
263: EPS_FEAST *ctx = (EPS_FEAST*)eps->data;
266: *npoints = ctx->npoints;
267: return(0);
268: }
270: /*@
271: EPSFEASTGetNumPoints - Gets the number of contour integration points for
272: the FEAST package.
274: Collective on EPS
276: Input Parameter:
277: . eps - the eigenproblem solver context
279: Output Parameter:
280: - npoints - number of contour integration points
282: Level: advanced
284: .seealso: EPSFEASTSetNumPoints()
285: @*/
286: PetscErrorCode EPSFEASTGetNumPoints(EPS eps,PetscInt *npoints)
287: {
293: PetscUseMethod(eps,"EPSFEASTGetNumPoints_C",(EPS,PetscInt*),(eps,npoints));
294: return(0);
295: }
297: SLEPC_EXTERN PetscErrorCode EPSCreate_FEAST(EPS eps)
298: {
299: EPS_FEAST *ctx;
303: PetscNewLog(eps,&ctx);
304: eps->data = (void*)ctx;
306: eps->categ = EPS_CATEGORY_CONTOUR;
308: eps->ops->solve = EPSSolve_FEAST;
309: eps->ops->setup = EPSSetUp_FEAST;
310: eps->ops->setfromoptions = EPSSetFromOptions_FEAST;
311: eps->ops->destroy = EPSDestroy_FEAST;
312: eps->ops->reset = EPSReset_FEAST;
313: eps->ops->view = EPSView_FEAST;
314: eps->ops->setdefaultst = EPSSetDefaultST_FEAST;
316: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",EPSFEASTSetNumPoints_FEAST);
317: PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",EPSFEASTGetNumPoints_FEAST);
318: return(0);
319: }