Actual source code: feast.c

slepc-3.20.2 2024-03-15
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: 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->ncv!=PETSC_DEFAULT) {
 60:     PetscCheck(eps->ncv>=eps->nev+2,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The value of ncv must be at least nev+2");
 61:   } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
 62:   if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
 63:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 20;
 64:   if (!eps->which) eps->which = EPS_ALL;
 65:   PetscCheck(eps->which==EPS_ALL && eps->inta!=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver must be used with a computational interval");
 66:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
 67:   EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);

 69:   if (!ctx->npoints) ctx->npoints = 8;

 71:   ncv = eps->ncv;
 72:   PetscCall(PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq));
 73:   PetscCall(PetscMalloc4(eps->nloc*ncv,&ctx->work1,eps->nloc*ncv,&ctx->work2,ncv*ncv,&ctx->Aq,ncv*ncv,&ctx->Bq));

 75:   PetscCall(EPSAllocateSolution(eps,0));
 76:   PetscCall(EPSSetWorkVecs(eps,2));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: static 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,ld;
 86:   PetscScalar    *pV,*pz,*X=NULL;
 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:   PetscFunctionBegin;
 96:   ncv = eps->ncv;
 97:   n   = eps->nloc;

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

112:   PetscCall(PetscMalloc1(eps->ncv,&evals));
113:   PetscCall(BVGetLeadingDimension(eps->V,&ld));
114:   PetscCall(BVGetArray(eps->V,&pV));
115:   if (ld==n) X = pV;
116:   else PetscCall(PetscMalloc1(eps->ncv*n,&X));

118:   ijob = -1;           /* first call to reverse communication interface */
119:   PetscCall(STGetNumMatrices(eps->st,&nmat));
120:   PetscCall(STGetMatrix(eps->st,0,&A));
121:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
122:   else B = NULL;
123:   PetscCall(MatCreateVecsEmpty(A,&x,&y));

125:   do {

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

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

175:   } while (ijob);

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

193:   for (i=0;i<eps->nconv;i++) eps->eigr[i] = evals[i];
194:   if (ld!=n) {
195:     for (i=0;i<eps->nconv;i++) PetscCall(PetscArraycpy(pV+i*ld,X+i*n,n));
196:     PetscCall(PetscFree(X));
197:   }
198:   PetscCall(BVRestoreArray(eps->V,&pV));
199:   PetscCall(VecDestroy(&x));
200:   PetscCall(VecDestroy(&y));
201:   PetscCall(PetscFree(evals));
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: static PetscErrorCode EPSReset_FEAST(EPS eps)
206: {
207:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;

209:   PetscFunctionBegin;
210:   PetscCall(PetscFree4(ctx->work1,ctx->work2,ctx->Aq,ctx->Bq));
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: static PetscErrorCode EPSDestroy_FEAST(EPS eps)
215: {
216:   PetscFunctionBegin;
217:   PetscCall(PetscFree(eps->data));
218:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",NULL));
219:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",NULL));
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: static PetscErrorCode EPSSetFromOptions_FEAST(EPS eps,PetscOptionItems *PetscOptionsObject)
224: {
225:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;
226:   PetscInt       n;
227:   PetscBool      flg;

229:   PetscFunctionBegin;
230:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS FEAST Options");

232:     n = ctx->npoints;
233:     PetscCall(PetscOptionsInt("-eps_feast_num_points","Number of contour integration points","EPSFEASTSetNumPoints",n,&n,&flg));
234:     if (flg) PetscCall(EPSFEASTSetNumPoints(eps,n));

236:   PetscOptionsHeadEnd();
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: static PetscErrorCode EPSView_FEAST(EPS eps,PetscViewer viewer)
241: {
242:   EPS_FEAST      *ctx = (EPS_FEAST*)eps->data;
243:   PetscBool      isascii;

245:   PetscFunctionBegin;
246:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
247:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer,"  number of contour integration points=%" PetscInt_FMT "\n",ctx->npoints));
248:   PetscFunctionReturn(PETSC_SUCCESS);
249: }

251: static PetscErrorCode EPSSetDefaultST_FEAST(EPS eps)
252: {
253:   PetscFunctionBegin;
254:   if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }

258: static PetscErrorCode EPSFEASTSetNumPoints_FEAST(EPS eps,PetscInt npoints)
259: {
260:   EPS_FEAST *ctx = (EPS_FEAST*)eps->data;

262:   PetscFunctionBegin;
263:   if (npoints == PETSC_DEFAULT) ctx->npoints = 8;
264:   else ctx->npoints = npoints;
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: /*@
269:    EPSFEASTSetNumPoints - Sets the number of contour integration points for
270:    the FEAST package.

272:    Logically Collective

274:    Input Parameters:
275: +  eps     - the eigenproblem solver context
276: -  npoints - number of contour integration points

278:    Options Database Key:
279: .  -eps_feast_num_points - Sets the number of points

281:    Level: advanced

283: .seealso: EPSFEASTGetNumPoints()
284: @*/
285: PetscErrorCode EPSFEASTSetNumPoints(EPS eps,PetscInt npoints)
286: {
287:   PetscFunctionBegin;
290:   PetscTryMethod(eps,"EPSFEASTSetNumPoints_C",(EPS,PetscInt),(eps,npoints));
291:   PetscFunctionReturn(PETSC_SUCCESS);
292: }

294: static PetscErrorCode EPSFEASTGetNumPoints_FEAST(EPS eps,PetscInt *npoints)
295: {
296:   EPS_FEAST *ctx = (EPS_FEAST*)eps->data;

298:   PetscFunctionBegin;
299:   *npoints = ctx->npoints;
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: /*@
304:    EPSFEASTGetNumPoints - Gets the number of contour integration points for
305:    the FEAST package.

307:    Not Collective

309:    Input Parameter:
310: .  eps     - the eigenproblem solver context

312:    Output Parameter:
313: .  npoints - number of contour integration points

315:    Level: advanced

317: .seealso: EPSFEASTSetNumPoints()
318: @*/
319: PetscErrorCode EPSFEASTGetNumPoints(EPS eps,PetscInt *npoints)
320: {
321:   PetscFunctionBegin;
323:   PetscAssertPointer(npoints,2);
324:   PetscUseMethod(eps,"EPSFEASTGetNumPoints_C",(EPS,PetscInt*),(eps,npoints));
325:   PetscFunctionReturn(PETSC_SUCCESS);
326: }

328: SLEPC_EXTERN PetscErrorCode EPSCreate_FEAST(EPS eps)
329: {
330:   EPS_FEAST      *ctx;

332:   PetscFunctionBegin;
333:   PetscCall(PetscNew(&ctx));
334:   eps->data = (void*)ctx;

336:   eps->categ = EPS_CATEGORY_CONTOUR;

338:   eps->ops->solve          = EPSSolve_FEAST;
339:   eps->ops->setup          = EPSSetUp_FEAST;
340:   eps->ops->setupsort      = EPSSetUpSort_Basic;
341:   eps->ops->setfromoptions = EPSSetFromOptions_FEAST;
342:   eps->ops->destroy        = EPSDestroy_FEAST;
343:   eps->ops->reset          = EPSReset_FEAST;
344:   eps->ops->view           = EPSView_FEAST;
345:   eps->ops->setdefaultst   = EPSSetDefaultST_FEAST;

347:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTSetNumPoints_C",EPSFEASTSetNumPoints_FEAST));
348:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSFEASTGetNumPoints_C",EPSFEASTGetNumPoints_FEAST));
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }