Actual source code: evsl.c
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: int evsl_ierr_ = func(__VA_ARGS__); \
20: PetscStackPop; \
21: PetscCheck(!evsl_ierr_,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling %s: error code %d",PetscStringize(func(__VA_ARGS__)),evsl_ierr_); \
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: PetscFunctionBegin;
52: PetscCallAbort(PetscObjectComm((PetscObject)A),VecPlaceArray(x,(PetscScalar*)xa));
53: PetscCallAbort(PetscObjectComm((PetscObject)A),VecPlaceArray(y,(PetscScalar*)ya));
54: PetscCallAbort(PetscObjectComm((PetscObject)A),MatMult(A,x,y));
55: PetscCallAbort(PetscObjectComm((PetscObject)A),VecResetArray(x));
56: PetscCallAbort(PetscObjectComm((PetscObject)A),VecResetArray(y));
57: PetscFunctionReturnVoid();
58: }
60: static PetscErrorCode EPSSetUp_EVSL(EPS eps)
61: {
62: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
63: PetscMPIInt size,rank;
64: PetscBool isshift;
65: PetscScalar *vinit;
66: PetscReal *mu,ecount,xintv[4],*xdos,*ydos;
67: Vec v0;
68: Mat A;
69: PetscRandom rnd;
71: PetscFunctionBegin;
72: EPSCheckStandard(eps);
73: EPSCheckHermitian(eps);
74: EPSCheckNotStructured(eps);
75: if (eps->nev==0) eps->nev = 1;
76: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
77: PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
79: if (ctx->initialized) EVSLFinish();
80: EVSLStart();
81: ctx->initialized=PETSC_TRUE;
83: /* get number of slices per process */
84: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
85: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
86: if (!ctx->nslices) ctx->nslices = size;
87: PetscCall(PetscLayoutDestroy(&ctx->map));
88: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)eps),PETSC_DECIDE,ctx->nslices,1,&ctx->map));
90: /* get matrix and prepare auxiliary vectors */
91: PetscCall(MatDestroy(&ctx->A));
92: PetscCall(STGetMatrix(eps->st,0,&A));
93: if (size==1) {
94: PetscCall(PetscObjectReference((PetscObject)A));
95: ctx->A = A;
96: } else PetscCall(MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&ctx->A));
97: SetAMatvec(eps->n,&AMatvec_EVSL,(void*)ctx);
98: if (!ctx->x) PetscCall(MatCreateVecsEmpty(ctx->A,&ctx->x,&ctx->y));
99: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
100: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE);
102: if (!eps->which) eps->which=EPS_ALL;
103: PetscCheck(eps->which==EPS_ALL && eps->inta!=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires setting an interval with EPSSetInterval()");
105: /* estimate numerical range */
106: if (ctx->estimrange || ctx->lmin == PETSC_MIN_REAL || ctx->lmax == PETSC_MAX_REAL) {
107: PetscCall(MatCreateVecs(ctx->A,&v0,NULL));
108: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
109: PetscCall(BVGetRandomContext(eps->V,&rnd));
110: PetscCall(VecSetRandom(v0,rnd));
111: PetscCall(VecGetArray(v0,&vinit));
112: PetscCallEVSL(LanTrbounds,50,200,eps->tol,vinit,1,&ctx->lmin,&ctx->lmax,NULL);
113: PetscCall(VecRestoreArray(v0,&vinit));
114: PetscCall(VecDestroy(&v0));
115: ctx->estimrange = PETSC_TRUE; /* estimate if called again with another matrix */
116: }
117: PetscCheck(ctx->lmin<=eps->inta && ctx->lmax>=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The requested interval [%g,%g] must be contained in the numerical range [%g,%g]",(double)eps->inta,(double)eps->intb,(double)ctx->lmin,(double)ctx->lmax);
118: xintv[0] = eps->inta;
119: xintv[1] = eps->intb;
120: xintv[2] = ctx->lmin;
121: xintv[3] = ctx->lmax;
123: /* estimate number of eigenvalues in the interval */
124: switch (ctx->dos) {
125: case EPS_EVSL_DOS_KPM:
126: PetscCall(PetscMalloc1(ctx->deg+1,&mu));
127: if (!rank) PetscCallEVSL(kpmdos,ctx->deg,(int)ctx->damping,ctx->nvec,xintv,mu,&ecount);
128: PetscCallMPI(MPI_Bcast(mu,ctx->deg+1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
129: break;
130: case EPS_EVSL_DOS_LANCZOS:
131: PetscCall(PetscMalloc2(ctx->npoints,&xdos,ctx->npoints,&ydos));
132: if (!rank) PetscCallEVSL(LanDos,ctx->nvec,PetscMin(ctx->steps,eps->n/2),ctx->npoints,xdos,ydos,&ecount,xintv);
133: PetscCallMPI(MPI_Bcast(xdos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
134: PetscCallMPI(MPI_Bcast(ydos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
135: break;
136: default:
137: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid DOS method");
138: }
139: PetscCallMPI(MPI_Bcast(&ecount,1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
141: PetscCall(PetscInfo(eps,"Estimated eigenvalue count in the interval: %g\n",ecount));
142: eps->ncv = (PetscInt)PetscCeilReal(1.5*ecount);
144: /* slice the spectrum */
145: PetscCall(PetscFree(ctx->sli));
146: PetscCall(PetscMalloc1(ctx->nslices+1,&ctx->sli));
147: if (ctx->dos == EPS_EVSL_DOS_KPM) {
148: PetscCallEVSL(spslicer,ctx->sli,mu,ctx->deg,xintv,ctx->nslices,10*(PetscInt)ecount);
149: PetscCall(PetscFree(mu));
150: } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
151: spslicer2(xdos,ydos,ctx->nslices,ctx->npoints,ctx->sli);
152: PetscCall(PetscFree2(xdos,ydos));
153: }
155: /* approximate number of eigenvalues wanted in each slice */
156: ctx->nev = (PetscInt)(1.0 + ecount/(PetscReal)ctx->nslices) + 2;
158: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
159: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
160: PetscCall(EPSAllocateSolution(eps,0));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode EPSSolve_EVSL(EPS eps)
165: {
166: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
167: PetscInt i,j,k=0,sl,mlan,nevout,*ind,nevmax,rstart,rend,*nevloc,*disp,N;
168: PetscReal *res,xintv[4],*errest;
169: PetscScalar *lam,*X,*Y,*vinit,*eigr;
170: PetscMPIInt size,rank;
171: PetscRandom rnd;
172: Vec v,w,v0,x;
173: VecScatter vs;
174: IS is;
175: polparams pol;
177: PetscFunctionBegin;
178: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
179: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
180: PetscCall(PetscLayoutGetRange(ctx->map,&rstart,&rend));
181: nevmax = (rend-rstart)*ctx->nev;
182: PetscCall(MatCreateVecs(ctx->A,&v0,NULL));
183: PetscCall(BVGetRandomContext(eps->V,&rnd));
184: PetscCall(VecSetRandom(v0,rnd));
185: PetscCall(VecGetArray(v0,&vinit));
186: PetscCall(PetscMalloc5(size,&nevloc,size+1,&disp,nevmax,&eigr,nevmax,&errest,nevmax*eps->n,&X));
187: mlan = PetscMin(PetscMax(5*ctx->nev,300),eps->n);
188: for (sl=rstart; sl<rend; sl++) {
189: xintv[0] = ctx->sli[sl];
190: xintv[1] = ctx->sli[sl+1];
191: xintv[2] = ctx->lmin;
192: xintv[3] = ctx->lmax;
193: PetscCall(PetscInfo(ctx->A,"Subinterval %" PetscInt_FMT ": [%.4e, %.4e]\n",sl+1,xintv[0],xintv[1]));
194: set_pol_def(&pol);
195: pol.max_deg = ctx->max_deg;
196: pol.damping = (int)ctx->damping;
197: pol.thresh_int = ctx->thresh;
198: find_pol(xintv,&pol);
199: PetscCall(PetscInfo(ctx->A,"Polynomial [type = %" PetscInt_FMT "], deg %" PetscInt_FMT ", bar %e gam %e\n",pol.type,pol.deg,pol.bar,pol.gam));
200: PetscCallEVSL(ChebLanNr,xintv,mlan,eps->tol,vinit,&pol,&nevout,&lam,&Y,&res,NULL);
201: PetscCheck(k+nevout<=nevmax,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
202: free_pol(&pol);
203: PetscCall(PetscInfo(ctx->A,"Computed %" PetscInt_FMT " eigenvalues\n",nevout));
204: PetscCall(PetscMalloc1(nevout,&ind));
205: sort_double(nevout,lam,ind);
206: for (i=0;i<nevout;i++) {
207: eigr[i+k] = lam[i];
208: errest[i+k] = res[ind[i]];
209: PetscCall(PetscArraycpy(X+(i+k)*eps->n,Y+ind[i]*eps->n,eps->n));
210: }
211: k += nevout;
212: if (lam) evsl_Free(lam);
213: if (Y) evsl_Free_device(Y);
214: if (res) evsl_Free(res);
215: PetscCall(PetscFree(ind));
216: }
217: PetscCall(VecRestoreArray(v0,&vinit));
218: PetscCall(VecDestroy(&v0));
220: /* gather eigenvalues computed by each MPI process */
221: PetscCallMPI(MPI_Allgather(&k,1,MPIU_INT,nevloc,1,MPIU_INT,PetscObjectComm((PetscObject)eps)));
222: eps->nev = nevloc[0];
223: disp[0] = 0;
224: for (i=1;i<size;i++) {
225: eps->nev += nevloc[i];
226: disp[i] = disp[i-1]+nevloc[i-1];
227: }
228: disp[size] = disp[size-1]+nevloc[size-1];
229: PetscCheck(eps->nev<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
230: PetscCallMPI(MPI_Allgatherv(eigr,k,MPIU_SCALAR,eps->eigr,nevloc,disp,MPIU_SCALAR,PetscObjectComm((PetscObject)eps)));
231: PetscCallMPI(MPI_Allgatherv(errest,k,MPIU_REAL,eps->errest,nevloc,disp,MPIU_REAL,PetscObjectComm((PetscObject)eps)));
232: eps->nconv = eps->nev;
233: eps->its = 1;
234: eps->reason = EPS_CONVERGED_TOL;
236: /* scatter computed eigenvectors and store them in eps->V */
237: PetscCall(BVCreateVec(eps->V,&w));
238: for (i=0;i<size;i++) {
239: N = (rank==i)? eps->n: 0;
240: PetscCall(VecCreateSeq(PETSC_COMM_SELF,N,&x));
241: PetscCall(VecSetFromOptions(x));
242: PetscCall(ISCreateStride(PETSC_COMM_SELF,N,0,1,&is));
243: PetscCall(VecScatterCreate(x,is,w,is,&vs));
244: PetscCall(ISDestroy(&is));
245: for (j=disp[i];j<disp[i+1];j++) {
246: PetscCall(BVGetColumn(eps->V,j,&v));
247: if (rank==i) PetscCall(VecPlaceArray(x,X+(j-disp[i])*eps->n));
248: PetscCall(VecScatterBegin(vs,x,v,INSERT_VALUES,SCATTER_FORWARD));
249: PetscCall(VecScatterEnd(vs,x,v,INSERT_VALUES,SCATTER_FORWARD));
250: if (rank==i) PetscCall(VecResetArray(x));
251: PetscCall(BVRestoreColumn(eps->V,j,&v));
252: }
253: PetscCall(VecScatterDestroy(&vs));
254: PetscCall(VecDestroy(&x));
255: }
256: PetscCall(VecDestroy(&w));
257: PetscCall(PetscFree5(nevloc,disp,eigr,errest,X));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode EPSEVSLSetSlices_EVSL(EPS eps,PetscInt nslices)
262: {
263: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
265: PetscFunctionBegin;
266: if (nslices == PETSC_DECIDE || nslices == PETSC_DEFAULT) nslices = 0;
267: else PetscCheck(nslices>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Number of slices must be 1 at least");
268: if (ctx->nslices != nslices) {
269: ctx->nslices = nslices;
270: eps->state = EPS_STATE_INITIAL;
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: EPSEVSLSetSlices - Set the number of slices in which the interval must be
277: subdivided.
279: Logically Collective
281: Input Parameters:
282: + eps - the linear eigensolver context
283: - nslices - the number of slices
285: Options Database Key:
286: . -eps_evsl_slices \<nslices\> - set the number of slices
288: Notes:
289: By default, one slice per MPI process is used. Depending on the number of
290: eigenvalues, using more slices may be beneficial, but very narrow subintervals
291: imply higher polynomial degree.
293: See the documentation of EVSL {cite:p}`Li19` for details.
295: Level: advanced
297: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLGetSlices()`
298: @*/
299: PetscErrorCode EPSEVSLSetSlices(EPS eps,PetscInt nslices)
300: {
301: PetscFunctionBegin;
304: PetscTryMethod(eps,"EPSEVSLSetSlices_C",(EPS,PetscInt),(eps,nslices));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode EPSEVSLGetSlices_EVSL(EPS eps,PetscInt *nslices)
309: {
310: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
312: PetscFunctionBegin;
313: *nslices = ctx->nslices;
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@
318: EPSEVSLGetSlices - Gets the number of slices in which the interval must be
319: subdivided.
321: Not Collective
323: Input Parameter:
324: . eps - the linear eigensolver context
326: Output Parameter:
327: . nslices - the number of slices
329: Level: advanced
331: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLSetSlices()`
332: @*/
333: PetscErrorCode EPSEVSLGetSlices(EPS eps,PetscInt *nslices)
334: {
335: PetscFunctionBegin;
337: PetscAssertPointer(nslices,2);
338: PetscUseMethod(eps,"EPSEVSLGetSlices_C",(EPS,PetscInt*),(eps,nslices));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode EPSEVSLSetRange_EVSL(EPS eps,PetscReal lmin,PetscReal lmax)
343: {
344: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
346: PetscFunctionBegin;
347: PetscCheck(lmin<lmax,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be lmin<lmax");
348: if (ctx->lmin != lmin || ctx->lmax != lmax) {
349: ctx->lmin = lmin;
350: ctx->lmax = lmax;
351: eps->state = EPS_STATE_INITIAL;
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: EPSEVSLSetRange - Defines the numerical range (or field of values) of the problem,
358: that is, the interval containing all eigenvalues.
360: Logically Collective
362: Input Parameters:
363: + eps - the linear eigensolver context
364: . lmin - left end of the interval
365: - lmax - right end of the interval
367: Options Database Key:
368: . -eps_evsl_range <lmin,lmax> - set $[\lambda_\mathrm{min},\lambda_\mathrm{max}]$ as the numerical range
370: Notes:
371: The filter will be most effective if the numerical range is tight, that is, `lmin`
372: and `lmax` are good approximations to the leftmost and rightmost eigenvalues,
373: respectively. If not set by the user, an approximation is computed internally.
375: The wanted computational interval specified via `EPSSetInterval()` must be
376: contained in the numerical range.
378: See the documentation of EVSL {cite:p}`Li19` for details.
380: Level: advanced
382: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLGetRange()`, `EPSSetInterval()`
383: @*/
384: PetscErrorCode EPSEVSLSetRange(EPS eps,PetscReal lmin,PetscReal lmax)
385: {
386: PetscFunctionBegin;
390: PetscTryMethod(eps,"EPSEVSLSetRange_C",(EPS,PetscReal,PetscReal),(eps,lmin,lmax));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode EPSEVSLGetRange_EVSL(EPS eps,PetscReal *lmin,PetscReal *lmax)
395: {
396: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
398: PetscFunctionBegin;
399: if (lmin) *lmin = ctx->lmin;
400: if (lmax) *lmax = ctx->lmax;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: EPSEVSLGetRange - Gets the interval containing all eigenvalues.
407: Not Collective
409: Input Parameter:
410: . eps - the linear eigensolver context
412: Output Parameters:
413: + lmin - left end of the interval
414: - lmax - right end of the interval
416: Level: advanced
418: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLSetRange()`
419: @*/
420: PetscErrorCode EPSEVSLGetRange(EPS eps,PetscReal *lmin,PetscReal *lmax)
421: {
422: PetscFunctionBegin;
424: PetscUseMethod(eps,"EPSEVSLGetRange_C",(EPS,PetscReal*,PetscReal*),(eps,lmin,lmax));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode EPSEVSLSetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
429: {
430: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
432: PetscFunctionBegin;
433: ctx->dos = dos;
434: if (nvec == PETSC_DETERMINE) ctx->nvec = 80;
435: else if (nvec != PETSC_CURRENT) {
436: PetscCheck(nvec>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The nvec argument must be > 0");
437: ctx->nvec = nvec;
438: }
439: switch (dos) {
440: case EPS_EVSL_DOS_KPM:
441: if (deg == PETSC_DETERMINE) ctx->deg = 300;
442: else if (deg != PETSC_CURRENT) {
443: PetscCheck(deg>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The deg argument must be > 0");
444: ctx->deg = deg;
445: }
446: break;
447: case EPS_EVSL_DOS_LANCZOS:
448: if (steps == PETSC_DETERMINE) ctx->steps = 40;
449: else if (steps != PETSC_CURRENT) {
450: PetscCheck(steps>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The steps argument must be > 0");
451: ctx->steps = steps;
452: }
453: if (npoints == PETSC_DETERMINE) ctx->npoints = 200;
454: else if (npoints != PETSC_CURRENT) {
455: PetscCheck(npoints>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npoints argument must be > 0");
456: ctx->npoints = npoints;
457: }
458: break;
459: }
460: eps->state = EPS_STATE_INITIAL;
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: EPSEVSLSetDOSParameters - Defines the parameters used for computing the
466: density of states (DOS) in the EVSL solver.
468: Logically Collective
470: Input Parameters:
471: + eps - the linear eigensolver context
472: . dos - DOS method, see `EPSEVSLDOSMethod` for possible values
473: . nvec - number of sample vectors
474: . deg - polynomial degree (KPM only)
475: . steps - number of Lanczos steps (Lanczos only)
476: - npoints - number of sample points (Lanczos only)
478: Options Database Keys:
479: + -eps_evsl_dos_method \<dos\> - set the DOS method, either `kpm` or `lanczos`
480: . -eps_evsl_dos_nvec \<nvec\> - set the number of sample vectors
481: . -eps_evsl_dos_degree \<deg\> - set the polynomial degree
482: . -eps_evsl_dos_steps \<steps\> - set the number of Lanczos steps
483: - -eps_evsl_dos_npoints \<npoints\> - set the number of sample points
485: Notes:
486: The density of states (or spectral density) can be approximated with two
487: methods, Kernel Polynomial Method (KPM) or Lanczos. Some parameters for
488: these methods can be set by the user with this function, with some of
489: them being relevant for one of the methods only.
491: For the integer argumens, you can use `PETSC_CURRENT` to keep the current
492: value, and `PETSC_DETERMINE` to set them to a reasonable default.
494: See the documentation of EVSL {cite:p}`Li19` for details.
496: Level: advanced
498: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLGetDOSParameters()`, `EPSEVSLDOSMethod`
499: @*/
500: PetscErrorCode EPSEVSLSetDOSParameters(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
501: {
502: PetscFunctionBegin;
509: PetscTryMethod(eps,"EPSEVSLSetDOSParameters_C",(EPS,EPSEVSLDOSMethod,PetscInt,PetscInt,PetscInt,PetscInt),(eps,dos,nvec,deg,steps,npoints));
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: static PetscErrorCode EPSEVSLGetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
514: {
515: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
517: PetscFunctionBegin;
518: if (dos) *dos = ctx->dos;
519: if (nvec) *nvec = ctx->nvec;
520: if (deg) *deg = ctx->deg;
521: if (steps) *steps = ctx->steps;
522: if (npoints) *npoints = ctx->npoints;
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@
527: EPSEVSLGetDOSParameters - Gets the parameters used for computing the
528: density of states (DOS) in the EVSL solver.
530: Not Collective
532: Input Parameter:
533: . eps - the linear eigensolver context
535: Output Parameters:
536: + dos - DOS method
537: . nvec - number of sample vectors
538: . deg - polynomial degree (KPM only)
539: . steps - number of Lanczos steps (Lanczos only)
540: - npoints - number of sample points (Lanczos only)
542: Level: advanced
544: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLSetDOSParameters()`
545: @*/
546: PetscErrorCode EPSEVSLGetDOSParameters(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
547: {
548: PetscFunctionBegin;
550: PetscUseMethod(eps,"EPSEVSLGetDOSParameters_C",(EPS,EPSEVSLDOSMethod*,PetscInt*,PetscInt*,PetscInt*,PetscInt*),(eps,dos,nvec,deg,steps,npoints));
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode EPSEVSLSetPolParameters_EVSL(EPS eps,PetscInt max_deg,PetscReal thresh)
555: {
556: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
558: PetscFunctionBegin;
559: if (max_deg == PETSC_DETERMINE) ctx->max_deg = 10000;
560: else if (max_deg != PETSC_CURRENT) {
561: PetscCheck(max_deg>2,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The max_deg argument must be > 2");
562: ctx->max_deg = max_deg;
563: }
564: if (thresh == (PetscReal)PETSC_DETERMINE) ctx->thresh = 0.8;
565: else if (thresh != (PetscReal)PETSC_CURRENT) {
566: PetscCheck(thresh>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The thresh argument must be > 0.0");
567: ctx->thresh = thresh;
568: }
569: eps->state = EPS_STATE_INITIAL;
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /*@
574: EPSEVSLSetPolParameters - Defines the parameters used for building the
575: building the polynomial in the EVSL solver.
577: Logically Collective
579: Input Parameters:
580: + eps - the linear eigensolver context
581: . max_deg - maximum degree allowed for the polynomial
582: - thresh - threshold for accepting a value
584: Options Database Keys:
585: + -eps_evsl_pol_max_deg \<max_deg\> - set maximum polynomial degree
586: - -eps_evsl_pol_thresh \<thresh\> - set the threshold
588: Notes:
589: `PETSC_CURRENT` can be used to preserve the current value of any of the
590: arguments, and `PETSC_DETERMINE` to set them to a default value.
592: See the documentation of EVSL {cite:p}`Li19` for details.
594: Level: advanced
596: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLGetPolParameters()`
597: @*/
598: PetscErrorCode EPSEVSLSetPolParameters(EPS eps,PetscInt max_deg,PetscReal thresh)
599: {
600: PetscFunctionBegin;
604: PetscTryMethod(eps,"EPSEVSLSetPolParameters_C",(EPS,PetscInt,PetscReal),(eps,max_deg,thresh));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: static PetscErrorCode EPSEVSLGetPolParameters_EVSL(EPS eps,PetscInt *max_deg,PetscReal *thresh)
609: {
610: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
612: PetscFunctionBegin;
613: if (max_deg) *max_deg = ctx->max_deg;
614: if (thresh) *thresh = ctx->thresh;
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: /*@
619: EPSEVSLGetPolParameters - Gets the parameters used for building the
620: polynomial in the EVSL solver.
622: Not Collective
624: Input Parameter:
625: . eps - the linear eigensolver context
627: Output Parameters:
628: + max_deg - the maximum degree of the polynomial
629: - thresh - the threshold
631: Level: advanced
633: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLSetPolParameters()`
634: @*/
635: PetscErrorCode EPSEVSLGetPolParameters(EPS eps,PetscInt *max_deg,PetscReal *thresh)
636: {
637: PetscFunctionBegin;
639: PetscUseMethod(eps,"EPSEVSLGetPolParameters_C",(EPS,PetscInt*,PetscReal*),(eps,max_deg,thresh));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode EPSEVSLSetDamping_EVSL(EPS eps,EPSEVSLDamping damping)
644: {
645: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
647: PetscFunctionBegin;
648: if (ctx->damping != damping) {
649: ctx->damping = damping;
650: eps->state = EPS_STATE_INITIAL;
651: }
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: EPSEVSLSetDamping - Set the type of damping to be used in EVSL.
658: Logically Collective
660: Input Parameters:
661: + eps - the linear eigensolver context
662: - damping - the type of damping, see `EPSEVSLDamping` for possible values
664: Options Database Key:
665: . -eps_evsl_damping \<damping\> - set the type of damping
667: Notes:
668: Damping is applied when building the polynomial to be used when solving the
669: eigenproblem, and also during estimation of the DOS with the KPM method.
671: See the documentation of EVSL {cite:p}`Li19` for details.
673: Level: advanced
675: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLGetDamping()`, `EPSEVSLSetDOSParameters()`
676: @*/
677: PetscErrorCode EPSEVSLSetDamping(EPS eps,EPSEVSLDamping damping)
678: {
679: PetscFunctionBegin;
682: PetscTryMethod(eps,"EPSEVSLSetDamping_C",(EPS,EPSEVSLDamping),(eps,damping));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: static PetscErrorCode EPSEVSLGetDamping_EVSL(EPS eps,EPSEVSLDamping *damping)
687: {
688: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
690: PetscFunctionBegin;
691: *damping = ctx->damping;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: EPSEVSLGetDamping - Gets the type of damping.
698: Not Collective
700: Input Parameter:
701: . eps - the linear eigensolver context
703: Output Parameter:
704: . damping - the type of damping
706: Level: advanced
708: .seealso: [](ch:eps), `EPSEVSL`, `EPSEVSLSetDamping()`
709: @*/
710: PetscErrorCode EPSEVSLGetDamping(EPS eps,EPSEVSLDamping *damping)
711: {
712: PetscFunctionBegin;
714: PetscAssertPointer(damping,2);
715: PetscUseMethod(eps,"EPSEVSLGetDamping_C",(EPS,EPSEVSLDamping*),(eps,damping));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: static PetscErrorCode EPSView_EVSL(EPS eps,PetscViewer viewer)
720: {
721: PetscBool isascii;
722: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
724: PetscFunctionBegin;
725: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
726: if (isascii) {
727: PetscCall(PetscViewerASCIIPrintf(viewer," numerical range = [%g,%g]\n",(double)ctx->lmin,(double)ctx->lmax));
728: PetscCall(PetscViewerASCIIPrintf(viewer," number of slices = %" PetscInt_FMT "\n",ctx->nslices));
729: PetscCall(PetscViewerASCIIPrintf(viewer," type of damping = %s\n",EPSEVSLDampings[ctx->damping]));
730: PetscCall(PetscViewerASCIIPrintf(viewer," computing DOS with %s: nvec=%" PetscInt_FMT ", ",EPSEVSLDOSMethods[ctx->dos],ctx->nvec));
731: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
732: switch (ctx->dos) {
733: case EPS_EVSL_DOS_KPM:
734: PetscCall(PetscViewerASCIIPrintf(viewer,"degree=%" PetscInt_FMT "\n",ctx->deg));
735: break;
736: case EPS_EVSL_DOS_LANCZOS:
737: PetscCall(PetscViewerASCIIPrintf(viewer,"steps=%" PetscInt_FMT ", npoints=%" PetscInt_FMT "\n",ctx->steps,ctx->npoints));
738: break;
739: }
740: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
741: PetscCall(PetscViewerASCIIPrintf(viewer," polynomial parameters: max degree = %" PetscInt_FMT ", threshold = %g\n",ctx->max_deg,(double)ctx->thresh));
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode EPSSetFromOptions_EVSL(EPS eps,PetscOptionItems PetscOptionsObject)
747: {
748: PetscReal array[2]={0,0},th;
749: PetscInt k,i1,i2,i3,i4;
750: PetscBool flg,flg1;
751: EPSEVSLDOSMethod dos;
752: EPSEVSLDamping damping;
753: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
755: PetscFunctionBegin;
756: PetscOptionsHeadBegin(PetscOptionsObject,"EPS EVSL Options");
758: k = 2;
759: PetscCall(PetscOptionsRealArray("-eps_evsl_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","EPSEVSLSetRange",array,&k,&flg));
760: if (flg) {
761: PetscCheck(k>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_SIZ,"Must pass two values in -eps_evsl_range (comma-separated without spaces)");
762: PetscCall(EPSEVSLSetRange(eps,array[0],array[1]));
763: }
765: PetscCall(PetscOptionsInt("-eps_evsl_slices","Number of slices","EPSEVSLSetSlices",ctx->nslices,&i1,&flg));
766: if (flg) PetscCall(EPSEVSLSetSlices(eps,i1));
768: PetscCall(PetscOptionsEnum("-eps_evsl_damping","Type of damping","EPSEVSLSetDamping",EPSEVSLDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg));
769: if (flg) PetscCall(EPSEVSLSetDamping(eps,damping));
771: PetscCall(EPSEVSLGetDOSParameters(eps,&dos,&i1,&i2,&i3,&i4));
772: PetscCall(PetscOptionsEnum("-eps_evsl_dos_method","Method to compute the DOS","EPSEVSLSetDOSParameters",EPSEVSLDOSMethods,(PetscEnum)ctx->dos,(PetscEnum*)&dos,&flg));
773: PetscCall(PetscOptionsInt("-eps_evsl_dos_nvec","Number of sample vectors for DOS","EPSEVSLSetDOSParameters",i1,&i1,&flg1));
774: if (flg1) flg = PETSC_TRUE;
775: PetscCall(PetscOptionsInt("-eps_evsl_dos_degree","Polynomial degree used for DOS","EPSEVSLSetDOSParameters",i2,&i2,&flg1));
776: if (flg1) flg = PETSC_TRUE;
777: PetscCall(PetscOptionsInt("-eps_evsl_dos_steps","Number of Lanczos steps in DOS","EPSEVSLSetDOSParameters",i3,&i3,&flg1));
778: if (flg1) flg = PETSC_TRUE;
779: PetscCall(PetscOptionsInt("-eps_evsl_dos_npoints","Number of sample points used for DOS","EPSEVSLSetDOSParameters",i4,&i4,&flg1));
780: if (flg || flg1) PetscCall(EPSEVSLSetDOSParameters(eps,dos,i1,i2,i3,i4));
782: PetscCall(EPSEVSLGetPolParameters(eps,&i1,&th));
783: PetscCall(PetscOptionsInt("-eps_evsl_pol_max_deg","Maximum degree allowed for the polynomial","EPSEVSLSetPolParameters",i1,&i1,&flg));
784: PetscCall(PetscOptionsReal("-eps_evsl_pol_threshold","Threshold for accepting polynomial","EPSEVSLSetPolParameters",th,&th,&flg1));
785: if (flg || flg1) PetscCall(EPSEVSLSetPolParameters(eps,i1,th));
787: PetscOptionsHeadEnd();
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: static PetscErrorCode EPSDestroy_EVSL(EPS eps)
792: {
793: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
795: PetscFunctionBegin;
796: if (ctx->initialized) EVSLFinish();
797: PetscCall(PetscLayoutDestroy(&ctx->map));
798: PetscCall(PetscFree(ctx->sli));
799: PetscCall(PetscFree(eps->data));
800: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",NULL));
801: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",NULL));
802: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",NULL));
803: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",NULL));
804: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",NULL));
805: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",NULL));
806: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",NULL));
807: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",NULL));
808: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",NULL));
809: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",NULL));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: static PetscErrorCode EPSReset_EVSL(EPS eps)
814: {
815: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
817: PetscFunctionBegin;
818: PetscCall(MatDestroy(&ctx->A));
819: PetscCall(VecDestroy(&ctx->x));
820: PetscCall(VecDestroy(&ctx->y));
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
824: /*MC
825: EPSEVSL - EPSEVSL = "evsl" - A wrapper to EVSL {cite:p}`Li19`.
827: Note:
828: The Eigenvalues Slicing Library {EVSL} is intended for the case that
829: all eigenvalues inside a given interval are requested, for Hermitian
830: problems. It is based on polynomial filters, and has similarities with
831: the algorithms implemented in `STFILTER`.
833: Level: beginner
835: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
836: M*/
837: SLEPC_EXTERN PetscErrorCode EPSCreate_EVSL(EPS eps)
838: {
839: EPS_EVSL *ctx;
841: PetscFunctionBegin;
842: PetscCall(PetscNew(&ctx));
843: eps->data = (void*)ctx;
845: ctx->nslices = 0;
846: ctx->lmin = PETSC_MIN_REAL;
847: ctx->lmax = PETSC_MAX_REAL;
848: ctx->dos = EPS_EVSL_DOS_KPM;
849: ctx->nvec = 80;
850: ctx->deg = 300;
851: ctx->steps = 40;
852: ctx->npoints = 200;
853: ctx->max_deg = 10000;
854: ctx->thresh = 0.8;
855: ctx->damping = EPS_EVSL_DAMPING_SIGMA;
857: eps->categ = EPS_CATEGORY_OTHER;
859: eps->ops->solve = EPSSolve_EVSL;
860: eps->ops->setup = EPSSetUp_EVSL;
861: eps->ops->setupsort = EPSSetUpSort_Basic;
862: eps->ops->setfromoptions = EPSSetFromOptions_EVSL;
863: eps->ops->destroy = EPSDestroy_EVSL;
864: eps->ops->reset = EPSReset_EVSL;
865: eps->ops->view = EPSView_EVSL;
866: eps->ops->backtransform = EPSBackTransform_Default;
867: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
869: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",EPSEVSLSetRange_EVSL));
870: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",EPSEVSLGetRange_EVSL));
871: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",EPSEVSLSetSlices_EVSL));
872: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",EPSEVSLGetSlices_EVSL));
873: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",EPSEVSLSetDOSParameters_EVSL));
874: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",EPSEVSLGetDOSParameters_EVSL));
875: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",EPSEVSLSetPolParameters_EVSL));
876: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",EPSEVSLGetPolParameters_EVSL));
877: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",EPSEVSLSetDamping_EVSL));
878: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",EPSEVSLGetDamping_EVSL));
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }