Actual source code: evsl.c
slepc-3.15.0 2021-03-31
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: typedef struct {
18: PetscBool initialized;
19: Mat A; /* problem matrix */
20: Vec x,y; /* auxiliary vectors */
21: PetscReal *sli; /* slice bounds */
22: PetscInt nev; /* approximate number of wanted eigenvalues in each slice */
23: PetscLayout map; /* used to distribute slices among MPI processes */
24: PetscBool estimrange; /* the filter range was not set by the user */
25: /* user parameters */
26: PetscInt nslices; /* number of slices */
27: PetscReal lmin,lmax; /* numerical range (min and max eigenvalue) */
28: EPSEVSLDOSMethod dos; /* DOS method, either KPM or Lanczos */
29: PetscInt nvec; /* number of sample vectors used for DOS */
30: PetscInt deg; /* polynomial degree used for DOS (KPM only) */
31: PetscInt steps; /* number of Lanczos steps used for DOS (Lanczos only) */
32: PetscInt npoints; /* number of sample points used for DOS (Lanczos only) */
33: PetscInt max_deg; /* maximum degree allowed for the polynomial */
34: PetscReal thresh; /* threshold for accepting polynomial */
35: EPSEVSLDamping damping; /* type of damping (for polynomial and for DOS-KPM) */
36: } EPS_EVSL;
38: static void AMatvec_EVSL(double *xa,double *ya,void *data)
39: {
41: EPS_EVSL *ctx = (EPS_EVSL*)data;
42: Vec x = ctx->x,y = ctx->y;
43: Mat A = ctx->A;
46: VecPlaceArray(x,(PetscScalar*)xa);CHKERRABORT(PetscObjectComm((PetscObject)A),ierr);
47: VecPlaceArray(y,(PetscScalar*)ya);CHKERRABORT(PetscObjectComm((PetscObject)A),ierr);
48: MatMult(A,x,y);CHKERRABORT(PetscObjectComm((PetscObject)A),ierr);
49: VecResetArray(x);CHKERRABORT(PetscObjectComm((PetscObject)A),ierr);
50: VecResetArray(y);CHKERRABORT(PetscObjectComm((PetscObject)A),ierr);
51: PetscFunctionReturnVoid();
52: }
54: PetscErrorCode EPSSetUp_EVSL(EPS eps)
55: {
57: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
58: PetscMPIInt size,rank;
59: PetscBool isshift;
60: PetscScalar *vinit;
61: PetscReal *mu,ecount,xintv[4],*xdos,*ydos;
62: Vec v0;
63: Mat A;
64: PetscRandom rnd;
67: EPSCheckStandard(eps);
68: EPSCheckHermitian(eps);
69: PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift);
70: if (!isshift) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
72: if (ctx->initialized) EVSLFinish();
73: EVSLStart();
74: ctx->initialized=PETSC_TRUE;
76: /* get number of slices per process */
77: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);CHKERRMPI(ierr);
78: MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank);CHKERRMPI(ierr);
79: if (!ctx->nslices) ctx->nslices = size;
80: PetscLayoutDestroy(&ctx->map);
81: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)eps),PETSC_DECIDE,ctx->nslices,1,&ctx->map);
83: /* get matrix and prepare auxiliary vectors */
84: MatDestroy(&ctx->A);
85: STGetMatrix(eps->st,0,&A);
86: if (size==1) {
87: PetscObjectReference((PetscObject)A);
88: ctx->A = A;
89: } else {
90: MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&ctx->A);
91: }
92: SetAMatvec(eps->n,&AMatvec_EVSL,(void*)ctx);
93: if (!ctx->x) {
94: MatCreateVecsEmpty(ctx->A,&ctx->x,&ctx->y);
95: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->x);
96: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->y);
97: }
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: if (eps->which!=EPS_ALL || eps->inta==eps->intb) SETERRQ(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: MatCreateVecs(ctx->A,&v0,NULL);
107: if (!eps->V) { EPSGetBV(eps,&eps->V); }
108: BVGetRandomContext(eps->V,&rnd);
109: VecSetRandom(v0,rnd);
110: VecGetArray(v0,&vinit);
111: LanTrbounds(50,200,eps->tol,vinit,1,&ctx->lmin,&ctx->lmax,NULL);
112: VecRestoreArray(v0,&vinit);
113: VecDestroy(&v0);
114: ctx->estimrange = PETSC_TRUE; /* estimate if called again with another matrix */
115: }
116: if (ctx->lmin > eps->inta || ctx->lmax < eps->intb) SETERRQ4(PetscObjectComm((PetscObject)eps),1,"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: if (ctx->dos == EPS_EVSL_DOS_KPM) {
124: PetscMalloc1(ctx->deg+1,&mu);
125: if (!rank) { kpmdos(ctx->deg,(int)ctx->damping,ctx->nvec,xintv,mu,&ecount); }
126: MPI_Bcast(mu,ctx->deg+1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
127: } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
128: PetscMalloc2(ctx->npoints,&xdos,ctx->npoints,&ydos);
129: if (!rank) { LanDos(ctx->nvec,PetscMin(ctx->steps,eps->n/2),ctx->npoints,xdos,ydos,&ecount,xintv); }
130: MPI_Bcast(xdos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
131: MPI_Bcast(ydos,ctx->npoints,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
132: } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid DOS method");
133: MPI_Bcast(&ecount,1,MPIU_REAL,0,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
135: PetscInfo1(eps,"Estimated eigenvalue count in the interval: %g\n",ecount);
136: eps->ncv = (PetscInt)PetscCeilReal(1.5*ecount);
138: /* slice the spectrum */
139: PetscFree(ctx->sli);
140: PetscMalloc1(ctx->nslices+1,&ctx->sli);
141: if (ctx->dos == EPS_EVSL_DOS_KPM) {
142: spslicer(ctx->sli,mu,ctx->deg,xintv,ctx->nslices,10*(PetscInt)ecount);
143: PetscFree(mu);
144: } else if (ctx->dos == EPS_EVSL_DOS_LANCZOS) {
145: spslicer2(xdos,ydos,ctx->nslices,ctx->npoints,ctx->sli);
146: PetscFree2(xdos,ydos);
147: }
149: /* approximate number of eigenvalues wanted in each slice */
150: ctx->nev = (PetscInt)(1.0 + ecount/(PetscReal)ctx->nslices) + 2;
152: if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
153: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
154: EPSAllocateSolution(eps,0);
155: return(0);
156: }
158: PetscErrorCode EPSSolve_EVSL(EPS eps)
159: {
161: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
162: PetscInt i,j,k=0,sl,mlan,nevout,*ind,nevmax,rstart,rend,*nevloc,*disp,N;
163: PetscReal *res,xintv[4],*errest;
164: PetscScalar *lam,*X,*Y,*vinit,*eigr;
165: PetscMPIInt size,rank;
166: PetscRandom rnd;
167: Vec v,w,v0,x;
168: VecScatter vs;
169: IS is;
170: polparams pol;
173: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);CHKERRMPI(ierr);
174: MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank);CHKERRMPI(ierr);
175: PetscLayoutGetRange(ctx->map,&rstart,&rend);
176: nevmax = (rend-rstart)*ctx->nev;
177: MatCreateVecs(ctx->A,&v0,NULL);
178: BVGetRandomContext(eps->V,&rnd);
179: VecSetRandom(v0,rnd);
180: VecGetArray(v0,&vinit);
181: PetscMalloc5(size,&nevloc,size+1,&disp,nevmax,&eigr,nevmax,&errest,nevmax*eps->n,&X);
182: mlan = PetscMin(PetscMax(5*ctx->nev,300),eps->n);
183: for (sl=rstart; sl<rend; sl++) {
184: xintv[0] = ctx->sli[sl];
185: xintv[1] = ctx->sli[sl+1];
186: xintv[2] = ctx->lmin;
187: xintv[3] = ctx->lmax;
188: PetscInfo3(ctx->A,"Subinterval %D: [%.4e, %.4e]\n",sl+1,xintv[0],xintv[1]);
189: set_pol_def(&pol);
190: pol.max_deg = ctx->max_deg;
191: pol.damping = (int)ctx->damping;
192: pol.thresh_int = ctx->thresh;
193: find_pol(xintv,&pol);
194: PetscInfo4(ctx->A,"Polynomial [type = %D], deg %D, bar %e gam %e\n",pol.type,pol.deg,pol.bar,pol.gam);
195: ChebLanNr(xintv,mlan,eps->tol,vinit,&pol,&nevout,&lam,&Y,&res,NULL);
196: if (k+nevout>nevmax) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
197: free_pol(&pol);
198: PetscInfo1(ctx->A,"Computed %D eigenvalues\n",nevout);
199: PetscMalloc1(nevout,&ind);
200: sort_double(nevout,lam,ind);
201: for (i=0;i<nevout;i++) {
202: eigr[i+k] = lam[i];
203: errest[i+k] = res[ind[i]];
204: PetscArraycpy(X+(i+k)*eps->n,Y+ind[i]*eps->n,eps->n);
205: }
206: k += nevout;
207: if (lam) evsl_Free(lam);
208: if (Y) evsl_Free_device(Y);
209: if (res) evsl_Free(res);
210: PetscFree(ind);
211: }
212: VecRestoreArray(v0,&vinit);
213: VecDestroy(&v0);
215: /* gather eigenvalues computed by each MPI process */
216: MPI_Allgather(&k,1,MPIU_INT,nevloc,1,MPIU_INT,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
217: eps->nev = nevloc[0];
218: disp[0] = 0;
219: for (i=1;i<size;i++) {
220: eps->nev += nevloc[i];
221: disp[i] = disp[i-1]+nevloc[i-1];
222: }
223: disp[size] = disp[size-1]+nevloc[size-1];
224: if (eps->nev>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Too low estimation of eigenvalue count, try modifying the sampling parameters");
225: MPI_Allgatherv(eigr,k,MPIU_SCALAR,eps->eigr,nevloc,disp,MPIU_SCALAR,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
226: MPI_Allgatherv(errest,k,MPIU_REAL,eps->errest,nevloc,disp,MPIU_REAL,PetscObjectComm((PetscObject)eps));CHKERRMPI(ierr);
227: eps->nconv = eps->nev;
228: eps->its = 1;
229: eps->reason = EPS_CONVERGED_TOL;
231: /* scatter computed eigenvectors and store them in eps->V */
232: BVCreateVec(eps->V,&w);
233: for (i=0;i<size;i++) {
234: N = (rank==i)? eps->n: 0;
235: VecCreateSeq(PETSC_COMM_SELF,N,&x);
236: VecSetFromOptions(x);
237: ISCreateStride(PETSC_COMM_SELF,N,0,1,&is);
238: VecScatterCreate(x,is,w,is,&vs);
239: ISDestroy(&is);
240: for (j=disp[i];j<disp[i+1];j++) {
241: BVGetColumn(eps->V,j,&v);
242: if (rank==i) { VecPlaceArray(x,X+(j-disp[i])*eps->n); }
243: VecScatterBegin(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
244: VecScatterEnd(vs,x,v,INSERT_VALUES,SCATTER_FORWARD);
245: if (rank==i) { VecResetArray(x); }
246: BVRestoreColumn(eps->V,j,&v);
247: }
248: VecScatterDestroy(&vs);
249: VecDestroy(&x);
250: }
251: VecDestroy(&w);
252: PetscFree5(nevloc,disp,eigr,errest,X);
253: return(0);
254: }
256: static PetscErrorCode EPSEVSLSetSlices_EVSL(EPS eps,PetscInt nslices)
257: {
258: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
261: if (nslices == PETSC_DECIDE || nslices == PETSC_DEFAULT) nslices = 0;
262: else if (nslices<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Number of slices must be 1 at least");
263: if (ctx->nslices != nslices) {
264: ctx->nslices = nslices;
265: eps->state = EPS_STATE_INITIAL;
266: }
267: return(0);
268: }
270: /*@
271: EPSEVSLSetSlices - Set the number of slices in which the interval must be
272: subdivided.
274: Logically Collective on eps
276: Input Parameters:
277: + eps - the eigensolver context
278: - nslices - the number of slices
280: Options Database Key:
281: . -eps_evsl_slices <n> - set the number of slices to n
283: Notes:
284: By default, one slice per MPI process is used. Depending on the number of
285: eigenvalues, using more slices may be beneficial, but very narrow subintervals
286: imply higher polynomial degree.
288: Level: intermediate
290: .seealso: EPSEVSLGetSlices()
291: @*/
292: PetscErrorCode EPSEVSLSetSlices(EPS eps,PetscInt nslices)
293: {
299: PetscTryMethod(eps,"EPSEVSLSetSlices_C",(EPS,PetscInt),(eps,nslices));
300: return(0);
301: }
303: static PetscErrorCode EPSEVSLGetSlices_EVSL(EPS eps,PetscInt *nslices)
304: {
305: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
308: *nslices = ctx->nslices;
309: return(0);
310: }
312: /*@
313: EPSEVSLGetSlices - Gets the number of slices in which the interval must be
314: subdivided.
316: Not Collective
318: Input Parameter:
319: . eps - the eigensolver context
321: Output Parameter:
322: . nslices - the number of slices
324: Level: intermediate
326: .seealso: EPSEVSLSetSlices()
327: @*/
328: PetscErrorCode EPSEVSLGetSlices(EPS eps,PetscInt *nslices)
329: {
335: PetscUseMethod(eps,"EPSEVSLGetSlices_C",(EPS,PetscInt*),(eps,nslices));
336: return(0);
337: }
339: static PetscErrorCode EPSEVSLSetRange_EVSL(EPS eps,PetscReal lmin,PetscReal lmax)
340: {
341: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
344: if (lmin>lmax) SETERRQ(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: return(0);
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 on eps
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: {
387: PetscTryMethod(eps,"EPSEVSLSetRange_C",(EPS,PetscReal,PetscReal),(eps,lmin,lmax));
388: return(0);
389: }
391: static PetscErrorCode EPSEVSLGetRange_EVSL(EPS eps,PetscReal *lmin,PetscReal *lmax)
392: {
393: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
396: if (lmin) *lmin = ctx->lmin;
397: if (lmax) *lmax = ctx->lmax;
398: return(0);
399: }
401: /*@
402: EPSEVSLGetRange - Gets the interval containing all eigenvalues.
404: Not Collective
406: Input Parameter:
407: . eps - the eigensolver context
409: Output Parameters:
410: + lmin - left end of the interval
411: - lmax - right end of the interval
413: Level: intermediate
415: .seealso: EPSEVSLSetRange()
416: @*/
417: PetscErrorCode EPSEVSLGetRange(EPS eps,PetscReal *lmin,PetscReal *lmax)
418: {
423: PetscUseMethod(eps,"EPSEVSLGetRange_C",(EPS,PetscReal*,PetscReal*),(eps,lmin,lmax));
424: return(0);
425: }
427: static PetscErrorCode EPSEVSLSetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
428: {
429: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
432: ctx->dos = dos;
433: if (nvec == PETSC_DECIDE || nvec == PETSC_DEFAULT) ctx->nvec = 80;
434: else if (nvec<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The nvec argument must be > 0");
435: else ctx->nvec = nvec;
436: switch (dos) {
437: case EPS_EVSL_DOS_KPM:
438: if (deg == PETSC_DECIDE || deg == PETSC_DEFAULT) ctx->deg = 300;
439: else if (deg<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The deg argument must be > 0");
440: else ctx->deg = deg;
441: break;
442: case EPS_EVSL_DOS_LANCZOS:
443: if (steps == PETSC_DECIDE || steps == PETSC_DEFAULT) ctx->steps = 40;
444: else if (steps<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The steps argument must be > 0");
445: else ctx->steps = steps;
446: if (npoints == PETSC_DECIDE || npoints == PETSC_DEFAULT) ctx->npoints = 200;
447: else if (npoints<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The npoints argument must be > 0");
448: else ctx->npoints = npoints;
449: break;
450: }
451: eps->state = EPS_STATE_INITIAL;
452: return(0);
453: }
455: /*@
456: EPSEVSLSetDOSParameters - Defines the parameters used for computing the
457: density of states (DOS) in the EVSL solver.
459: Logically Collective on eps
461: Input Parameters:
462: + eps - the eigensolver context
463: . dos - DOS method, either KPM or Lanczos
464: . nvec - number of sample vectors
465: . deg - polynomial degree (KPM only)
466: . steps - number of Lanczos steps (Lanczos only)
467: - npoints - number of sample points (Lanczos only)
469: Options Database Keys:
470: + -eps_evsl_dos_method <dos> - set the DOS method, either kpm or lanczos
471: . -eps_evsl_dos_nvec <n> - set the number of sample vectors
472: . -eps_evsl_dos_degree <n> - set the polynomial degree
473: . -eps_evsl_dos_steps <n> - set the number of Lanczos steps
474: - -eps_evsl_dos_npoints <n> - set the number of sample points
476: Notes:
477: The density of states (or spectral density) can be approximated with two
478: methods: kernel polynomial method (KPM) or Lanczos. Some parameters for
479: these methods can be set by the user with this function, with some of
480: them being relevant for one of the methods only.
482: Level: intermediate
484: .seealso: EPSEVSLGetDOSParameters()
485: @*/
486: PetscErrorCode EPSEVSLSetDOSParameters(EPS eps,EPSEVSLDOSMethod dos,PetscInt nvec,PetscInt deg,PetscInt steps,PetscInt npoints)
487: {
497: PetscTryMethod(eps,"EPSEVSLSetDOSParameters_C",(EPS,EPSEVSLDOSMethod,PetscInt,PetscInt,PetscInt,PetscInt),(eps,dos,nvec,deg,steps,npoints));
498: return(0);
499: }
501: static PetscErrorCode EPSEVSLGetDOSParameters_EVSL(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
502: {
503: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
506: if (dos) *dos = ctx->dos;
507: if (nvec) *nvec = ctx->nvec;
508: if (deg) *deg = ctx->deg;
509: if (steps) *steps = ctx->steps;
510: if (npoints) *npoints = ctx->npoints;
511: return(0);
512: }
514: /*@
515: EPSEVSLGetDOSParameters - Gets the parameters used for computing the
516: density of states (DOS) in the EVSL solver.
518: Not Collective
520: Input Parameter:
521: . eps - the eigensolver context
523: Output Parameters:
524: + dos - DOS method, either KPM or Lanczos
525: . nvec - number of sample vectors
526: . deg - polynomial degree (KPM only)
527: . steps - number of Lanczos steps (Lanczos only)
528: - npoints - number of sample points (Lanczos only)
530: Level: intermediate
532: .seealso: EPSEVSLSetDOSParameters()
533: @*/
534: PetscErrorCode EPSEVSLGetDOSParameters(EPS eps,EPSEVSLDOSMethod *dos,PetscInt *nvec,PetscInt *deg,PetscInt *steps,PetscInt *npoints)
535: {
540: PetscUseMethod(eps,"EPSEVSLGetDOSParameters_C",(EPS,EPSEVSLDOSMethod*,PetscInt*,PetscInt*,PetscInt*,PetscInt*),(eps,dos,nvec,deg,steps,npoints));
541: return(0);
542: }
544: static PetscErrorCode EPSEVSLSetPolParameters_EVSL(EPS eps,PetscInt max_deg,PetscReal thresh)
545: {
546: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
549: if (max_deg == PETSC_DECIDE || max_deg == PETSC_DEFAULT) ctx->max_deg = 10000;
550: else if (max_deg<3) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The max_deg argument must be > 2");
551: else ctx->max_deg = max_deg;
552: if (thresh == PETSC_DECIDE || thresh == PETSC_DEFAULT) ctx->thresh = 0.8;
553: else if (thresh<0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The thresh argument must be > 0.0");
554: else ctx->thresh = thresh;
555: eps->state = EPS_STATE_INITIAL;
556: return(0);
557: }
559: /*@
560: EPSEVSLSetPolParameters - Defines the parameters used for building the
561: building the polynomial in the EVSL solver.
563: Logically Collective on eps
565: Input Parameters:
566: + eps - the eigensolver context
567: . max_deg - maximum degree allowed for the polynomial
568: - thresh - threshold for accepting polynomial
570: Options Database Keys:
571: + -eps_evsl_pol_max_deg <d> - set maximum polynomial degree
572: - -eps_evsl_pol_thresh <t> - set the threshold
574: Level: intermediate
576: .seealso: EPSEVSLGetPolParameters()
577: @*/
578: PetscErrorCode EPSEVSLSetPolParameters(EPS eps,PetscInt max_deg,PetscReal thresh)
579: {
586: PetscTryMethod(eps,"EPSEVSLSetPolParameters_C",(EPS,PetscInt,PetscReal),(eps,max_deg,thresh));
587: return(0);
588: }
590: static PetscErrorCode EPSEVSLGetPolParameters_EVSL(EPS eps,PetscInt *max_deg,PetscReal *thresh)
591: {
592: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
595: if (max_deg) *max_deg = ctx->max_deg;
596: if (thresh) *thresh = ctx->thresh;
597: return(0);
598: }
600: /*@
601: EPSEVSLGetPolParameters - Gets the parameters used for building the
602: polynomial in the EVSL solver.
604: Not Collective
606: Input Parameter:
607: . eps - the eigensolver context
609: Output Parameters:
610: + max_deg - the maximum degree of the polynomial
611: - thresh - the threshold
613: Level: intermediate
615: .seealso: EPSEVSLSetPolParameters()
616: @*/
617: PetscErrorCode EPSEVSLGetPolParameters(EPS eps,PetscInt *max_deg,PetscReal *thresh)
618: {
623: PetscUseMethod(eps,"EPSEVSLGetPolParameters_C",(EPS,PetscInt*,PetscReal*),(eps,max_deg,thresh));
624: return(0);
625: }
627: static PetscErrorCode EPSEVSLSetDamping_EVSL(EPS eps,EPSEVSLDamping damping)
628: {
629: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
632: if (ctx->damping != damping) {
633: ctx->damping = damping;
634: eps->state = EPS_STATE_INITIAL;
635: }
636: return(0);
637: }
639: /*@
640: EPSEVSLSetDamping - Set the type of damping to be used in EVSL.
642: Logically Collective on eps
644: Input Parameters:
645: + eps - the eigensolver context
646: - damping - the type of damping
648: Options Database Key:
649: . -eps_evsl_damping <n> - set the type of damping
651: Notes:
652: Damping is applied when building the polynomial to be used when solving the
653: eigenproblem, and also during estimation of DOS with the KPM method.
655: Level: intermediate
657: .seealso: EPSEVSLGetDamping(), EPSEVSLSetDOSParameters()
658: @*/
659: PetscErrorCode EPSEVSLSetDamping(EPS eps,EPSEVSLDamping damping)
660: {
666: PetscTryMethod(eps,"EPSEVSLSetDamping_C",(EPS,EPSEVSLDamping),(eps,damping));
667: return(0);
668: }
670: static PetscErrorCode EPSEVSLGetDamping_EVSL(EPS eps,EPSEVSLDamping *damping)
671: {
672: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
675: *damping = ctx->damping;
676: return(0);
677: }
679: /*@
680: EPSEVSLGetDamping - Gets the type of damping.
682: Not Collective
684: Input Parameter:
685: . eps - the eigensolver context
687: Output Parameter:
688: . damping - the type of damping
690: Level: intermediate
692: .seealso: EPSEVSLSetDamping()
693: @*/
694: PetscErrorCode EPSEVSLGetDamping(EPS eps,EPSEVSLDamping *damping)
695: {
701: PetscUseMethod(eps,"EPSEVSLGetDamping_C",(EPS,EPSEVSLDamping*),(eps,damping));
702: return(0);
703: }
705: PetscErrorCode EPSView_EVSL(EPS eps,PetscViewer viewer)
706: {
708: PetscBool isascii;
709: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
712: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
713: if (isascii) {
714: PetscViewerASCIIPrintf(viewer," numerical range = [%g,%g]\n",(double)ctx->lmin,(double)ctx->lmax);
715: PetscViewerASCIIPrintf(viewer," number of slices = %D\n",ctx->nslices);
716: PetscViewerASCIIPrintf(viewer," type of damping = %s\n",EPSEVSLDampings[ctx->damping]);
717: PetscViewerASCIIPrintf(viewer," computing DOS with %s: nvec=%D, ",EPSEVSLDOSMethods[ctx->dos],ctx->nvec);
718: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
719: switch (ctx->dos) {
720: case EPS_EVSL_DOS_KPM:
721: PetscViewerASCIIPrintf(viewer,"degree=%D\n",ctx->deg);
722: break;
723: case EPS_EVSL_DOS_LANCZOS:
724: PetscViewerASCIIPrintf(viewer,"steps=%D, npoints=%D\n",ctx->steps,ctx->npoints);
725: break;
726: }
727: PetscViewerASCIIPrintf(viewer," polynomial parameters: max degree = %D, threshold = %g\n",ctx->max_deg,(double)ctx->thresh);
728: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
729: }
730: return(0);
731: }
733: PetscErrorCode EPSSetFromOptions_EVSL(PetscOptionItems *PetscOptionsObject,EPS eps)
734: {
735: PetscErrorCode ierr;
736: PetscReal array[2]={0,0},th;
737: PetscInt k,i1,i2,i3,i4;
738: PetscBool flg,flg1;
739: EPSEVSLDOSMethod dos;
740: EPSEVSLDamping damping;
741: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
744: PetscOptionsHead(PetscOptionsObject,"EPS EVSL Options");
746: k = 2;
747: PetscOptionsRealArray("-eps_evsl_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","EPSEVSLSetRange",array,&k,&flg);
748: if (flg) {
749: if (k<2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_SIZ,"Must pass two values in -eps_evsl_range (comma-separated without spaces)");
750: EPSEVSLSetRange(eps,array[0],array[1]);
751: }
753: PetscOptionsInt("-eps_evsl_slices","Number of slices","EPSEVSLSetSlices",ctx->nslices,&i1,&flg);
754: if (flg) { EPSEVSLSetSlices(eps,i1); }
756: PetscOptionsEnum("-eps_evsl_damping","Type of damping","EPSEVSLSetDamping",EPSEVSLDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg);
757: if (flg) { EPSEVSLSetDamping(eps,damping); }
759: EPSEVSLGetDOSParameters(eps,&dos,&i1,&i2,&i3,&i4);
760: PetscOptionsEnum("-eps_evsl_dos_method","Method to compute the DOS","EPSEVSLSetDOSParameters",EPSEVSLDOSMethods,(PetscEnum)ctx->dos,(PetscEnum*)&dos,&flg);
761: PetscOptionsInt("-eps_evsl_dos_nvec","Number of sample vectors for DOS","EPSEVSLSetDOSParameters",i1,&i1,&flg1);
762: flg = flg || flg1;
763: PetscOptionsInt("-eps_evsl_dos_degree","Polynomial degree used for DOS","EPSEVSLSetDOSParameters",i2,&i2,&flg1);
764: flg = flg || flg1;
765: PetscOptionsInt("-eps_evsl_dos_steps","Number of Lanczos steps in DOS","EPSEVSLSetDOSParameters",i3,&i3,&flg1);
766: flg = flg || flg1;
767: PetscOptionsInt("-eps_evsl_dos_npoints","Number of sample points used for DOS","EPSEVSLSetDOSParameters",i4,&i4,&flg1);
768: if (flg || flg1) { EPSEVSLSetDOSParameters(eps,dos,i1,i2,i3,i4); }
770: EPSEVSLGetPolParameters(eps,&i1,&th);
771: PetscOptionsInt("-eps_evsl_pol_max_deg","Maximum degree allowed for the polynomial","EPSEVSLSetPolParameters",i1,&i1,&flg);
772: PetscOptionsReal("-eps_evsl_pol_threshold","Threshold for accepting polynomial","EPSEVSLSetPolParameters",th,&th,&flg1);
773: if (flg || flg1) { EPSEVSLSetPolParameters(eps,i1,th); }
775: PetscOptionsTail();
776: return(0);
777: }
779: PetscErrorCode EPSDestroy_EVSL(EPS eps)
780: {
782: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
785: if (ctx->initialized) EVSLFinish();
786: PetscLayoutDestroy(&ctx->map);
787: PetscFree(ctx->sli);
788: PetscFree(eps->data);
789: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",NULL);
790: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",NULL);
791: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",NULL);
792: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",NULL);
793: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",NULL);
794: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",NULL);
795: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",NULL);
796: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",NULL);
797: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",NULL);
798: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",NULL);
799: return(0);
800: }
802: PetscErrorCode EPSReset_EVSL(EPS eps)
803: {
805: EPS_EVSL *ctx = (EPS_EVSL*)eps->data;
808: MatDestroy(&ctx->A);
809: VecDestroy(&ctx->x);
810: VecDestroy(&ctx->y);
811: return(0);
812: }
814: SLEPC_EXTERN PetscErrorCode EPSCreate_EVSL(EPS eps)
815: {
816: EPS_EVSL *ctx;
820: PetscNewLog(eps,&ctx);
821: eps->data = (void*)ctx;
823: ctx->nslices = 0;
824: ctx->lmin = PETSC_MIN_REAL;
825: ctx->lmax = PETSC_MAX_REAL;
826: ctx->dos = EPS_EVSL_DOS_KPM;
827: ctx->nvec = 80;
828: ctx->deg = 300;
829: ctx->steps = 40;
830: ctx->npoints = 200;
831: ctx->max_deg = 10000;
832: ctx->thresh = 0.8;
833: ctx->damping = EPS_EVSL_DAMPING_SIGMA;
835: eps->categ = EPS_CATEGORY_OTHER;
837: eps->ops->solve = EPSSolve_EVSL;
838: eps->ops->setup = EPSSetUp_EVSL;
839: eps->ops->setupsort = EPSSetUpSort_Basic;
840: eps->ops->setfromoptions = EPSSetFromOptions_EVSL;
841: eps->ops->destroy = EPSDestroy_EVSL;
842: eps->ops->reset = EPSReset_EVSL;
843: eps->ops->view = EPSView_EVSL;
844: eps->ops->backtransform = EPSBackTransform_Default;
845: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
847: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetRange_C",EPSEVSLSetRange_EVSL);
848: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetRange_C",EPSEVSLGetRange_EVSL);
849: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetSlices_C",EPSEVSLSetSlices_EVSL);
850: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetSlices_C",EPSEVSLGetSlices_EVSL);
851: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDOSParameters_C",EPSEVSLSetDOSParameters_EVSL);
852: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDOSParameters_C",EPSEVSLGetDOSParameters_EVSL);
853: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetPolParameters_C",EPSEVSLSetPolParameters_EVSL);
854: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetPolParameters_C",EPSEVSLGetPolParameters_EVSL);
855: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLSetDamping_C",EPSEVSLSetDamping_EVSL);
856: PetscObjectComposeFunction((PetscObject)eps,"EPSEVSLGetDamping_C",EPSEVSLGetDamping_EVSL);
857: return(0);
858: }