Actual source code: feast.c

slepc-3.18.0 2022-10-01
Report Typos and Errors
  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);

 74:   EPSAllocateSolution(eps,0);
 75:   EPSSetWorkVecs(eps,2);
 76:   return 0;
 77: }

 79: PetscErrorCode EPSSolve_FEAST(EPS eps)
 80: {
 81:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;
 82:   MKL_INT        fpm[128],ijob,n,ncv,nconv,loop,info;
 83:   PetscReal      *evals,epsout=0.0;
 84:   PetscInt       i,k,nmat;
 85:   PetscScalar    *pV,*pz;
 86:   Vec            x,y,w=eps->work[0],z=eps->work[1];
 87:   Mat            A,B;
 88: #if defined(PETSC_USE_REAL_SINGLE)
 89:   MKL_Complex8   Ze;
 90: #else
 91:   MKL_Complex16  Ze;
 92: #endif

 94:   ncv = eps->ncv;
 95:   n   = eps->nloc;

 97:   /* parameters */
 98:   feastinit(fpm);
 99:   fpm[0] = (eps->numbermonitors>0)? 1: 0;   /* runtime comments */
100:   fpm[1] = ctx->npoints;                    /* contour points */
101: #if !defined(PETSC_USE_REAL_SINGLE)
102:   fpm[2] = -PetscLog10Real(eps->tol);       /* tolerance for trace */
103: #endif
104:   fpm[3] = eps->max_it;                     /* refinement loops */
105:   fpm[5] = 1;                               /* second stopping criterion */
106: #if defined(PETSC_USE_REAL_SINGLE)
107:   fpm[6] = -PetscLog10Real(eps->tol);       /* tolerance for trace */
108: #endif

110:   PetscMalloc1(eps->ncv,&evals);
111:   BVGetArray(eps->V,&pV);

113:   ijob = -1;           /* first call to reverse communication interface */
114:   STGetNumMatrices(eps->st,&nmat);
115:   STGetMatrix(eps->st,0,&A);
116:   if (nmat>1) STGetMatrix(eps->st,1,&B);
117:   else B = NULL;
118:   MatCreateVecsEmpty(A,&x,&y);

120:   do {

122:     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);

125:     if (ijob == 10) {
126:       /* set new quadrature point */
127:       STSetShift(eps->st,Ze.real);
128:     } else if (ijob == 20) {
129:       /* use same quadrature point and factorization for transpose solve */
130:     } else if (ijob == 11 || ijob == 21) {
131:       /* linear solve (A-sigma*B)\work2, overwrite work2 */
132:       for (k=0;k<ncv;k++) {
133:         VecGetArray(z,&pz);
134: #if defined(PETSC_USE_COMPLEX)
135:         for (i=0;i<eps->nloc;i++) pz[i] = PetscCMPLX(ctx->work2[eps->nloc*k+i].real,ctx->work2[eps->nloc*k+i].imag);
136: #else
137:         for (i=0;i<eps->nloc;i++) pz[i] = ctx->work2[eps->nloc*k+i].real;
138: #endif
139:         VecRestoreArray(z,&pz);
140:         if (ijob == 11) STMatSolve(eps->st,z,w);
141:         else {
142:           VecConjugate(z);
143:           STMatSolveTranspose(eps->st,z,w);
144:           VecConjugate(w);
145:         }
146:         VecGetArray(w,&pz);
147: #if defined(PETSC_USE_COMPLEX)
148:         for (i=0;i<eps->nloc;i++) {
149:           ctx->work2[eps->nloc*k+i].real = PetscRealPart(pz[i]);
150:           ctx->work2[eps->nloc*k+i].imag = PetscImaginaryPart(pz[i]);
151:         }
152: #else
153:         for (i=0;i<eps->nloc;i++) ctx->work2[eps->nloc*k+i].real = pz[i];
154: #endif
155:         VecRestoreArray(w,&pz);
156:       }
157:     } else if (ijob == 30 || ijob == 40) {
158:       /* multiplication A*V or B*V, result in work1 */
159:       for (k=fpm[23]-1;k<fpm[23]+fpm[24]-1;k++) {
160:         VecPlaceArray(x,&pV[k*eps->nloc]);
161:         VecPlaceArray(y,&ctx->work1[k*eps->nloc]);
162:         if (ijob == 30) MatMult(A,x,y);
163:         else if (nmat>1) MatMult(B,x,y);
164:         else VecCopy(x,y);
165:         VecResetArray(x);
166:         VecResetArray(y);
167:       }

170:   } while (ijob);

172:   eps->reason = EPS_CONVERGED_TOL;
173:   eps->its    = loop;
174:   eps->nconv  = nconv;
175:   if (info) {
176:     switch (info) {
177:       case 1:  /* No eigenvalue has been found in the proposed search interval */
178:         eps->nconv = 0;
179:         break;
180:       case 2:   /* FEAST did not converge "yet" */
181:         eps->reason = EPS_DIVERGED_ITS;
182:         break;
183:       default:
184:         SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by FEAST (%d)",info);
185:     }
186:   }

188:   for (i=0;i<eps->nconv;i++) eps->eigr[i] = evals[i];

190:   BVRestoreArray(eps->V,&pV);
191:   VecDestroy(&x);
192:   VecDestroy(&y);
193:   PetscFree(evals);
194:   return 0;
195: }

197: PetscErrorCode EPSReset_FEAST(EPS eps)
198: {
199:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;

201:   PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq);
202:   return 0;
203: }

205: PetscErrorCode EPSDestroy_FEAST(EPS eps)
206: {
207:   PetscFree(eps->data);
208:   PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",NULL);
209:   PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",NULL);
210:   return 0;
211: }

213: PetscErrorCode EPSSetFromOptions_FEAST(EPS eps,PetscOptionItems *PetscOptionsObject)
214: {
215:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;
216:   PetscInt       n;
217:   PetscBool      flg;

219:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS FEAST Options");

221:     n = ctx->npoints;
222:     PetscOptionsInt("-eps_feast_num_points","Number of contour integration points","EPSFEASTSetNumPoints",n,&n,&flg);
223:     if (flg) EPSFEASTSetNumPoints(eps,n);

225:   PetscOptionsHeadEnd();
226:   return 0;
227: }

229: PetscErrorCode EPSView_FEAST(EPS eps,PetscViewer viewer)
230: {
231:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;
232:   PetscBool      isascii;

234:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
235:   if (isascii) PetscViewerASCIIPrintf(viewer,"  number of contour integration points=%" PetscInt_FMT "\n",ctx->npoints);
236:   return 0;
237: }

239: PetscErrorCode EPSSetDefaultST_FEAST(EPS eps)
240: {
241:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
242:   return 0;
243: }

245: static PetscErrorCode EPSFEASTSetNumPoints_FEAST(EPS eps,PetscInt npoints)
246: {
247:   EPS_FEAST *ctx = (EPS_FEAST*)eps->data;

249:   if (npoints == PETSC_DEFAULT) ctx->npoints = 8;
250:   else ctx->npoints = npoints;
251:   return 0;
252: }

254: /*@
255:    EPSFEASTSetNumPoints - Sets the number of contour integration points for
256:    the FEAST package.

258:    Collective on EPS

260:    Input Parameters:
261: +  eps     - the eigenproblem solver context
262: -  npoints - number of contour integration points

264:    Options Database Key:
265: .  -eps_feast_num_points - Sets the number of points

267:    Level: advanced

269: .seealso: EPSFEASTGetNumPoints()
270: @*/
271: PetscErrorCode EPSFEASTSetNumPoints(EPS eps,PetscInt npoints)
272: {
275:   PetscTryMethod(eps,"EPSFEASTSetNumPoints_C",(EPS,PetscInt),(eps,npoints));
276:   return 0;
277: }

279: static PetscErrorCode EPSFEASTGetNumPoints_FEAST(EPS eps,PetscInt *npoints)
280: {
281:   EPS_FEAST *ctx = (EPS_FEAST*)eps->data;

283:   *npoints = ctx->npoints;
284:   return 0;
285: }

287: /*@
288:    EPSFEASTGetNumPoints - Gets the number of contour integration points for
289:    the FEAST package.

291:    Collective on EPS

293:    Input Parameter:
294: .  eps     - the eigenproblem solver context

296:    Output Parameter:
297: .  npoints - number of contour integration points

299:    Level: advanced

301: .seealso: EPSFEASTSetNumPoints()
302: @*/
303: PetscErrorCode EPSFEASTGetNumPoints(EPS eps,PetscInt *npoints)
304: {
307:   PetscUseMethod(eps,"EPSFEASTGetNumPoints_C",(EPS,PetscInt*),(eps,npoints));
308:   return 0;
309: }

311: SLEPC_EXTERN PetscErrorCode EPSCreate_FEAST(EPS eps)
312: {
313:   EPS_FEAST      *ctx;

315:   PetscNew(&ctx);
316:   eps->data = (void*)ctx;

318:   eps->categ = EPS_CATEGORY_CONTOUR;

320:   eps->ops->solve          = EPSSolve_FEAST;
321:   eps->ops->setup          = EPSSetUp_FEAST;
322:   eps->ops->setupsort      = EPSSetUpSort_Basic;
323:   eps->ops->setfromoptions = EPSSetFromOptions_FEAST;
324:   eps->ops->destroy        = EPSDestroy_FEAST;
325:   eps->ops->reset          = EPSReset_FEAST;
326:   eps->ops->view           = EPSView_FEAST;
327:   eps->ops->setdefaultst   = EPSSetDefaultST_FEAST;

329:   PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",EPSFEASTSetNumPoints_FEAST);
330:   PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",EPSFEASTGetNumPoints_FEAST);
331:   return 0;
332: }