Actual source code: chase.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 ChASE.
12: */
14: #include <slepc/private/epsimpl.h>
15: #include <petsc/private/petscscalapack.h>
17: #if !defined(PETSC_USE_COMPLEX)
18: #if defined(PETSC_USE_REAL_SINGLE)
19: /* s */
20: #define CHASEchase_init_blockcyclic pschase_init_blockcyclic_
21: #define CHASEchase pschase_
22: #define CHASEchase_finalize pschase_finalize_
23: #elif defined(PETSC_USE_REAL_DOUBLE)
24: /* d */
25: #define CHASEchase_init_blockcyclic pdchase_init_blockcyclic_
26: #define CHASEchase pdchase_
27: #define CHASEchase_finalize pdchase_finalize_
28: #endif
29: #else
30: #if defined(PETSC_USE_REAL_SINGLE)
31: /* c */
32: #define CHASEchase_init_blockcyclic pcchase_init_blockcyclic_
33: #define CHASEchase pcchase_
34: #define CHASEchase_finalize pcchase_finalize_
35: #elif defined(PETSC_USE_REAL_DOUBLE)
36: /* z */
37: #define CHASEchase_init_blockcyclic pzchase_init_blockcyclic_
38: #define CHASEchase pzchase_
39: #define CHASEchase_finalize pzchase_finalize_
40: #endif
41: #endif
43: SLEPC_EXTERN void CHASEchase_init_blockcyclic(PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscScalar*,PetscInt*,PetscScalar*,PetscReal*,PetscInt*,PetscInt*,char*,PetscInt*,PetscInt*,MPI_Comm*,PetscInt*);
44: SLEPC_EXTERN void CHASEchase(PetscInt*,PetscReal*,char*,char*,char*);
45: SLEPC_EXTERN void CHASEchase_finalize(PetscInt*);
47: typedef struct {
48: Mat As; /* converted matrix */
49: PetscInt nex; /* extra searching space size */
50: PetscInt deg; /* initial degree of Chebyshev polynomial filter */
51: PetscBool opt; /* internal optimization of polynomial degree */
52: } EPS_ChASE;
54: static PetscErrorCode EPSSetUp_ChASE(EPS eps)
55: {
56: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
57: Mat A;
58: PetscBool isshift;
59: PetscScalar shift;
61: PetscFunctionBegin;
62: EPSCheckStandard(eps);
63: EPSCheckHermitian(eps);
64: EPSCheckNotStructured(eps);
65: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
66: PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
67: if (eps->nev==0) eps->nev = 1;
68: if (eps->ncv==PETSC_DETERMINE) eps->ncv = PetscMin(eps->n,PetscMax(2*eps->nev,eps->nev+15));
69: else PetscCheck(eps->ncv>=eps->nev+1 || (eps->ncv==eps->nev && eps->ncv==eps->n),PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
70: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
71: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
72: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
73: PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports computation of leftmost eigenvalues");
74: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION);
75: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
76: PetscCall(EPSAllocateSolution(eps,0));
78: ctx->nex = eps->ncv-eps->nev;
79: if (ctx->deg==PETSC_DEFAULT) ctx->deg = 10;
81: /* convert matrices */
82: PetscCall(MatDestroy(&ctx->As));
83: PetscCall(STGetMatrix(eps->st,0,&A));
84: PetscCall(MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
85: PetscCall(STGetShift(eps->st,&shift));
86: if (shift != 0.0) PetscCall(MatShift(ctx->As,-shift));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode EPSSolve_ChASE(EPS eps)
91: {
92: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
93: Mat As = ctx->As,Vs,V;
94: Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)As->data,*v;
95: MPI_Comm comm;
96: PetscReal *w = eps->errest; /* used to store real eigenvalues */
97: PetscInt i,init=0,flg=1;
98: char grid_major='R',mode='X',opt=ctx->opt?'S':'X',qr='C';
100: PetscFunctionBegin;
101: PetscCheck(a->grid->npcol==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Current implementation of the ChASE interface only supports process grids with one column; use -mat_scalapack_grid_height <p> where <p> is the number of processes");
102: PetscCall(BVGetMat(eps->V,&V));
103: PetscCall(MatConvert(V,MATSCALAPACK,MAT_INITIAL_MATRIX,&Vs));
104: PetscCall(BVRestoreMat(eps->V,&V));
105: v = (Mat_ScaLAPACK*)Vs->data;
107: /* initialization */
108: PetscCall(PetscObjectGetComm((PetscObject)As,&comm));
109: CHASEchase_init_blockcyclic(&a->N,&eps->nev,&ctx->nex,&a->mb,&a->nb,a->loc,&a->lld,v->loc,w,&a->grid->nprow,&a->grid->npcol,&grid_major,&a->rsrc,&a->csrc,&comm,&init);
110: PetscCheck(init==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Problem initializing ChASE");
112: /* solve */
113: CHASEchase(&ctx->deg,&eps->tol,&mode,&opt,&qr);
114: CHASEchase_finalize(&flg);
115: PetscCheck(flg==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Problem solving with ChASE");
117: for (i=0;i<eps->ncv;i++) {
118: eps->eigr[i] = eps->errest[i];
119: eps->errest[i] = eps->tol;
120: }
122: PetscCall(BVGetMat(eps->V,&V));
123: PetscCall(MatConvert(Vs,MATDENSE,MAT_REUSE_MATRIX,&V));
124: PetscCall(BVRestoreMat(eps->V,&V));
125: PetscCall(MatDestroy(&Vs));
127: eps->nconv = eps->nev;
128: eps->its = 1;
129: eps->reason = EPS_CONVERGED_TOL;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode EPSCHASESetDegree_ChASE(EPS eps,PetscInt deg,PetscBool opt)
134: {
135: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
137: PetscFunctionBegin;
138: if (deg==PETSC_DEFAULT || deg==PETSC_DECIDE) {
139: ctx->deg = PETSC_DEFAULT;
140: eps->state = EPS_STATE_INITIAL;
141: } else {
142: PetscCheck(deg>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"The degree must be >0");
143: ctx->deg = deg;
144: }
145: ctx->opt = opt;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: EPSCHASESetDegree - Sets the degree of the Chebyshev polynomial filter in the ChASE solver.
152: Logically Collective
154: Input Parameters:
155: + eps - the linear eigensolver context
156: . deg - initial degree of Chebyshev polynomial filter
157: - opt - internal optimization of polynomial degree
159: Options Database Keys:
160: + -eps_chase_degree \<deg\> - set the initial degree
161: - -eps_chase_degree_opt - toggle the optimization
163: Note:
164: See the documentation of ChASE {cite:p}`Win19` for details.
166: Level: advanced
168: .seealso: [](ch:eps), `EPSCHASE`, `EPSCHASEGetDegree()`
169: @*/
170: PetscErrorCode EPSCHASESetDegree(EPS eps,PetscInt deg,PetscBool opt)
171: {
172: PetscFunctionBegin;
176: PetscTryMethod(eps,"EPSCHASESetDegree_C",(EPS,PetscInt,PetscBool),(eps,deg,opt));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode EPSCHASEGetDegree_ChASE(EPS eps,PetscInt *deg,PetscBool *opt)
181: {
182: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
184: PetscFunctionBegin;
185: if (deg) *deg = ctx->deg;
186: if (opt) *opt = ctx->opt;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*@
191: EPSCHASEGetDegree - Gets the degree of the Chebyshev polynomial filter used in the ChASE solver.
193: Not Collective
195: Input Parameter:
196: . eps - the linear eigensolver context
198: Output Parameters:
199: + deg - initial degree of Chebyshev polynomial filter
200: - opt - internal optimization of polynomial degree
202: Level: advanced
204: .seealso: [](ch:eps), `EPSCHASE`, `EPSCHASESetDegree()`
205: @*/
206: PetscErrorCode EPSCHASEGetDegree(EPS eps,PetscInt *deg,PetscBool *opt)
207: {
208: PetscFunctionBegin;
210: PetscUseMethod(eps,"EPSCHASEGetDegree_C",(EPS,PetscInt*,PetscBool*),(eps,deg,opt));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode EPSSetFromOptions_ChASE(EPS eps,PetscOptionItems PetscOptionsObject)
215: {
216: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
217: PetscBool flg1,flg2;
218: PetscInt deg;
219: PetscBool opt;
221: PetscFunctionBegin;
222: PetscOptionsHeadBegin(PetscOptionsObject,"EPS ChASE Options");
224: deg = ctx->deg;
225: PetscCall(PetscOptionsInt("-eps_chase_degree","Initial degree of Chebyshev polynomial filter","EPSCHASESetDegree",deg,°,&flg1));
226: opt = ctx->opt;
227: PetscCall(PetscOptionsBool("-eps_chase_degree_opt","Internal optimization of polynomial degree","EPSCHASESetDegree",opt,&opt,&flg2));
228: if (flg1 || flg2) PetscCall(EPSCHASESetDegree(eps,deg,opt));
230: PetscOptionsHeadEnd();
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: static PetscErrorCode EPSDestroy_ChASE(EPS eps)
235: {
236: PetscFunctionBegin;
237: PetscCall(PetscFree(eps->data));
238: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",NULL));
239: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",NULL));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode EPSReset_ChASE(EPS eps)
244: {
245: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
247: PetscFunctionBegin;
248: PetscCall(MatDestroy(&ctx->As));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode EPSView_ChASE(EPS eps,PetscViewer viewer)
253: {
254: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
255: PetscBool isascii;
257: PetscFunctionBegin;
258: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
259: if (isascii) {
260: PetscCall(PetscViewerASCIIPrintf(viewer," initial degree of Chebyshev polynomial filter: %" PetscInt_FMT "\n",ctx->deg));
261: if (ctx->opt) PetscCall(PetscViewerASCIIPrintf(viewer," internal optimization of polynomial degree enabled\n"));
262: }
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*MC
267: EPSCHASE - EPSCHASE = "chase" - A wrapper to CHASE {cite:p}`Win19`.
269: Notes:
270: The Chebyshev Accelerated Subspace Iteration Eigensolver (ChASE) is
271: expected to be particularly efficient when computing many eigenvalues
272: at the lower end of the spectrum, in sequences of eigenproblems.
274: The current implementation of the interface in SLEPc is restricted
275: to the case where the matrix has a block cyclic distribution, so
276: SLEPc will redistribute the matrix to ScaLAPACK format.
278: Level: beginner
280: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
281: M*/
282: SLEPC_EXTERN PetscErrorCode EPSCreate_ChASE(EPS eps)
283: {
284: EPS_ChASE *ctx;
286: PetscFunctionBegin;
287: PetscCall(PetscNew(&ctx));
288: eps->data = (void*)ctx;
289: ctx->deg = PETSC_DEFAULT;
290: ctx->opt = PETSC_TRUE;
292: eps->categ = EPS_CATEGORY_OTHER;
294: eps->ops->solve = EPSSolve_ChASE;
295: eps->ops->setup = EPSSetUp_ChASE;
296: eps->ops->setupsort = EPSSetUpSort_Basic;
297: eps->ops->setfromoptions = EPSSetFromOptions_ChASE;
298: eps->ops->destroy = EPSDestroy_ChASE;
299: eps->ops->reset = EPSReset_ChASE;
300: eps->ops->view = EPSView_ChASE;
301: eps->ops->backtransform = EPSBackTransform_Default;
302: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
304: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",EPSCHASESetDegree_ChASE));
305: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",EPSCHASEGetDegree_ChASE));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }