Actual source code: feast.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }