Actual source code: feast.c

slepc-3.17.2 2022-08-09
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);
 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: }