Actual source code: evsl.c
slepc-main 2024-11-09
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: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
76: PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
78: if (ctx->initialized) EVSLFinish();
79: EVSLStart();
80: ctx->initialized=PETSC_TRUE;
82: /* get number of slices per process */
83: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
84: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
85: if (!ctx->nslices) ctx->nslices = size;
86: PetscCall(PetscLayoutDestroy(&ctx->map));
87: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)eps),PETSC_DECIDE,ctx->nslices,1,&ctx->map));
89: /* get matrix and prepare auxiliary vectors */
90: PetscCall(MatDestroy(&ctx->A));
91: PetscCall(STGetMatrix(eps->st,0,&A));
92: if (size==1) {
93: PetscCall(PetscObjectReference((PetscObject)A));
94: ctx->A = A;
95: } else PetscCall(MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&ctx->A));
96: SetAMatvec(eps->n,&AMatvec_EVSL,(void*)ctx);
97: if (!ctx->x) PetscCall(MatCreateVecsEmpty(ctx->A,&ctx->x,&ctx->y));
98: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
99: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE);
101: if (!eps->which) eps->which=EPS_ALL;
102: PetscCheck(eps->which==EPS_ALL && eps->inta!=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires setting an interval with EPSSetInterval()");
104: /* estimate numerical range */
105: if (ctx->estimrange || ctx->lmin == PETSC_MIN_REAL || ctx->lmax == PETSC_MAX_REAL) {
106: PetscCall(MatCreateVecs(ctx->A,&v0,NULL));
107: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
108: PetscCall(BVGetRandomContext(eps->V,&rnd));
109: PetscCall(VecSetRandom(v0,rnd));
110: PetscCall(VecGetArray(v0,&vinit));
111: PetscCallEVSL(LanTrbounds,50,200,eps->tol,vinit,1,&ctx->lmin,&ctx->lmax,NULL);
112: PetscCall(VecRestoreArray(v0,&vinit));
113: PetscCall(VecDestroy(&v0));
114: ctx->estimrange = PETSC_TRUE; /* estimate if called again with another matrix */
115: }
116: 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);
117: xintv[0] = eps->inta;
118: xintv[1] = eps->intb;
119: xintv[2] = ctx->lmin;
120: xintv[3] = ctx->lmax;
122: /* estimate number of eigenvalues in the interval */
123: switch (ctx->dos) {
124: case EPS_EVSL_DOS_KPM:
125: PetscCall(PetscMalloc1(ctx->deg+1,&mu));
126: if (!rank) PetscCallEVSL(kpmdos,ctx->deg,(int)ctx->damping,ctx->nvec,xintv,mu,&ecount);
127: PetscCallMPI(MPI_Bcast(mu,ctx->deg+1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
128: break;
129: case EPS_EVSL_DOS_LANCZOS:
130: PetscCall(PetscMalloc2(ctx->npoints,&xdos,ctx->npoints,&ydos));
131: if (!rank) PetscCallEVSL(LanDos,ctx->nvec,PetscMin(ctx->steps,eps->n/2),ctx->npoints,xdos,ydos,&ecount,xintv);
132: PetscCallMPI(MPI_Bcast(xdos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
133: PetscCallMPI(MPI_Bcast(ydos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
134: break;
135: default:
136: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid DOS method");
137: }
138: PetscCallMPI(MPI_Bcast(&ecount,1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps)));
140: PetscCall(PetscInfo(eps,"Estimated eigenvalue count in the interval: %g\n",ecount));
141: eps->ncv = (PetscInt)PetscCeilReal(1.5*ecount);
143: /* slice the spectrum */
144: PetscCall(PetscFree(ctx->sli));
145: PetscCall(PetscMalloc1(ctx->nslices+1,&ctx->sli));
146: if (ctx->dos == EPS_EVSL_DOS_KPM) {
147: PetscCallEVSL(spslicer,ctx->sli,mu,ctx->deg,xintv,ctx->nslices,10*(PetscInt)ecount);
148: PetscCall(PetscFree(mu));
149: } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
150: spslicer2(xdos,ydos,ctx->nslices,ctx->npoints,ctx->sli);
151: PetscCall(PetscFree2(xdos,ydos));
152: }
154: /* approximate number of eigenvalues wanted in each slice */
155: ctx->nev = (PetscInt)(1.0 + ecount/(PetscReal)ctx->nslices) + 2;
157: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
158: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
159: PetscCall(EPSAllocateSolution(eps,0));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode EPSSolve_EVSL(EPS eps)
164: {
165: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
166: PetscInt i,j,k=0,sl,mlan,nevout,*ind,nevmax,rstart,rend,*nevloc,*disp,N;
167: PetscReal *res,xintv[4],*errest;
168: PetscScalar *lam,*X,*Y,*vinit,*eigr;
169: PetscMPIInt size,rank;
170: PetscRandom rnd;
171: Vec v,w,v0,x;
172: VecScatter vs;
173: IS is;
174: polparams pol;
176: PetscFunctionBegin;
177: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
178: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
179: PetscCall(PetscLayoutGetRange(ctx->map,&rstart,&rend));
180: nevmax = (rend-rstart)*ctx->nev;
181: PetscCall(MatCreateVecs(ctx->A,&v0,NULL));
182: PetscCall(BVGetRandomContext(eps->V,&rnd));
183: PetscCall(VecSetRandom(v0,rnd));
184: PetscCall(VecGetArray(v0,&vinit));
185: PetscCall(PetscMalloc5(size,&nevloc,size+1,&disp,nevmax,&eigr,nevmax,&errest,nevmax*eps->n,&X));
186: mlan = PetscMin(PetscMax(5*ctx->nev,300),eps->n);
187: for (sl=rstart; sl<rend; sl++) {
188: xintv[0] = ctx->sli[sl];
189: xintv[1] = ctx->sli[sl+1];
190: xintv[2] = ctx->lmin;
191: xintv[3] = ctx->lmax;
192: PetscCall(PetscInfo(ctx->A,"Subinterval %" PetscInt_FMT ": [%.4e, %.4e]\n",sl+1,xintv[0],xintv[1]));
193: set_pol_def(&pol);
194: pol.max_deg = ctx->max_deg;
195: pol.damping = (int)ctx->damping;
196: pol.thresh_int = ctx->thresh;
197: find_pol(xintv,&pol);
198: PetscCall(PetscInfo(ctx->A,"Polynomial [type = %" PetscInt_FMT "], deg %" PetscInt_FMT ", bar %e gam %e\n",pol.type,pol.deg,pol.bar,pol.gam));
199: PetscCallEVSL(ChebLanNr,xintv,mlan,eps->tol,vinit,&pol,&nevout,&lam,&Y,&res,NULL);
200: PetscCheck(k+nevout<=nevmax,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
201: free_pol(&pol);
202: PetscCall(PetscInfo(ctx->A,"Computed %" PetscInt_FMT " eigenvalues\n",nevout));
203: PetscCall(PetscMalloc1(nevout,&ind));
204: sort_double(nevout,lam,ind);
205: for (i=0;i<nevout;i++) {
206: eigr[i+k] = lam[i];
207: errest[i+k] = res[ind[i]];
208: PetscCall(PetscArraycpy(X+(i+k)*eps->n,Y+ind[i]*eps->n,eps->n));
209: }
210: k += nevout;
211: if (lam) evsl_Free(lam);
212: if (Y) evsl_Free_device(Y);
213: if (res) evsl_Free(res);
214: PetscCall(PetscFree(ind));
215: }
216: PetscCall(VecRestoreArray(v0,&vinit));
217: PetscCall(VecDestroy(&v0));
219: /* gather eigenvalues computed by each MPI process */
220: PetscCallMPI(MPI_Allgather(&k,1,MPIU_INT,nevloc,1,MPIU_INT,PetscObjectComm((PetscObject)eps)));
221: eps->nev = nevloc[0];
222: disp[0] = 0;
223: for (i=1;i<size;i++) {
224: eps->nev += nevloc[i];
225: disp[i] = disp[i-1]+nevloc[i-1];
226: }
227: disp[size] = disp[size-1]+nevloc[size-1];
228: PetscCheck(eps->nev<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
229: PetscCallMPI(MPI_Allgatherv(eigr,k,MPIU_SCALAR,eps->eigr,nevloc,disp,MPIU_SCALAR,PetscObjectComm((PetscObject)eps)));
230: PetscCallMPI(MPI_Allgatherv(errest,k,MPIU_REAL,eps->errest,nevloc,disp,MPIU_REAL,PetscObjectComm((PetscObject)eps)));
231: eps->nconv = eps->nev;
232: eps->its = 1;
233: eps->reason = EPS_CONVERGED_TOL;
235: /* scatter computed eigenvectors and store them in eps->V */
236: PetscCall(BVCreateVec(eps->V,&w));
237: for (i=0;i<size;i++) {
238: N = (rank==i)? eps->n: 0;
239: PetscCall(VecCreateSeq(PETSC_COMM_SELF,N,&x));
240: PetscCall(VecSetFromOptions(x));
241: PetscCall(ISCreateStride(PETSC_COMM_SELF,N,0,1,&is));
242: PetscCall(VecScatterCreate(x,is,w,is,&vs));
243: PetscCall(ISDestroy(&is));
244: for (j=disp[i];j<disp[i+1];j++) {
245: PetscCall(BVGetColumn(eps->V,j,&v));
246: if (rank==i) PetscCall(VecPlaceArray(x,X+(j-disp[i])*eps->n));
247: PetscCall(VecScatterBegin(vs,x,v,INSERT_VALUES,SCATTER_FORWARD));
248: PetscCall(VecScatterEnd(vs,x,v,INSERT_VALUES,SCATTER_FORWARD));
249: if (rank==i) PetscCall(VecResetArray(x));
250: PetscCall(BVRestoreColumn(eps->V,j,&v));
251: }
252: PetscCall(VecScatterDestroy(&vs));
253: PetscCall(VecDestroy(&x));
254: }
255: PetscCall(VecDestroy(&w));
256: PetscCall(PetscFree5(nevloc,disp,eigr,errest,X));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode EPSEVSLSetSlices_EVSL(EPS eps,PetscInt nslices)
261: {
262: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
264: PetscFunctionBegin;
265: if (nslices == PETSC_DECIDE || nslices == PETSC_DEFAULT) nslices = 0;
266: else PetscCheck(nslices>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Number of slices must be 1 at least");
267: if (ctx->nslices != nslices) {
268: ctx->nslices = nslices;
269: eps->state = EPS_STATE_INITIAL;
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: EPSEVSLSetSlices - Set the number of slices in which the interval must be
276: subdivided.
278: Logically Collective
280: Input Parameters:
281: + eps - the eigensolver context
282: - nslices - the number of slices
284: Options Database Key:
285: . -eps_evsl_slices <n> - set the number of slices to n
287: Notes:
288: By default, one slice per MPI process is used. Depending on the number of
289: eigenvalues, using more slices may be beneficial, but very narrow subintervals
290: imply higher polynomial degree.
292: Level: intermediate
294: .seealso: EPSEVSLGetSlices()
295: @*/
296: PetscErrorCode EPSEVSLSetSlices(EPS eps,PetscInt nslices)
297: {
298: PetscFunctionBegin;
301: PetscTryMethod(eps,"EPSEVSLSetSlices_C",(EPS,PetscInt),(eps,nslices));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode EPSEVSLGetSlices_EVSL(EPS eps,PetscInt *nslices)
306: {
307: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
309: PetscFunctionBegin;
310: *nslices = ctx->nslices;
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@
315: EPSEVSLGetSlices - Gets the number of slices in which the interval must be
316: subdivided.
318: Not Collective
320: Input Parameter:
321: . eps - the eigensolver context
323: Output Parameter:
324: . nslices - the number of slices
326: Level: intermediate
328: .seealso: EPSEVSLSetSlices()
329: @*/
330: PetscErrorCode EPSEVSLGetSlices(EPS eps,PetscInt *nslices)
331: {
332: PetscFunctionBegin;
334: PetscAssertPointer(nslices,2);
335: PetscUseMethod(eps,"EPSEVSLGetSlices_C",(EPS,PetscInt*),(eps,nslices));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode EPSEVSLSetRange_EVSL(EPS eps,PetscReal lmin,PetscReal lmax)
340: {
341: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
343: PetscFunctionBegin;
344: PetscCheck(lmin<lmax,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be lmin<lmax");
345: if (ctx->lmin != lmin || ctx->lmax != lmax) {
346: ctx->lmin = lmin;
347: ctx->lmax = lmax;
348: eps->state = EPS_STATE_INITIAL;
349: }
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*@
354: EPSEVSLSetRange - Defines the numerical range (or field of values) of the problem,
355: that is, the interval containing all eigenvalues.
357: Logically Collective
359: Input Parameters:
360: + eps - the eigensolver context
361: . lmin - left end of the interval
362: - lmax - right end of the interval
364: Options Database Key:
365: . -eps_evsl_range <a,b> - set [a,b] as the numerical range
367: Notes:
368: The filter will be most effective if the numerical range is tight, that is, lmin
369: and lmax are good approximations to the leftmost and rightmost eigenvalues,
370: respectively. If not set by the user, an approximation is computed internally.
372: The wanted computational interval specified via EPSSetInterval() must be
373: contained in the numerical range.
375: Level: intermediate
377: .seealso: EPSEVSLGetRange(), EPSSetInterval()
378: @*/
379: PetscErrorCode EPSEVSLSetRange(EPS eps,PetscReal lmin,PetscReal lmax)
380: {
381: PetscFunctionBegin;
385: PetscTryMethod(eps,"EPSEVSLSetRange_C",(EPS,PetscReal,PetscReal),(eps,lmin,lmax));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode EPSEVSLGetRange_EVSL(EPS eps,PetscReal *lmin,PetscReal *lmax)
390: {
391: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
393: PetscFunctionBegin;
394: if (lmin) *lmin = ctx->lmin;
395: if (lmax) *lmax = ctx->lmax;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /*@
400: EPSEVSLGetRange - Gets the interval containing all eigenvalues.
402: Not Collective
404: Input Parameter:
405: . eps - the eigensolver context
407: Output Parameters:
408: + lmin - left end of the interval
409: - lmax - right end of the interval
411: Level: intermediate
413: .seealso: EPSEVSLSetRange()
414: @*/
415: PetscErrorCode EPSEVSLGetRange(EPS eps,PetscReal *lmin,PetscReal *lmax)
416: {
417: PetscFunctionBegin;
419: PetscUseMethod(eps,"EPSEVSLGetRange_C",(EPS,PetscReal*,PetscReal*),(eps,lmin,lmax));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: static PetscErrorCode EPSEVSLSetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
424: {
425: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
427: PetscFunctionBegin;
428: ctx->dos = dos;
429: if (nvec == PETSC_DETERMINE) ctx->nvec = 80;
430: else if (nvec != PETSC_CURRENT) {
431: PetscCheck(nvec>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The nvec argument must be > 0");
432: ctx->nvec = nvec;
433: }
434: switch (dos) {
435: case EPS_EVSL_DOS_KPM:
436: if (deg == PETSC_DETERMINE) ctx->deg = 300;
437: else if (deg != PETSC_CURRENT) {
438: PetscCheck(deg>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The deg argument must be > 0");
439: ctx->deg = deg;
440: }
441: break;
442: case EPS_EVSL_DOS_LANCZOS:
443: if (steps == PETSC_DETERMINE) ctx->steps = 40;
444: else if (steps != PETSC_CURRENT) {
445: PetscCheck(steps>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The steps argument must be > 0");
446: ctx->steps = steps;
447: }
448: if (npoints == PETSC_DETERMINE) ctx->npoints = 200;
449: else if (npoints != PETSC_CURRENT) {
450: PetscCheck(npoints>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npoints argument must be > 0");
451: ctx->npoints = npoints;
452: }
453: break;
454: }
455: eps->state = EPS_STATE_INITIAL;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /*@
460: EPSEVSLSetDOSParameters - Defines the parameters used for computing the
461: density of states (DOS) in the EVSL solver.
463: Logically Collective
465: Input Parameters:
466: + eps - the eigensolver context
467: . dos - DOS method, either KPM or Lanczos
468: . nvec - number of sample vectors
469: . deg - polynomial degree (KPM only)
470: . steps - number of Lanczos steps (Lanczos only)
471: - npoints - number of sample points (Lanczos only)
473: Options Database Keys:
474: + -eps_evsl_dos_method <dos> - set the DOS method, either kpm or lanczos
475: . -eps_evsl_dos_nvec <n> - set the number of sample vectors
476: . -eps_evsl_dos_degree <n> - set the polynomial degree
477: . -eps_evsl_dos_steps <n> - set the number of Lanczos steps
478: - -eps_evsl_dos_npoints <n> - set the number of sample points
480: Notes:
481: The density of states (or spectral density) can be approximated with two
482: methods, kernel polynomial method (KPM) or Lanczos. Some parameters for
483: these methods can be set by the user with this function, with some of
484: them being relevant for one of the methods only.
486: For the integer argumens, you can use PETSC_CURRENT to keep the current
487: value, and PETSC_DETERMINE to set them to a reasonable default.
489: Level: intermediate
491: .seealso: EPSEVSLGetDOSParameters()
492: @*/
493: PetscErrorCode EPSEVSLSetDOSParameters(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
494: {
495: PetscFunctionBegin;
502: PetscTryMethod(eps,"EPSEVSLSetDOSParameters_C",(EPS,EPSEVSLDOSMethod,PetscInt,PetscInt,PetscInt,PetscInt),(eps,dos,nvec,deg,steps,npoints));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode EPSEVSLGetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
507: {
508: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
510: PetscFunctionBegin;
511: if (dos) *dos = ctx->dos;
512: if (nvec) *nvec = ctx->nvec;
513: if (deg) *deg = ctx->deg;
514: if (steps) *steps = ctx->steps;
515: if (npoints) *npoints = ctx->npoints;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*@
520: EPSEVSLGetDOSParameters - Gets the parameters used for computing the
521: density of states (DOS) in the EVSL solver.
523: Not Collective
525: Input Parameter:
526: . eps - the eigensolver context
528: Output Parameters:
529: + dos - DOS method, either KPM or Lanczos
530: . nvec - number of sample vectors
531: . deg - polynomial degree (KPM only)
532: . steps - number of Lanczos steps (Lanczos only)
533: - npoints - number of sample points (Lanczos only)
535: Level: intermediate
537: .seealso: EPSEVSLSetDOSParameters()
538: @*/
539: PetscErrorCode EPSEVSLGetDOSParameters(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
540: {
541: PetscFunctionBegin;
543: PetscUseMethod(eps,"EPSEVSLGetDOSParameters_C",(EPS,EPSEVSLDOSMethod*,PetscInt*,PetscInt*,PetscInt*,PetscInt*),(eps,dos,nvec,deg,steps,npoints));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: static PetscErrorCode EPSEVSLSetPolParameters_EVSL(EPS eps,PetscInt max_deg,PetscReal thresh)
548: {
549: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
551: PetscFunctionBegin;
552: if (max_deg == PETSC_DETERMINE) ctx->max_deg = 10000;
553: else if (max_deg != PETSC_CURRENT) {
554: PetscCheck(max_deg>2,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The max_deg argument must be > 2");
555: ctx->max_deg = max_deg;
556: }
557: if (thresh == (PetscReal)PETSC_DETERMINE) ctx->thresh = 0.8;
558: else if (thresh != (PetscReal)PETSC_CURRENT) {
559: PetscCheck(thresh>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The thresh argument must be > 0.0");
560: ctx->thresh = thresh;
561: }
562: eps->state = EPS_STATE_INITIAL;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*@
567: EPSEVSLSetPolParameters - Defines the parameters used for building the
568: building the polynomial in the EVSL solver.
570: Logically Collective
572: Input Parameters:
573: + eps - the eigensolver context
574: . max_deg - maximum degree allowed for the polynomial
575: - thresh - threshold for accepting polynomial
577: Options Database Keys:
578: + -eps_evsl_pol_max_deg <d> - set maximum polynomial degree
579: - -eps_evsl_pol_thresh <t> - set the threshold
581: Note:
582: PETSC_CURRENT can be used to preserve the current value of any of the
583: arguments, and PETSC_DETERMINE to set them to a default value.
585: Level: intermediate
587: .seealso: EPSEVSLGetPolParameters()
588: @*/
589: PetscErrorCode EPSEVSLSetPolParameters(EPS eps,PetscInt max_deg,PetscReal thresh)
590: {
591: PetscFunctionBegin;
595: PetscTryMethod(eps,"EPSEVSLSetPolParameters_C",(EPS,PetscInt,PetscReal),(eps,max_deg,thresh));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode EPSEVSLGetPolParameters_EVSL(EPS eps,PetscInt *max_deg,PetscReal *thresh)
600: {
601: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
603: PetscFunctionBegin;
604: if (max_deg) *max_deg = ctx->max_deg;
605: if (thresh) *thresh = ctx->thresh;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: EPSEVSLGetPolParameters - Gets the parameters used for building the
611: polynomial in the EVSL solver.
613: Not Collective
615: Input Parameter:
616: . eps - the eigensolver context
618: Output Parameters:
619: + max_deg - the maximum degree of the polynomial
620: - thresh - the threshold
622: Level: intermediate
624: .seealso: EPSEVSLSetPolParameters()
625: @*/
626: PetscErrorCode EPSEVSLGetPolParameters(EPS eps,PetscInt *max_deg,PetscReal *thresh)
627: {
628: PetscFunctionBegin;
630: PetscUseMethod(eps,"EPSEVSLGetPolParameters_C",(EPS,PetscInt*,PetscReal*),(eps,max_deg,thresh));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode EPSEVSLSetDamping_EVSL(EPS eps,EPSEVSLDamping damping)
635: {
636: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
638: PetscFunctionBegin;
639: if (ctx->damping != damping) {
640: ctx->damping = damping;
641: eps->state = EPS_STATE_INITIAL;
642: }
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: /*@
647: EPSEVSLSetDamping - Set the type of damping to be used in EVSL.
649: Logically Collective
651: Input Parameters:
652: + eps - the eigensolver context
653: - damping - the type of damping
655: Options Database Key:
656: . -eps_evsl_damping <n> - set the type of damping
658: Notes:
659: Damping is applied when building the polynomial to be used when solving the
660: eigenproblem, and also during estimation of DOS with the KPM method.
662: Level: intermediate
664: .seealso: EPSEVSLGetDamping(), EPSEVSLSetDOSParameters()
665: @*/
666: PetscErrorCode EPSEVSLSetDamping(EPS eps,EPSEVSLDamping damping)
667: {
668: PetscFunctionBegin;
671: PetscTryMethod(eps,"EPSEVSLSetDamping_C",(EPS,EPSEVSLDamping),(eps,damping));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: static PetscErrorCode EPSEVSLGetDamping_EVSL(EPS eps,EPSEVSLDamping *damping)
676: {
677: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
679: PetscFunctionBegin;
680: *damping = ctx->damping;
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: /*@
685: EPSEVSLGetDamping - Gets the type of damping.
687: Not Collective
689: Input Parameter:
690: . eps - the eigensolver context
692: Output Parameter:
693: . damping - the type of damping
695: Level: intermediate
697: .seealso: EPSEVSLSetDamping()
698: @*/
699: PetscErrorCode EPSEVSLGetDamping(EPS eps,EPSEVSLDamping *damping)
700: {
701: PetscFunctionBegin;
703: PetscAssertPointer(damping,2);
704: PetscUseMethod(eps,"EPSEVSLGetDamping_C",(EPS,EPSEVSLDamping*),(eps,damping));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode EPSView_EVSL(EPS eps,PetscViewer viewer)
709: {
710: PetscBool isascii;
711: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
713: PetscFunctionBegin;
714: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
715: if (isascii) {
716: PetscCall(PetscViewerASCIIPrintf(viewer," numerical range = [%g,%g]\n",(double)ctx->lmin,(double)ctx->lmax));
717: PetscCall(PetscViewerASCIIPrintf(viewer," number of slices = %" PetscInt_FMT "\n",ctx->nslices));
718: PetscCall(PetscViewerASCIIPrintf(viewer," type of damping = %s\n",EPSEVSLDampings[ctx->damping]));
719: PetscCall(PetscViewerASCIIPrintf(viewer," computing DOS with %s: nvec=%" PetscInt_FMT ", ",EPSEVSLDOSMethods[ctx->dos],ctx->nvec));
720: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
721: switch (ctx->dos) {
722: case EPS_EVSL_DOS_KPM:
723: PetscCall(PetscViewerASCIIPrintf(viewer,"degree=%" PetscInt_FMT "\n",ctx->deg));
724: break;
725: case EPS_EVSL_DOS_LANCZOS:
726: PetscCall(PetscViewerASCIIPrintf(viewer,"steps=%" PetscInt_FMT ", npoints=%" PetscInt_FMT "\n",ctx->steps,ctx->npoints));
727: break;
728: }
729: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
730: PetscCall(PetscViewerASCIIPrintf(viewer," polynomial parameters: max degree = %" PetscInt_FMT ", threshold = %g\n",ctx->max_deg,(double)ctx->thresh));
731: }
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: static PetscErrorCode EPSSetFromOptions_EVSL(EPS eps,PetscOptionItems *PetscOptionsObject)
736: {
737: PetscReal array[2]={0,0},th;
738: PetscInt k,i1,i2,i3,i4;
739: PetscBool flg,flg1;
740: EPSEVSLDOSMethod dos;
741: EPSEVSLDamping damping;
742: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
744: PetscFunctionBegin;
745: PetscOptionsHeadBegin(PetscOptionsObject,"EPS EVSL Options");
747: k = 2;
748: PetscCall(PetscOptionsRealArray("-eps_evsl_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","EPSEVSLSetRange",array,&k,&flg));
749: if (flg) {
750: PetscCheck(k>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_SIZ,"Must pass two values in -eps_evsl_range (comma-separated without spaces)");
751: PetscCall(EPSEVSLSetRange(eps,array[0],array[1]));
752: }
754: PetscCall(PetscOptionsInt("-eps_evsl_slices","Number of slices","EPSEVSLSetSlices",ctx->nslices,&i1,&flg));
755: if (flg) PetscCall(EPSEVSLSetSlices(eps,i1));
757: PetscCall(PetscOptionsEnum("-eps_evsl_damping","Type of damping","EPSEVSLSetDamping",EPSEVSLDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg));
758: if (flg) PetscCall(EPSEVSLSetDamping(eps,damping));
760: PetscCall(EPSEVSLGetDOSParameters(eps,&dos,&i1,&i2,&i3,&i4));
761: PetscCall(PetscOptionsEnum("-eps_evsl_dos_method","Method to compute the DOS","EPSEVSLSetDOSParameters",EPSEVSLDOSMethods,(PetscEnum)ctx->dos,(PetscEnum*)&dos,&flg));
762: PetscCall(PetscOptionsInt("-eps_evsl_dos_nvec","Number of sample vectors for DOS","EPSEVSLSetDOSParameters",i1,&i1,&flg1));
763: if (flg1) flg = PETSC_TRUE;
764: PetscCall(PetscOptionsInt("-eps_evsl_dos_degree","Polynomial degree used for DOS","EPSEVSLSetDOSParameters",i2,&i2,&flg1));
765: if (flg1) flg = PETSC_TRUE;
766: PetscCall(PetscOptionsInt("-eps_evsl_dos_steps","Number of Lanczos steps in DOS","EPSEVSLSetDOSParameters",i3,&i3,&flg1));
767: if (flg1) flg = PETSC_TRUE;
768: PetscCall(PetscOptionsInt("-eps_evsl_dos_npoints","Number of sample points used for DOS","EPSEVSLSetDOSParameters",i4,&i4,&flg1));
769: if (flg || flg1) PetscCall(EPSEVSLSetDOSParameters(eps,dos,i1,i2,i3,i4));
771: PetscCall(EPSEVSLGetPolParameters(eps,&i1,&th));
772: PetscCall(PetscOptionsInt("-eps_evsl_pol_max_deg","Maximum degree allowed for the polynomial","EPSEVSLSetPolParameters",i1,&i1,&flg));
773: PetscCall(PetscOptionsReal("-eps_evsl_pol_threshold","Threshold for accepting polynomial","EPSEVSLSetPolParameters",th,&th,&flg1));
774: if (flg || flg1) PetscCall(EPSEVSLSetPolParameters(eps,i1,th));
776: PetscOptionsHeadEnd();
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: static PetscErrorCode EPSDestroy_EVSL(EPS eps)
781: {
782: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
784: PetscFunctionBegin;
785: if (ctx->initialized) EVSLFinish();
786: PetscCall(PetscLayoutDestroy(&ctx->map));
787: PetscCall(PetscFree(ctx->sli));
788: PetscCall(PetscFree(eps->data));
789: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",NULL));
790: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",NULL));
791: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",NULL));
792: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",NULL));
793: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",NULL));
794: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",NULL));
795: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",NULL));
796: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",NULL));
797: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",NULL));
798: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",NULL));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: static PetscErrorCode EPSReset_EVSL(EPS eps)
803: {
804: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
806: PetscFunctionBegin;
807: PetscCall(MatDestroy(&ctx->A));
808: PetscCall(VecDestroy(&ctx->x));
809: PetscCall(VecDestroy(&ctx->y));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: SLEPC_EXTERN PetscErrorCode EPSCreate_EVSL(EPS eps)
814: {
815: EPS_EVSL *ctx;
817: PetscFunctionBegin;
818: PetscCall(PetscNew(&ctx));
819: eps->data = (void*)ctx;
821: ctx->nslices = 0;
822: ctx->lmin = PETSC_MIN_REAL;
823: ctx->lmax = PETSC_MAX_REAL;
824: ctx->dos = EPS_EVSL_DOS_KPM;
825: ctx->nvec = 80;
826: ctx->deg = 300;
827: ctx->steps = 40;
828: ctx->npoints = 200;
829: ctx->max_deg = 10000;
830: ctx->thresh = 0.8;
831: ctx->damping = EPS_EVSL_DAMPING_SIGMA;
833: eps->categ = EPS_CATEGORY_OTHER;
835: eps->ops->solve = EPSSolve_EVSL;
836: eps->ops->setup = EPSSetUp_EVSL;
837: eps->ops->setupsort = EPSSetUpSort_Basic;
838: eps->ops->setfromoptions = EPSSetFromOptions_EVSL;
839: eps->ops->destroy = EPSDestroy_EVSL;
840: eps->ops->reset = EPSReset_EVSL;
841: eps->ops->view = EPSView_EVSL;
842: eps->ops->backtransform = EPSBackTransform_Default;
843: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
845: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",EPSEVSLSetRange_EVSL));
846: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",EPSEVSLGetRange_EVSL));
847: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",EPSEVSLSetSlices_EVSL));
848: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",EPSEVSLGetSlices_EVSL));
849: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",EPSEVSLSetDOSParameters_EVSL));
850: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",EPSEVSLGetDOSParameters_EVSL));
851: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",EPSEVSLSetPolParameters_EVSL));
852: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",EPSEVSLGetPolParameters_EVSL));
853: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",EPSEVSLSetDamping_EVSL));
854: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",EPSEVSLGetDamping_EVSL));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }