Actual source code: evsl.c

slepc-3.18.1 2022-11-02
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 eigensolvers in EVSL.
 12: */

 14: #include <slepc/private/epsimpl.h>
 15: #include <evsl.h>

 17: #define PetscCallEVSL(func, ...) do {                                                   \
 18:     PetscStackPushExternal(PetscStringize(func));                                                      \
 19:     PetscErrorCode evsl_ierr_ = func(__VA_ARGS__);                                              \
 20:     PetscStackPop;                                                                             \
 22:   } while (0)

 24: typedef struct {
 25:   PetscBool         initialized;
 26:   Mat               A;           /* problem matrix */
 27:   Vec               x,y;         /* auxiliary vectors */
 28:   PetscReal         *sli;        /* slice bounds */
 29:   PetscInt          nev;         /* approximate number of wanted eigenvalues in each slice */
 30:   PetscLayout       map;         /* used to distribute slices among MPI processes */
 31:   PetscBool         estimrange;  /* the filter range was not set by the user */
 32:   /* user parameters */
 33:   PetscInt          nslices;     /* number of slices */
 34:   PetscReal         lmin,lmax;   /* numerical range (min and max eigenvalue) */
 35:   EPSEVSLDOSMethod  dos;         /* DOS method, either KPM or Lanczos */
 36:   PetscInt          nvec;        /* number of sample vectors used for DOS */
 37:   PetscInt          deg;         /* polynomial degree used for DOS (KPM only) */
 38:   PetscInt          steps;       /* number of Lanczos steps used for DOS (Lanczos only) */
 39:   PetscInt          npoints;     /* number of sample points used for DOS (Lanczos only) */
 40:   PetscInt          max_deg;     /* maximum degree allowed for the polynomial */
 41:   PetscReal         thresh;      /* threshold for accepting polynomial */
 42:   EPSEVSLDamping    damping;     /* type of damping (for polynomial and for DOS-KPM) */
 43: } EPS_EVSL;

 45: static void AMatvec_EVSL(double *xa,double *ya,void *data)
 46: {
 47:   EPS_EVSL       *ctx = (EPS_EVSL*)data;
 48:   Vec            x = ctx->x,y = ctx->y;
 49:   Mat            A = ctx->A;

 51:   PetscObjectComm((PetscObject)A),VecPlaceArray(x,(PetscScalar*)xa);
 52:   PetscObjectComm((PetscObject)A),VecPlaceArray(y,(PetscScalar*)ya);
 53:   PetscObjectComm((PetscObject)A),MatMult(A,x,y);
 54:   PetscObjectComm((PetscObject)A),VecResetArray(x);
 55:   PetscObjectComm((PetscObject)A),VecResetArray(y);
 56:   return;
 57: }

 59: PetscErrorCode EPSSetUp_EVSL(EPS eps)
 60: {
 61:   EPS_EVSL       *ctx = (EPS_EVSL*)eps->data;
 62:   PetscMPIInt    size,rank;
 63:   PetscBool      isshift;
 64:   PetscScalar    *vinit;
 65:   PetscReal      *mu,ecount,xintv[4],*xdos,*ydos;
 66:   Vec            v0;
 67:   Mat            A;
 68:   PetscRandom    rnd;

 70:   EPSCheckStandard(eps);
 71:   EPSCheckHermitian(eps);
 72:   PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift);

 75:   if (ctx->initialized) EVSLFinish();
 76:   EVSLStart();
 77:   ctx->initialized=PETSC_TRUE;

 79:   /* get number of slices per process */
 80:   MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
 81:   MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank);
 82:   if (!ctx->nslices) ctx->nslices = size;
 83:   PetscLayoutDestroy(&ctx->map);
 84:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)eps),PETSC_DECIDE,ctx->nslices,1,&ctx->map);

 86:   /* get matrix and prepare auxiliary vectors */
 87:   MatDestroy(&ctx->A);
 88:   STGetMatrix(eps->st,0,&A);
 89:   if (size==1) {
 90:     PetscObjectReference((PetscObject)A);
 91:     ctx->A = A;
 92:   } else MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&ctx->A);
 93:   SetAMatvec(eps->n,&AMatvec_EVSL,(void*)ctx);
 94:   if (!ctx->x) MatCreateVecsEmpty(ctx->A,&ctx->x,&ctx->y);
 95:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
 96:   EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE);

 98:   if (!eps->which) eps->which=EPS_ALL;

101:   /* estimate numerical range */
102:   if (ctx->estimrange || ctx->lmin == PETSC_MIN_REAL || ctx->lmax == PETSC_MAX_REAL) {
103:     MatCreateVecs(ctx->A,&v0,NULL);
104:     if (!eps->V) EPSGetBV(eps,&eps->V);
105:     BVGetRandomContext(eps->V,&rnd);
106:     VecSetRandom(v0,rnd);
107:     VecGetArray(v0,&vinit);
108:     PetscCallEVSL(LanTrbounds,50,200,eps->tol,vinit,1,&ctx->lmin,&ctx->lmax,NULL);
109:     VecRestoreArray(v0,&vinit);
110:     VecDestroy(&v0);
111:     ctx->estimrange = PETSC_TRUE;   /* estimate if called again with another matrix */
112:   }
114:   xintv[0] = eps->inta;
115:   xintv[1] = eps->intb;
116:   xintv[2] = ctx->lmin;
117:   xintv[3] = ctx->lmax;

119:   /* estimate number of eigenvalues in the interval */
120:   switch (ctx->dos) {
121:     case EPS_EVSL_DOS_KPM:
122:       PetscMalloc1(ctx->deg+1,&mu);
123:       if (!rank) PetscCallEVSL(kpmdos,ctx->deg,(int)ctx->damping,ctx->nvec,xintv,mu,&ecount);
124:       MPI_Bcast(mu,ctx->deg+1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
125:       break;
126:     case EPS_EVSL_DOS_LANCZOS:
127:       PetscMalloc2(ctx->npoints,&xdos,ctx->npoints,&ydos);
128:       if (!rank) PetscCallEVSL(LanDos,ctx->nvec,PetscMin(ctx->steps,eps->n/2),ctx->npoints,xdos,ydos,&ecount,xintv);
129:       MPI_Bcast(xdos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
130:       MPI_Bcast(ydos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
131:       break;
132:     default:
133:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid DOS method");
134:   }
135:   MPI_Bcast(&ecount,1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));

137:   PetscInfo(eps,"Estimated eigenvalue count in the interval: %g\n",ecount);
138:   eps->ncv = (PetscInt)PetscCeilReal(1.5*ecount);

140:   /* slice the spectrum */
141:   PetscFree(ctx->sli);
142:   PetscMalloc1(ctx->nslices+1,&ctx->sli);
143:   if (ctx->dos == EPS_EVSL_DOS_KPM) {
144:     PetscCallEVSL(spslicer,ctx->sli,mu,ctx->deg,xintv,ctx->nslices,10*(PetscInt)ecount);
145:     PetscFree(mu);
146:   } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
147:     spslicer2(xdos,ydos,ctx->nslices,ctx->npoints,ctx->sli);
148:     PetscFree2(xdos,ydos);
149:   }

151:   /* approximate number of eigenvalues wanted in each slice */
152:   ctx->nev = (PetscInt)(1.0 + ecount/(PetscReal)ctx->nslices) + 2;

154:   if (eps->mpd!=PETSC_DEFAULT) PetscInfo(eps,"Warning: parameter mpd ignored\n");
155:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
156:   EPSAllocateSolution(eps,0);
157:   return 0;
158: }

160: PetscErrorCode EPSSolve_EVSL(EPS eps)
161: {
162:   EPS_EVSL       *ctx = (EPS_EVSL*)eps->data;
163:   PetscInt       i,j,k=0,sl,mlan,nevout,*ind,nevmax,rstart,rend,*nevloc,*disp,N;
164:   PetscReal      *res,xintv[4],*errest;
165:   PetscScalar    *lam,*X,*Y,*vinit,*eigr;
166:   PetscMPIInt    size,rank;
167:   PetscRandom    rnd;
168:   Vec            v,w,v0,x;
169:   VecScatter     vs;
170:   IS             is;
171:   polparams      pol;

173:   MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
174:   MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank);
175:   PetscLayoutGetRange(ctx->map,&rstart,&rend);
176:   nevmax = (rend-rstart)*ctx->nev;
177:   MatCreateVecs(ctx->A,&v0,NULL);
178:   BVGetRandomContext(eps->V,&rnd);
179:   VecSetRandom(v0,rnd);
180:   VecGetArray(v0,&vinit);
181:   PetscMalloc5(size,&nevloc,size+1,&disp,nevmax,&eigr,nevmax,&errest,nevmax*eps->n,&X);
182:   mlan = PetscMin(PetscMax(5*ctx->nev,300),eps->n);
183:   for (sl=rstart; sl<rend; sl++) {
184:     xintv[0] = ctx->sli[sl];
185:     xintv[1] = ctx->sli[sl+1];
186:     xintv[2] = ctx->lmin;
187:     xintv[3] = ctx->lmax;
188:     PetscInfo(ctx->A,"Subinterval %" PetscInt_FMT ": [%.4e, %.4e]\n",sl+1,xintv[0],xintv[1]);
189:     set_pol_def(&pol);
190:     pol.max_deg    = ctx->max_deg;
191:     pol.damping    = (int)ctx->damping;
192:     pol.thresh_int = ctx->thresh;
193:     find_pol(xintv,&pol);
194:     PetscInfo(ctx->A,"Polynomial [type = %" PetscInt_FMT "], deg %" PetscInt_FMT ", bar %e gam %e\n",pol.type,pol.deg,pol.bar,pol.gam);
195:     PetscCallEVSL(ChebLanNr,xintv,mlan,eps->tol,vinit,&pol,&nevout,&lam,&Y,&res,NULL);
197:     free_pol(&pol);
198:     PetscInfo(ctx->A,"Computed %" PetscInt_FMT " eigenvalues\n",nevout);
199:     PetscMalloc1(nevout,&ind);
200:     sort_double(nevout,lam,ind);
201:     for (i=0;i<nevout;i++) {
202:       eigr[i+k]   = lam[i];
203:       errest[i+k] = res[ind[i]];
204:       PetscArraycpy(X+(i+k)*eps->n,Y+ind[i]*eps->n,eps->n);
205:     }
206:     k += nevout;
207:     if (lam) evsl_Free(lam);
208:     if (Y)   evsl_Free_device(Y);
209:     if (res) evsl_Free(res);
210:     PetscFree(ind);
211:   }
212:   VecRestoreArray(v0,&vinit);
213:   VecDestroy(&v0);

215:   /* gather eigenvalues computed by each MPI process */
216:   MPI_Allgather(&k,1,MPIU_INT,nevloc,1,MPIU_INT,PetscObjectComm((PetscObject)eps));
217:   eps->nev = nevloc[0];
218:   disp[0]  = 0;
219:   for (i=1;i<size;i++) {
220:     eps->nev += nevloc[i];
221:     disp[i]   = disp[i-1]+nevloc[i-1];
222:   }
223:   disp[size] = disp[size-1]+nevloc[size-1];
225:   MPI_Allgatherv(eigr,k,MPIU_SCALAR,eps->eigr,nevloc,disp,MPIU_SCALAR,PetscObjectComm((PetscObject)eps));
226:   MPI_Allgatherv(errest,k,MPIU_REAL,eps->errest,nevloc,disp,MPIU_REAL,PetscObjectComm((PetscObject)eps));
227:   eps->nconv  = eps->nev;
228:   eps->its    = 1;
229:   eps->reason = EPS_CONVERGED_TOL;

231:   /* scatter computed eigenvectors and store them in eps->V */
232:   BVCreateVec(eps->V,&w);
233:   for (i=0;i<size;i++) {
234:     N = (rank==i)? eps->n: 0;
235:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
236:     VecSetFromOptions(x);
237:     ISCreateStride(PETSC_COMM_SELF,N,0,1,&is);
238:     VecScatterCreate(x,is,w,is,&vs);
239:     ISDestroy(&is);
240:     for (j=disp[i];j<disp[i+1];j++) {
241:       BVGetColumn(eps->V,j,&v);
242:       if (rank==i) VecPlaceArray(x,X+(j-disp[i])*eps->n);
243:       VecScatterBegin(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
244:       VecScatterEnd(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
245:       if (rank==i) VecResetArray(x);
246:       BVRestoreColumn(eps->V,j,&v);
247:     }
248:     VecScatterDestroy(&vs);
249:     VecDestroy(&x);
250:   }
251:   VecDestroy(&w);
252:   PetscFree5(nevloc,disp,eigr,errest,X);
253:   return 0;
254: }

256: static PetscErrorCode EPSEVSLSetSlices_EVSL(EPS eps,PetscInt nslices)
257: {
258:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

260:   if (nslices == PETSC_DECIDE || nslices == PETSC_DEFAULT) nslices = 0;
262:   if (ctx->nslices != nslices) {
263:     ctx->nslices = nslices;
264:     eps->state   = EPS_STATE_INITIAL;
265:   }
266:   return 0;
267: }

269: /*@
270:    EPSEVSLSetSlices - Set the number of slices in which the interval must be
271:    subdivided.

273:    Logically Collective on eps

275:    Input Parameters:
276: +  eps     - the eigensolver context
277: -  nslices - the number of slices

279:    Options Database Key:
280: .  -eps_evsl_slices <n> - set the number of slices to n

282:    Notes:
283:    By default, one slice per MPI process is used. Depending on the number of
284:    eigenvalues, using more slices may be beneficial, but very narrow subintervals
285:    imply higher polynomial degree.

287:    Level: intermediate

289: .seealso: EPSEVSLGetSlices()
290: @*/
291: PetscErrorCode EPSEVSLSetSlices(EPS eps,PetscInt nslices)
292: {
295:   PetscTryMethod(eps,"EPSEVSLSetSlices_C",(EPS,PetscInt),(eps,nslices));
296:   return 0;
297: }

299: static PetscErrorCode EPSEVSLGetSlices_EVSL(EPS eps,PetscInt *nslices)
300: {
301:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

303:   *nslices = ctx->nslices;
304:   return 0;
305: }

307: /*@
308:    EPSEVSLGetSlices - Gets the number of slices in which the interval must be
309:    subdivided.

311:    Not Collective

313:    Input Parameter:
314: .  eps - the eigensolver context

316:    Output Parameter:
317: .  nslices - the number of slices

319:    Level: intermediate

321: .seealso: EPSEVSLSetSlices()
322: @*/
323: PetscErrorCode EPSEVSLGetSlices(EPS eps,PetscInt *nslices)
324: {
327:   PetscUseMethod(eps,"EPSEVSLGetSlices_C",(EPS,PetscInt*),(eps,nslices));
328:   return 0;
329: }

331: static PetscErrorCode EPSEVSLSetRange_EVSL(EPS eps,PetscReal lmin,PetscReal lmax)
332: {
333:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

336:   if (ctx->lmin != lmin || ctx->lmax != lmax) {
337:     ctx->lmin  = lmin;
338:     ctx->lmax  = lmax;
339:     eps->state = EPS_STATE_INITIAL;
340:   }
341:   return 0;
342: }

344: /*@
345:    EPSEVSLSetRange - Defines the numerical range (or field of values) of the problem,
346:    that is, the interval containing all eigenvalues.

348:    Logically Collective on eps

350:    Input Parameters:
351: +  eps  - the eigensolver context
352: .  lmin - left end of the interval
353: -  lmax - right end of the interval

355:    Options Database Key:
356: .  -eps_evsl_range <a,b> - set [a,b] as the numerical range

358:    Notes:
359:    The filter will be most effective if the numerical range is tight, that is, lmin
360:    and lmax are good approximations to the leftmost and rightmost eigenvalues,
361:    respectively. If not set by the user, an approximation is computed internally.

363:    The wanted computational interval specified via EPSSetInterval() must be
364:    contained in the numerical range.

366:    Level: intermediate

368: .seealso: EPSEVSLGetRange(), EPSSetInterval()
369: @*/
370: PetscErrorCode EPSEVSLSetRange(EPS eps,PetscReal lmin,PetscReal lmax)
371: {
375:   PetscTryMethod(eps,"EPSEVSLSetRange_C",(EPS,PetscReal,PetscReal),(eps,lmin,lmax));
376:   return 0;
377: }

379: static PetscErrorCode EPSEVSLGetRange_EVSL(EPS eps,PetscReal *lmin,PetscReal *lmax)
380: {
381:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

383:   if (lmin) *lmin = ctx->lmin;
384:   if (lmax) *lmax = ctx->lmax;
385:   return 0;
386: }

388: /*@
389:    EPSEVSLGetRange - Gets the interval containing all eigenvalues.

391:    Not Collective

393:    Input Parameter:
394: .  eps - the eigensolver context

396:    Output Parameters:
397: +  lmin - left end of the interval
398: -  lmax - right end of the interval

400:    Level: intermediate

402: .seealso: EPSEVSLSetRange()
403: @*/
404: PetscErrorCode EPSEVSLGetRange(EPS eps,PetscReal *lmin,PetscReal *lmax)
405: {
407:   PetscUseMethod(eps,"EPSEVSLGetRange_C",(EPS,PetscReal*,PetscReal*),(eps,lmin,lmax));
408:   return 0;
409: }

411: static PetscErrorCode EPSEVSLSetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
412: {
413:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

415:   ctx->dos = dos;
416:   if (nvec == PETSC_DECIDE || nvec == PETSC_DEFAULT) ctx->nvec = 80;
417:   else {
419:     ctx->nvec = nvec;
420:   }
421:   switch (dos) {
422:     case EPS_EVSL_DOS_KPM:
423:       if (deg == PETSC_DECIDE || deg == PETSC_DEFAULT) ctx->deg = 300;
424:       else {
426:         ctx->deg = deg;
427:       }
428:       break;
429:     case EPS_EVSL_DOS_LANCZOS:
430:       if (steps == PETSC_DECIDE || steps == PETSC_DEFAULT) ctx->steps = 40;
431:       else {
433:         ctx->steps = steps;
434:       }
435:       if (npoints == PETSC_DECIDE || npoints == PETSC_DEFAULT) ctx->npoints = 200;
436:       else {
438:         ctx->npoints = npoints;
439:       }
440:       break;
441:   }
442:   eps->state = EPS_STATE_INITIAL;
443:   return 0;
444: }

446: /*@
447:    EPSEVSLSetDOSParameters - Defines the parameters used for computing the
448:    density of states (DOS) in the EVSL solver.

450:    Logically Collective on eps

452:    Input Parameters:
453: +  eps     - the eigensolver context
454: .  dos     - DOS method, either KPM or Lanczos
455: .  nvec    - number of sample vectors
456: .  deg     - polynomial degree (KPM only)
457: .  steps   - number of Lanczos steps (Lanczos only)
458: -  npoints - number of sample points (Lanczos only)

460:    Options Database Keys:
461: +  -eps_evsl_dos_method <dos> - set the DOS method, either kpm or lanczos
462: .  -eps_evsl_dos_nvec <n> - set the number of sample vectors
463: .  -eps_evsl_dos_degree <n> - set the polynomial degree
464: .  -eps_evsl_dos_steps <n> - set the number of Lanczos steps
465: -  -eps_evsl_dos_npoints <n> - set the number of sample points

467:    Notes:
468:    The density of states (or spectral density) can be approximated with two
469:    methods, kernel polynomial method (KPM) or Lanczos. Some parameters for
470:    these methods can be set by the user with this function, with some of
471:    them being relevant for one of the methods only.

473:    Level: intermediate

475: .seealso: EPSEVSLGetDOSParameters()
476: @*/
477: PetscErrorCode EPSEVSLSetDOSParameters(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
478: {
485:   PetscTryMethod(eps,"EPSEVSLSetDOSParameters_C",(EPS,EPSEVSLDOSMethod,PetscInt,PetscInt,PetscInt,PetscInt),(eps,dos,nvec,deg,steps,npoints));
486:   return 0;
487: }

489: static PetscErrorCode EPSEVSLGetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
490: {
491:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

493:   if (dos)     *dos     = ctx->dos;
494:   if (nvec)    *nvec    = ctx->nvec;
495:   if (deg)     *deg     = ctx->deg;
496:   if (steps)   *steps   = ctx->steps;
497:   if (npoints) *npoints = ctx->npoints;
498:   return 0;
499: }

501: /*@
502:    EPSEVSLGetDOSParameters - Gets the parameters used for computing the
503:    density of states (DOS) in the EVSL solver.

505:    Not Collective

507:    Input Parameter:
508: .  eps - the eigensolver context

510:    Output Parameters:
511: +  dos     - DOS method, either KPM or Lanczos
512: .  nvec    - number of sample vectors
513: .  deg     - polynomial degree (KPM only)
514: .  steps   - number of Lanczos steps (Lanczos only)
515: -  npoints - number of sample points (Lanczos only)

517:    Level: intermediate

519: .seealso: EPSEVSLSetDOSParameters()
520: @*/
521: PetscErrorCode EPSEVSLGetDOSParameters(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
522: {
524:   PetscUseMethod(eps,"EPSEVSLGetDOSParameters_C",(EPS,EPSEVSLDOSMethod*,PetscInt*,PetscInt*,PetscInt*,PetscInt*),(eps,dos,nvec,deg,steps,npoints));
525:   return 0;
526: }

528: static PetscErrorCode EPSEVSLSetPolParameters_EVSL(EPS eps,PetscInt max_deg,PetscReal thresh)
529: {
530:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

532:   if (max_deg == PETSC_DECIDE || max_deg == PETSC_DEFAULT) ctx->max_deg = 10000;
533:   else {
535:     ctx->max_deg = max_deg;
536:   }
537:   if (thresh == PETSC_DECIDE || thresh == PETSC_DEFAULT) ctx->thresh = 0.8;
538:   else {
540:     ctx->thresh = thresh;
541:   }
542:   eps->state = EPS_STATE_INITIAL;
543:   return 0;
544: }

546: /*@
547:    EPSEVSLSetPolParameters - Defines the parameters used for building the
548:    building the polynomial in the EVSL solver.

550:    Logically Collective on eps

552:    Input Parameters:
553: +  eps     - the eigensolver context
554: .  max_deg - maximum degree allowed for the polynomial
555: -  thresh  - threshold for accepting polynomial

557:    Options Database Keys:
558: +  -eps_evsl_pol_max_deg <d> - set maximum polynomial degree
559: -  -eps_evsl_pol_thresh <t> - set the threshold

561:    Level: intermediate

563: .seealso: EPSEVSLGetPolParameters()
564: @*/
565: PetscErrorCode EPSEVSLSetPolParameters(EPS eps,PetscInt max_deg,PetscReal thresh)
566: {
570:   PetscTryMethod(eps,"EPSEVSLSetPolParameters_C",(EPS,PetscInt,PetscReal),(eps,max_deg,thresh));
571:   return 0;
572: }

574: static PetscErrorCode EPSEVSLGetPolParameters_EVSL(EPS eps,PetscInt *max_deg,PetscReal *thresh)
575: {
576:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

578:   if (max_deg) *max_deg = ctx->max_deg;
579:   if (thresh)  *thresh  = ctx->thresh;
580:   return 0;
581: }

583: /*@
584:    EPSEVSLGetPolParameters - Gets the parameters used for building the
585:    polynomial in the EVSL solver.

587:    Not Collective

589:    Input Parameter:
590: .  eps - the eigensolver context

592:    Output Parameters:
593: +  max_deg - the maximum degree of the polynomial
594: -  thresh  - the threshold

596:    Level: intermediate

598: .seealso: EPSEVSLSetPolParameters()
599: @*/
600: PetscErrorCode EPSEVSLGetPolParameters(EPS eps,PetscInt *max_deg,PetscReal *thresh)
601: {
603:   PetscUseMethod(eps,"EPSEVSLGetPolParameters_C",(EPS,PetscInt*,PetscReal*),(eps,max_deg,thresh));
604:   return 0;
605: }

607: static PetscErrorCode EPSEVSLSetDamping_EVSL(EPS eps,EPSEVSLDamping damping)
608: {
609:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

611:   if (ctx->damping != damping) {
612:     ctx->damping = damping;
613:     eps->state   = EPS_STATE_INITIAL;
614:   }
615:   return 0;
616: }

618: /*@
619:    EPSEVSLSetDamping - Set the type of damping to be used in EVSL.

621:    Logically Collective on eps

623:    Input Parameters:
624: +  eps     - the eigensolver context
625: -  damping - the type of damping

627:    Options Database Key:
628: .  -eps_evsl_damping <n> - set the type of damping

630:    Notes:
631:    Damping is applied when building the polynomial to be used when solving the
632:    eigenproblem, and also during estimation of DOS with the KPM method.

634:    Level: intermediate

636: .seealso: EPSEVSLGetDamping(), EPSEVSLSetDOSParameters()
637: @*/
638: PetscErrorCode EPSEVSLSetDamping(EPS eps,EPSEVSLDamping damping)
639: {
642:   PetscTryMethod(eps,"EPSEVSLSetDamping_C",(EPS,EPSEVSLDamping),(eps,damping));
643:   return 0;
644: }

646: static PetscErrorCode EPSEVSLGetDamping_EVSL(EPS eps,EPSEVSLDamping *damping)
647: {
648:   EPS_EVSL *ctx = (EPS_EVSL*)eps->data;

650:   *damping = ctx->damping;
651:   return 0;
652: }

654: /*@
655:    EPSEVSLGetDamping - Gets the type of damping.

657:    Not Collective

659:    Input Parameter:
660: .  eps - the eigensolver context

662:    Output Parameter:
663: .  damping - the type of damping

665:    Level: intermediate

667: .seealso: EPSEVSLSetDamping()
668: @*/
669: PetscErrorCode EPSEVSLGetDamping(EPS eps,EPSEVSLDamping *damping)
670: {
673:   PetscUseMethod(eps,"EPSEVSLGetDamping_C",(EPS,EPSEVSLDamping*),(eps,damping));
674:   return 0;
675: }

677: PetscErrorCode EPSView_EVSL(EPS eps,PetscViewer viewer)
678: {
679:   PetscBool      isascii;
680:   EPS_EVSL       *ctx = (EPS_EVSL*)eps->data;

682:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
683:   if (isascii) {
684:     PetscViewerASCIIPrintf(viewer,"  numerical range = [%g,%g]\n",(double)ctx->lmin,(double)ctx->lmax);
685:     PetscViewerASCIIPrintf(viewer,"  number of slices = %" PetscInt_FMT "\n",ctx->nslices);
686:     PetscViewerASCIIPrintf(viewer,"  type of damping = %s\n",EPSEVSLDampings[ctx->damping]);
687:     PetscViewerASCIIPrintf(viewer,"  computing DOS with %s: nvec=%" PetscInt_FMT ", ",EPSEVSLDOSMethods[ctx->dos],ctx->nvec);
688:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
689:     switch (ctx->dos) {
690:       case EPS_EVSL_DOS_KPM:
691:         PetscViewerASCIIPrintf(viewer,"degree=%" PetscInt_FMT "\n",ctx->deg);
692:         break;
693:       case EPS_EVSL_DOS_LANCZOS:
694:         PetscViewerASCIIPrintf(viewer,"steps=%" PetscInt_FMT ", npoints=%" PetscInt_FMT "\n",ctx->steps,ctx->npoints);
695:         break;
696:     }
697:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
698:     PetscViewerASCIIPrintf(viewer,"  polynomial parameters: max degree = %" PetscInt_FMT ", threshold = %g\n",ctx->max_deg,(double)ctx->thresh);
699:   }
700:   return 0;
701: }

703: PetscErrorCode EPSSetFromOptions_EVSL(EPS eps,PetscOptionItems *PetscOptionsObject)
704: {
705:   PetscReal        array[2]={0,0},th;
706:   PetscInt         k,i1,i2,i3,i4;
707:   PetscBool        flg,flg1;
708:   EPSEVSLDOSMethod dos;
709:   EPSEVSLDamping   damping;
710:   EPS_EVSL         *ctx = (EPS_EVSL*)eps->data;

712:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS EVSL Options");

714:     k = 2;
715:     PetscOptionsRealArray("-eps_evsl_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","EPSEVSLSetRange",array,&k,&flg);
716:     if (flg) {
718:       EPSEVSLSetRange(eps,array[0],array[1]);
719:     }

721:     PetscOptionsInt("-eps_evsl_slices","Number of slices","EPSEVSLSetSlices",ctx->nslices,&i1,&flg);
722:     if (flg) EPSEVSLSetSlices(eps,i1);

724:     PetscOptionsEnum("-eps_evsl_damping","Type of damping","EPSEVSLSetDamping",EPSEVSLDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg);
725:     if (flg) EPSEVSLSetDamping(eps,damping);

727:     EPSEVSLGetDOSParameters(eps,&dos,&i1,&i2,&i3,&i4);
728:     PetscOptionsEnum("-eps_evsl_dos_method","Method to compute the DOS","EPSEVSLSetDOSParameters",EPSEVSLDOSMethods,(PetscEnum)ctx->dos,(PetscEnum*)&dos,&flg);
729:     PetscOptionsInt("-eps_evsl_dos_nvec","Number of sample vectors for DOS","EPSEVSLSetDOSParameters",i1,&i1,&flg1);
730:     if (flg1) flg = PETSC_TRUE;
731:     PetscOptionsInt("-eps_evsl_dos_degree","Polynomial degree used for DOS","EPSEVSLSetDOSParameters",i2,&i2,&flg1);
732:     if (flg1) flg = PETSC_TRUE;
733:     PetscOptionsInt("-eps_evsl_dos_steps","Number of Lanczos steps in DOS","EPSEVSLSetDOSParameters",i3,&i3,&flg1);
734:     if (flg1) flg = PETSC_TRUE;
735:     PetscOptionsInt("-eps_evsl_dos_npoints","Number of sample points used for DOS","EPSEVSLSetDOSParameters",i4,&i4,&flg1);
736:     if (flg || flg1) EPSEVSLSetDOSParameters(eps,dos,i1,i2,i3,i4);

738:     EPSEVSLGetPolParameters(eps,&i1,&th);
739:     PetscOptionsInt("-eps_evsl_pol_max_deg","Maximum degree allowed for the polynomial","EPSEVSLSetPolParameters",i1,&i1,&flg);
740:     PetscOptionsReal("-eps_evsl_pol_threshold","Threshold for accepting polynomial","EPSEVSLSetPolParameters",th,&th,&flg1);
741:     if (flg || flg1) EPSEVSLSetPolParameters(eps,i1,th);

743:   PetscOptionsHeadEnd();
744:   return 0;
745: }

747: PetscErrorCode EPSDestroy_EVSL(EPS eps)
748: {
749:   EPS_EVSL       *ctx = (EPS_EVSL*)eps->data;

751:   if (ctx->initialized) EVSLFinish();
752:   PetscLayoutDestroy(&ctx->map);
753:   PetscFree(ctx->sli);
754:   PetscFree(eps->data);
755:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",NULL);
756:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",NULL);
757:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",NULL);
758:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",NULL);
759:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",NULL);
760:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",NULL);
761:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",NULL);
762:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",NULL);
763:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",NULL);
764:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",NULL);
765:   return 0;
766: }

768: PetscErrorCode EPSReset_EVSL(EPS eps)
769: {
770:   EPS_EVSL       *ctx = (EPS_EVSL*)eps->data;

772:   MatDestroy(&ctx->A);
773:   VecDestroy(&ctx->x);
774:   VecDestroy(&ctx->y);
775:   return 0;
776: }

778: SLEPC_EXTERN PetscErrorCode EPSCreate_EVSL(EPS eps)
779: {
780:   EPS_EVSL       *ctx;

782:   PetscNew(&ctx);
783:   eps->data = (void*)ctx;

785:   ctx->nslices = 0;
786:   ctx->lmin    = PETSC_MIN_REAL;
787:   ctx->lmax    = PETSC_MAX_REAL;
788:   ctx->dos     = EPS_EVSL_DOS_KPM;
789:   ctx->nvec    = 80;
790:   ctx->deg     = 300;
791:   ctx->steps   = 40;
792:   ctx->npoints = 200;
793:   ctx->max_deg = 10000;
794:   ctx->thresh  = 0.8;
795:   ctx->damping = EPS_EVSL_DAMPING_SIGMA;

797:   eps->categ = EPS_CATEGORY_OTHER;

799:   eps->ops->solve          = EPSSolve_EVSL;
800:   eps->ops->setup          = EPSSetUp_EVSL;
801:   eps->ops->setupsort      = EPSSetUpSort_Basic;
802:   eps->ops->setfromoptions = EPSSetFromOptions_EVSL;
803:   eps->ops->destroy        = EPSDestroy_EVSL;
804:   eps->ops->reset          = EPSReset_EVSL;
805:   eps->ops->view           = EPSView_EVSL;
806:   eps->ops->backtransform  = EPSBackTransform_Default;
807:   eps->ops->setdefaultst   = EPSSetDefaultST_NoFactor;

809:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",EPSEVSLSetRange_EVSL);
810:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",EPSEVSLGetRange_EVSL);
811:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",EPSEVSLSetSlices_EVSL);
812:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",EPSEVSLGetSlices_EVSL);
813:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",EPSEVSLSetDOSParameters_EVSL);
814:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",EPSEVSLGetDOSParameters_EVSL);
815:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",EPSEVSLSetPolParameters_EVSL);
816:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",EPSEVSLGetPolParameters_EVSL);
817:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",EPSEVSLSetDamping_EVSL);
818:   PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",EPSEVSLGetDamping_EVSL);
819:   return 0;
820: }