Actual source code: evsl.c
slepc-3.17.2 2022-08-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: PetscStackPush(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 {
93: MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&ctx->A);
94: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->A);
95: }
96: SetAMatvec(eps->n,&AMatvec_EVSL,(void*)ctx);
97: if (!ctx->x) {
98: MatCreateVecsEmpty(ctx->A,&ctx->x,&ctx->y);
99: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->x);
100: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->y);
101: }
102: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
103: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE);
105: if (!eps->which) eps->which=EPS_ALL;
108: /* estimate numerical range */
109: if (ctx->estimrange || ctx->lmin == PETSC_MIN_REAL || ctx->lmax == PETSC_MAX_REAL) {
110: MatCreateVecs(ctx->A,&v0,NULL);
111: if (!eps->V) EPSGetBV(eps,&eps->V);
112: BVGetRandomContext(eps->V,&rnd);
113: VecSetRandom(v0,rnd);
114: VecGetArray(v0,&vinit);
115: LanTrbounds,50,200,eps->tol,vinit,1,&ctx->lmin,&ctx->lmax,NULL;
116: VecRestoreArray(v0,&vinit);
117: VecDestroy(&v0);
118: ctx->estimrange = PETSC_TRUE; /* estimate if called again with another matrix */
119: }
121: xintv[0] = eps->inta;
122: xintv[1] = eps->intb;
123: xintv[2] = ctx->lmin;
124: xintv[3] = ctx->lmax;
126: /* estimate number of eigenvalues in the interval */
127: switch (ctx->dos) {
128: case EPS_EVSL_DOS_KPM:
129: PetscMalloc1(ctx->deg+1,&mu);
130: if (!rank) kpmdos,ctx->deg,(int)ctx->damping,ctx->nvec,xintv,mu,&ecount;
131: MPI_Bcast(mu,ctx->deg+1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
132: break;
133: case EPS_EVSL_DOS_LANCZOS:
134: PetscMalloc2(ctx->npoints,&xdos,ctx->npoints,&ydos);
135: if (!rank) LanDos,ctx->nvec,PetscMin(ctx->steps,eps->n/2),ctx->npoints,xdos,ydos,&ecount,xintv;
136: MPI_Bcast(xdos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
137: MPI_Bcast(ydos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
138: break;
139: default:
140: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid DOS method");
141: }
142: MPI_Bcast(&ecount,1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));
144: PetscInfo(eps,"Estimated eigenvalue count in the interval: %g\n",ecount);
145: eps->ncv = (PetscInt)PetscCeilReal(1.5*ecount);
147: /* slice the spectrum */
148: PetscFree(ctx->sli);
149: PetscMalloc1(ctx->nslices+1,&ctx->sli);
150: if (ctx->dos == EPS_EVSL_DOS_KPM) {
151: spslicer,ctx->sli,mu,ctx->deg,xintv,ctx->nslices,10*(PetscInt)ecount;
152: PetscFree(mu);
153: } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
154: spslicer2(xdos,ydos,ctx->nslices,ctx->npoints,ctx->sli);
155: PetscFree2(xdos,ydos);
156: }
158: /* approximate number of eigenvalues wanted in each slice */
159: ctx->nev = (PetscInt)(1.0 + ecount/(PetscReal)ctx->nslices) + 2;
161: if (eps->mpd!=PETSC_DEFAULT) PetscInfo(eps,"Warning: parameter mpd ignored\n");
162: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
163: EPSAllocateSolution(eps,0);
164: PetscFunctionReturn(0);
165: }
167: PetscErrorCode EPSSolve_EVSL(EPS eps)
168: {
169: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
170: PetscInt i,j,k=0,sl,mlan,nevout,*ind,nevmax,rstart,rend,*nevloc,*disp,N;
171: PetscReal *res,xintv[4],*errest;
172: PetscScalar *lam,*X,*Y,*vinit,*eigr;
173: PetscMPIInt size,rank;
174: PetscRandom rnd;
175: Vec v,w,v0,x;
176: VecScatter vs;
177: IS is;
178: polparams pol;
180: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
181: MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank);
182: PetscLayoutGetRange(ctx->map,&rstart,&rend);
183: nevmax = (rend-rstart)*ctx->nev;
184: MatCreateVecs(ctx->A,&v0,NULL);
185: BVGetRandomContext(eps->V,&rnd);
186: VecSetRandom(v0,rnd);
187: VecGetArray(v0,&vinit);
188: PetscMalloc5(size,&nevloc,size+1,&disp,nevmax,&eigr,nevmax,&errest,nevmax*eps->n,&X);
189: mlan = PetscMin(PetscMax(5*ctx->nev,300),eps->n);
190: for (sl=rstart; sl<rend; sl++) {
191: xintv[0] = ctx->sli[sl];
192: xintv[1] = ctx->sli[sl+1];
193: xintv[2] = ctx->lmin;
194: xintv[3] = ctx->lmax;
195: PetscInfo(ctx->A,"Subinterval %" PetscInt_FMT ": [%.4e, %.4e]\n",sl+1,xintv[0],xintv[1]);
196: set_pol_def(&pol);
197: pol.max_deg = ctx->max_deg;
198: pol.damping = (int)ctx->damping;
199: pol.thresh_int = ctx->thresh;
200: find_pol(xintv,&pol);
201: PetscInfo(ctx->A,"Polynomial [type = %" PetscInt_FMT "], deg %" PetscInt_FMT ", bar %e gam %e\n",pol.type,pol.deg,pol.bar,pol.gam);
202: ChebLanNr,xintv,mlan,eps->tol,vinit,&pol,&nevout,&lam,&Y,&res,NULL;
204: free_pol(&pol);
205: PetscInfo(ctx->A,"Computed %" PetscInt_FMT " eigenvalues\n",nevout);
206: PetscMalloc1(nevout,&ind);
207: sort_double(nevout,lam,ind);
208: for (i=0;i<nevout;i++) {
209: eigr[i+k] = lam[i];
210: errest[i+k] = res[ind[i]];
211: PetscArraycpy(X+(i+k)*eps->n,Y+ind[i]*eps->n,eps->n);
212: }
213: k += nevout;
214: if (lam) evsl_Free(lam);
215: if (Y) evsl_Free_device(Y);
216: if (res) evsl_Free(res);
217: PetscFree(ind);
218: }
219: VecRestoreArray(v0,&vinit);
220: VecDestroy(&v0);
222: /* gather eigenvalues computed by each MPI process */
223: MPI_Allgather(&k,1,MPIU_INT,nevloc,1,MPIU_INT,PetscObjectComm((PetscObject)eps));
224: eps->nev = nevloc[0];
225: disp[0] = 0;
226: for (i=1;i<size;i++) {
227: eps->nev += nevloc[i];
228: disp[i] = disp[i-1]+nevloc[i-1];
229: }
230: disp[size] = disp[size-1]+nevloc[size-1];
232: MPI_Allgatherv(eigr,k,MPIU_SCALAR,eps->eigr,nevloc,disp,MPIU_SCALAR,PetscObjectComm((PetscObject)eps));
233: MPI_Allgatherv(errest,k,MPIU_REAL,eps->errest,nevloc,disp,MPIU_REAL,PetscObjectComm((PetscObject)eps));
234: eps->nconv = eps->nev;
235: eps->its = 1;
236: eps->reason = EPS_CONVERGED_TOL;
238: /* scatter computed eigenvectors and store them in eps->V */
239: BVCreateVec(eps->V,&w);
240: for (i=0;i<size;i++) {
241: N = (rank==i)? eps->n: 0;
242: VecCreateSeq(PETSC_COMM_SELF,N,&x);
243: VecSetFromOptions(x);
244: ISCreateStride(PETSC_COMM_SELF,N,0,1,&is);
245: VecScatterCreate(x,is,w,is,&vs);
246: ISDestroy(&is);
247: for (j=disp[i];j<disp[i+1];j++) {
248: BVGetColumn(eps->V,j,&v);
249: if (rank==i) VecPlaceArray(x,X+(j-disp[i])*eps->n);
250: VecScatterBegin(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
251: VecScatterEnd(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
252: if (rank==i) VecResetArray(x);
253: BVRestoreColumn(eps->V,j,&v);
254: }
255: VecScatterDestroy(&vs);
256: VecDestroy(&x);
257: }
258: VecDestroy(&w);
259: PetscFree5(nevloc,disp,eigr,errest,X);
260: PetscFunctionReturn(0);
261: }
263: static PetscErrorCode EPSEVSLSetSlices_EVSL(EPS eps,PetscInt nslices)
264: {
265: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
267: if (nslices == PETSC_DECIDE || nslices == PETSC_DEFAULT) nslices = 0;
269: if (ctx->nslices != nslices) {
270: ctx->nslices = nslices;
271: eps->state = EPS_STATE_INITIAL;
272: }
273: PetscFunctionReturn(0);
274: }
276: /*@
277: EPSEVSLSetSlices - Set the number of slices in which the interval must be
278: subdivided.
280: Logically Collective on eps
282: Input Parameters:
283: + eps - the eigensolver context
284: - nslices - the number of slices
286: Options Database Key:
287: . -eps_evsl_slices <n> - set the number of slices to n
289: Notes:
290: By default, one slice per MPI process is used. Depending on the number of
291: eigenvalues, using more slices may be beneficial, but very narrow subintervals
292: imply higher polynomial degree.
294: Level: intermediate
296: .seealso: EPSEVSLGetSlices()
297: @*/
298: PetscErrorCode EPSEVSLSetSlices(EPS eps,PetscInt nslices)
299: {
302: PetscTryMethod(eps,"EPSEVSLSetSlices_C",(EPS,PetscInt),(eps,nslices));
303: PetscFunctionReturn(0);
304: }
306: static PetscErrorCode EPSEVSLGetSlices_EVSL(EPS eps,PetscInt *nslices)
307: {
308: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
310: *nslices = ctx->nslices;
311: PetscFunctionReturn(0);
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: {
334: PetscUseMethod(eps,"EPSEVSLGetSlices_C",(EPS,PetscInt*),(eps,nslices));
335: PetscFunctionReturn(0);
336: }
338: static PetscErrorCode EPSEVSLSetRange_EVSL(EPS eps,PetscReal lmin,PetscReal lmax)
339: {
340: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
343: if (ctx->lmin != lmin || ctx->lmax != lmax) {
344: ctx->lmin = lmin;
345: ctx->lmax = lmax;
346: eps->state = EPS_STATE_INITIAL;
347: }
348: PetscFunctionReturn(0);
349: }
351: /*@
352: EPSEVSLSetRange - Defines the numerical range (or field of values) of the problem,
353: that is, the interval containing all eigenvalues.
355: Logically Collective on eps
357: Input Parameters:
358: + eps - the eigensolver context
359: . lmin - left end of the interval
360: - lmax - right end of the interval
362: Options Database Key:
363: . -eps_evsl_range <a,b> - set [a,b] as the numerical range
365: Notes:
366: The filter will be most effective if the numerical range is tight, that is, lmin
367: and lmax are good approximations to the leftmost and rightmost eigenvalues,
368: respectively. If not set by the user, an approximation is computed internally.
370: The wanted computational interval specified via EPSSetInterval() must be
371: contained in the numerical range.
373: Level: intermediate
375: .seealso: EPSEVSLGetRange(), EPSSetInterval()
376: @*/
377: PetscErrorCode EPSEVSLSetRange(EPS eps,PetscReal lmin,PetscReal lmax)
378: {
382: PetscTryMethod(eps,"EPSEVSLSetRange_C",(EPS,PetscReal,PetscReal),(eps,lmin,lmax));
383: PetscFunctionReturn(0);
384: }
386: static PetscErrorCode EPSEVSLGetRange_EVSL(EPS eps,PetscReal *lmin,PetscReal *lmax)
387: {
388: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
390: if (lmin) *lmin = ctx->lmin;
391: if (lmax) *lmax = ctx->lmax;
392: PetscFunctionReturn(0);
393: }
395: /*@
396: EPSEVSLGetRange - Gets the interval containing all eigenvalues.
398: Not Collective
400: Input Parameter:
401: . eps - the eigensolver context
403: Output Parameters:
404: + lmin - left end of the interval
405: - lmax - right end of the interval
407: Level: intermediate
409: .seealso: EPSEVSLSetRange()
410: @*/
411: PetscErrorCode EPSEVSLGetRange(EPS eps,PetscReal *lmin,PetscReal *lmax)
412: {
414: PetscUseMethod(eps,"EPSEVSLGetRange_C",(EPS,PetscReal*,PetscReal*),(eps,lmin,lmax));
415: PetscFunctionReturn(0);
416: }
418: static PetscErrorCode EPSEVSLSetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
419: {
420: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
422: ctx->dos = dos;
423: if (nvec == PETSC_DECIDE || nvec == PETSC_DEFAULT) ctx->nvec = 80;
424: else {
426: ctx->nvec = nvec;
427: }
428: switch (dos) {
429: case EPS_EVSL_DOS_KPM:
430: if (deg == PETSC_DECIDE || deg == PETSC_DEFAULT) ctx->deg = 300;
431: else {
433: ctx->deg = deg;
434: }
435: break;
436: case EPS_EVSL_DOS_LANCZOS:
437: if (steps == PETSC_DECIDE || steps == PETSC_DEFAULT) ctx->steps = 40;
438: else {
440: ctx->steps = steps;
441: }
442: if (npoints == PETSC_DECIDE || npoints == PETSC_DEFAULT) ctx->npoints = 200;
443: else {
445: ctx->npoints = npoints;
446: }
447: break;
448: }
449: eps->state = EPS_STATE_INITIAL;
450: PetscFunctionReturn(0);
451: }
453: /*@
454: EPSEVSLSetDOSParameters - Defines the parameters used for computing the
455: density of states (DOS) in the EVSL solver.
457: Logically Collective on eps
459: Input Parameters:
460: + eps - the eigensolver context
461: . dos - DOS method, either KPM or Lanczos
462: . nvec - number of sample vectors
463: . deg - polynomial degree (KPM only)
464: . steps - number of Lanczos steps (Lanczos only)
465: - npoints - number of sample points (Lanczos only)
467: Options Database Keys:
468: + -eps_evsl_dos_method <dos> - set the DOS method, either kpm or lanczos
469: . -eps_evsl_dos_nvec <n> - set the number of sample vectors
470: . -eps_evsl_dos_degree <n> - set the polynomial degree
471: . -eps_evsl_dos_steps <n> - set the number of Lanczos steps
472: - -eps_evsl_dos_npoints <n> - set the number of sample points
474: Notes:
475: The density of states (or spectral density) can be approximated with two
476: methods, kernel polynomial method (KPM) or Lanczos. Some parameters for
477: these methods can be set by the user with this function, with some of
478: them being relevant for one of the methods only.
480: Level: intermediate
482: .seealso: EPSEVSLGetDOSParameters()
483: @*/
484: PetscErrorCode EPSEVSLSetDOSParameters(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
485: {
492: PetscTryMethod(eps,"EPSEVSLSetDOSParameters_C",(EPS,EPSEVSLDOSMethod,PetscInt,PetscInt,PetscInt,PetscInt),(eps,dos,nvec,deg,steps,npoints));
493: PetscFunctionReturn(0);
494: }
496: static PetscErrorCode EPSEVSLGetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
497: {
498: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
500: if (dos) *dos = ctx->dos;
501: if (nvec) *nvec = ctx->nvec;
502: if (deg) *deg = ctx->deg;
503: if (steps) *steps = ctx->steps;
504: if (npoints) *npoints = ctx->npoints;
505: PetscFunctionReturn(0);
506: }
508: /*@
509: EPSEVSLGetDOSParameters - Gets the parameters used for computing the
510: density of states (DOS) in the EVSL solver.
512: Not Collective
514: Input Parameter:
515: . eps - the eigensolver context
517: Output Parameters:
518: + dos - DOS method, either KPM or Lanczos
519: . nvec - number of sample vectors
520: . deg - polynomial degree (KPM only)
521: . steps - number of Lanczos steps (Lanczos only)
522: - npoints - number of sample points (Lanczos only)
524: Level: intermediate
526: .seealso: EPSEVSLSetDOSParameters()
527: @*/
528: PetscErrorCode EPSEVSLGetDOSParameters(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
529: {
531: PetscUseMethod(eps,"EPSEVSLGetDOSParameters_C",(EPS,EPSEVSLDOSMethod*,PetscInt*,PetscInt*,PetscInt*,PetscInt*),(eps,dos,nvec,deg,steps,npoints));
532: PetscFunctionReturn(0);
533: }
535: static PetscErrorCode EPSEVSLSetPolParameters_EVSL(EPS eps,PetscInt max_deg,PetscReal thresh)
536: {
537: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
539: if (max_deg == PETSC_DECIDE || max_deg == PETSC_DEFAULT) ctx->max_deg = 10000;
540: else {
542: ctx->max_deg = max_deg;
543: }
544: if (thresh == PETSC_DECIDE || thresh == PETSC_DEFAULT) ctx->thresh = 0.8;
545: else {
547: ctx->thresh = thresh;
548: }
549: eps->state = EPS_STATE_INITIAL;
550: PetscFunctionReturn(0);
551: }
553: /*@
554: EPSEVSLSetPolParameters - Defines the parameters used for building the
555: building the polynomial in the EVSL solver.
557: Logically Collective on eps
559: Input Parameters:
560: + eps - the eigensolver context
561: . max_deg - maximum degree allowed for the polynomial
562: - thresh - threshold for accepting polynomial
564: Options Database Keys:
565: + -eps_evsl_pol_max_deg <d> - set maximum polynomial degree
566: - -eps_evsl_pol_thresh <t> - set the threshold
568: Level: intermediate
570: .seealso: EPSEVSLGetPolParameters()
571: @*/
572: PetscErrorCode EPSEVSLSetPolParameters(EPS eps,PetscInt max_deg,PetscReal thresh)
573: {
577: PetscTryMethod(eps,"EPSEVSLSetPolParameters_C",(EPS,PetscInt,PetscReal),(eps,max_deg,thresh));
578: PetscFunctionReturn(0);
579: }
581: static PetscErrorCode EPSEVSLGetPolParameters_EVSL(EPS eps,PetscInt *max_deg,PetscReal *thresh)
582: {
583: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
585: if (max_deg) *max_deg = ctx->max_deg;
586: if (thresh) *thresh = ctx->thresh;
587: PetscFunctionReturn(0);
588: }
590: /*@
591: EPSEVSLGetPolParameters - Gets the parameters used for building the
592: polynomial in the EVSL solver.
594: Not Collective
596: Input Parameter:
597: . eps - the eigensolver context
599: Output Parameters:
600: + max_deg - the maximum degree of the polynomial
601: - thresh - the threshold
603: Level: intermediate
605: .seealso: EPSEVSLSetPolParameters()
606: @*/
607: PetscErrorCode EPSEVSLGetPolParameters(EPS eps,PetscInt *max_deg,PetscReal *thresh)
608: {
610: PetscUseMethod(eps,"EPSEVSLGetPolParameters_C",(EPS,PetscInt*,PetscReal*),(eps,max_deg,thresh));
611: PetscFunctionReturn(0);
612: }
614: static PetscErrorCode EPSEVSLSetDamping_EVSL(EPS eps,EPSEVSLDamping damping)
615: {
616: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
618: if (ctx->damping != damping) {
619: ctx->damping = damping;
620: eps->state = EPS_STATE_INITIAL;
621: }
622: PetscFunctionReturn(0);
623: }
625: /*@
626: EPSEVSLSetDamping - Set the type of damping to be used in EVSL.
628: Logically Collective on eps
630: Input Parameters:
631: + eps - the eigensolver context
632: - damping - the type of damping
634: Options Database Key:
635: . -eps_evsl_damping <n> - set the type of damping
637: Notes:
638: Damping is applied when building the polynomial to be used when solving the
639: eigenproblem, and also during estimation of DOS with the KPM method.
641: Level: intermediate
643: .seealso: EPSEVSLGetDamping(), EPSEVSLSetDOSParameters()
644: @*/
645: PetscErrorCode EPSEVSLSetDamping(EPS eps,EPSEVSLDamping damping)
646: {
649: PetscTryMethod(eps,"EPSEVSLSetDamping_C",(EPS,EPSEVSLDamping),(eps,damping));
650: PetscFunctionReturn(0);
651: }
653: static PetscErrorCode EPSEVSLGetDamping_EVSL(EPS eps,EPSEVSLDamping *damping)
654: {
655: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
657: *damping = ctx->damping;
658: PetscFunctionReturn(0);
659: }
661: /*@
662: EPSEVSLGetDamping - Gets the type of damping.
664: Not Collective
666: Input Parameter:
667: . eps - the eigensolver context
669: Output Parameter:
670: . damping - the type of damping
672: Level: intermediate
674: .seealso: EPSEVSLSetDamping()
675: @*/
676: PetscErrorCode EPSEVSLGetDamping(EPS eps,EPSEVSLDamping *damping)
677: {
680: PetscUseMethod(eps,"EPSEVSLGetDamping_C",(EPS,EPSEVSLDamping*),(eps,damping));
681: PetscFunctionReturn(0);
682: }
684: PetscErrorCode EPSView_EVSL(EPS eps,PetscViewer viewer)
685: {
686: PetscBool isascii;
687: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
689: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
690: if (isascii) {
691: PetscViewerASCIIPrintf(viewer," numerical range = [%g,%g]\n",(double)ctx->lmin,(double)ctx->lmax);
692: PetscViewerASCIIPrintf(viewer," number of slices = %" PetscInt_FMT "\n",ctx->nslices);
693: PetscViewerASCIIPrintf(viewer," type of damping = %s\n",EPSEVSLDampings[ctx->damping]);
694: PetscViewerASCIIPrintf(viewer," computing DOS with %s: nvec=%" PetscInt_FMT ", ",EPSEVSLDOSMethods[ctx->dos],ctx->nvec);
695: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
696: switch (ctx->dos) {
697: case EPS_EVSL_DOS_KPM:
698: PetscViewerASCIIPrintf(viewer,"degree=%" PetscInt_FMT "\n",ctx->deg);
699: break;
700: case EPS_EVSL_DOS_LANCZOS:
701: PetscViewerASCIIPrintf(viewer,"steps=%" PetscInt_FMT ", npoints=%" PetscInt_FMT "\n",ctx->steps,ctx->npoints);
702: break;
703: }
704: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
705: PetscViewerASCIIPrintf(viewer," polynomial parameters: max degree = %" PetscInt_FMT ", threshold = %g\n",ctx->max_deg,(double)ctx->thresh);
706: }
707: PetscFunctionReturn(0);
708: }
710: PetscErrorCode EPSSetFromOptions_EVSL(PetscOptionItems *PetscOptionsObject,EPS eps)
711: {
712: PetscReal array[2]={0,0},th;
713: PetscInt k,i1,i2,i3,i4;
714: PetscBool flg,flg1;
715: EPSEVSLDOSMethod dos;
716: EPSEVSLDamping damping;
717: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
719: PetscOptionsHead(PetscOptionsObject,"EPS EVSL Options");
721: k = 2;
722: PetscOptionsRealArray("-eps_evsl_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","EPSEVSLSetRange",array,&k,&flg);
723: if (flg) {
725: EPSEVSLSetRange(eps,array[0],array[1]);
726: }
728: PetscOptionsInt("-eps_evsl_slices","Number of slices","EPSEVSLSetSlices",ctx->nslices,&i1,&flg);
729: if (flg) EPSEVSLSetSlices(eps,i1);
731: PetscOptionsEnum("-eps_evsl_damping","Type of damping","EPSEVSLSetDamping",EPSEVSLDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg);
732: if (flg) EPSEVSLSetDamping(eps,damping);
734: EPSEVSLGetDOSParameters(eps,&dos,&i1,&i2,&i3,&i4);
735: PetscOptionsEnum("-eps_evsl_dos_method","Method to compute the DOS","EPSEVSLSetDOSParameters",EPSEVSLDOSMethods,(PetscEnum)ctx->dos,(PetscEnum*)&dos,&flg);
736: PetscOptionsInt("-eps_evsl_dos_nvec","Number of sample vectors for DOS","EPSEVSLSetDOSParameters",i1,&i1,&flg1);
737: if (flg1) flg = PETSC_TRUE;
738: PetscOptionsInt("-eps_evsl_dos_degree","Polynomial degree used for DOS","EPSEVSLSetDOSParameters",i2,&i2,&flg1);
739: if (flg1) flg = PETSC_TRUE;
740: PetscOptionsInt("-eps_evsl_dos_steps","Number of Lanczos steps in DOS","EPSEVSLSetDOSParameters",i3,&i3,&flg1);
741: if (flg1) flg = PETSC_TRUE;
742: PetscOptionsInt("-eps_evsl_dos_npoints","Number of sample points used for DOS","EPSEVSLSetDOSParameters",i4,&i4,&flg1);
743: if (flg || flg1) EPSEVSLSetDOSParameters(eps,dos,i1,i2,i3,i4);
745: EPSEVSLGetPolParameters(eps,&i1,&th);
746: PetscOptionsInt("-eps_evsl_pol_max_deg","Maximum degree allowed for the polynomial","EPSEVSLSetPolParameters",i1,&i1,&flg);
747: PetscOptionsReal("-eps_evsl_pol_threshold","Threshold for accepting polynomial","EPSEVSLSetPolParameters",th,&th,&flg1);
748: if (flg || flg1) EPSEVSLSetPolParameters(eps,i1,th);
750: PetscOptionsTail();
751: PetscFunctionReturn(0);
752: }
754: PetscErrorCode EPSDestroy_EVSL(EPS eps)
755: {
756: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
758: if (ctx->initialized) EVSLFinish();
759: PetscLayoutDestroy(&ctx->map);
760: PetscFree(ctx->sli);
761: PetscFree(eps->data);
762: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",NULL);
763: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",NULL);
764: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",NULL);
765: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",NULL);
766: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",NULL);
767: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",NULL);
768: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",NULL);
769: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",NULL);
770: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",NULL);
771: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",NULL);
772: PetscFunctionReturn(0);
773: }
775: PetscErrorCode EPSReset_EVSL(EPS eps)
776: {
777: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
779: MatDestroy(&ctx->A);
780: VecDestroy(&ctx->x);
781: VecDestroy(&ctx->y);
782: PetscFunctionReturn(0);
783: }
785: SLEPC_EXTERN PetscErrorCode EPSCreate_EVSL(EPS eps)
786: {
787: EPS_EVSL *ctx;
789: PetscNewLog(eps,&ctx);
790: eps->data = (void*)ctx;
792: ctx->nslices = 0;
793: ctx->lmin = PETSC_MIN_REAL;
794: ctx->lmax = PETSC_MAX_REAL;
795: ctx->dos = EPS_EVSL_DOS_KPM;
796: ctx->nvec = 80;
797: ctx->deg = 300;
798: ctx->steps = 40;
799: ctx->npoints = 200;
800: ctx->max_deg = 10000;
801: ctx->thresh = 0.8;
802: ctx->damping = EPS_EVSL_DAMPING_SIGMA;
804: eps->categ = EPS_CATEGORY_OTHER;
806: eps->ops->solve = EPSSolve_EVSL;
807: eps->ops->setup = EPSSetUp_EVSL;
808: eps->ops->setupsort = EPSSetUpSort_Basic;
809: eps->ops->setfromoptions = EPSSetFromOptions_EVSL;
810: eps->ops->destroy = EPSDestroy_EVSL;
811: eps->ops->reset = EPSReset_EVSL;
812: eps->ops->view = EPSView_EVSL;
813: eps->ops->backtransform = EPSBackTransform_Default;
814: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
816: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",EPSEVSLSetRange_EVSL);
817: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",EPSEVSLGetRange_EVSL);
818: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",EPSEVSLSetSlices_EVSL);
819: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",EPSEVSLGetSlices_EVSL);
820: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",EPSEVSLSetDOSParameters_EVSL);
821: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",EPSEVSLGetDOSParameters_EVSL);
822: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",EPSEVSLSetPolParameters_EVSL);
823: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",EPSEVSLGetPolParameters_EVSL);
824: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",EPSEVSLSetDamping_EVSL);
825: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",EPSEVSLGetDamping_EVSL);
826: PetscFunctionReturn(0);
827: }