Actual source code: chase.c
slepc-main 2024-12-17
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 eigenproblem solver context
156: . deg - initial degree of Chebyshev polynomial filter
157: - opt - internal optimization of polynomial degree
159: Options Database Keys:
160: + -eps_chase_degree - Sets the initial degree
161: - -eps_chase_degree_opt - Enables/disables the optimization
163: Level: advanced
165: .seealso: EPSCHASEGetDegree()
166: @*/
167: PetscErrorCode EPSCHASESetDegree(EPS eps,PetscInt deg,PetscBool opt)
168: {
169: PetscFunctionBegin;
173: PetscTryMethod(eps,"EPSCHASESetDegree_C",(EPS,PetscInt,PetscBool),(eps,deg,opt));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode EPSCHASEGetDegree_ChASE(EPS eps,PetscInt *deg,PetscBool *opt)
178: {
179: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
181: PetscFunctionBegin;
182: if (deg) *deg = ctx->deg;
183: if (opt) *opt = ctx->opt;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*@
188: EPSCHASEGetDegree - Gets the degree of the Chebyshev polynomial filter used in the ChASE solver.
190: Not Collective
192: Input Parameter:
193: . eps - the eigenproblem solver context
195: Output Parameters:
196: + deg - initial degree of Chebyshev polynomial filter
197: - opt - internal optimization of polynomial degree
199: Level: advanced
201: .seealso: EPSCHASESetDegree()
202: @*/
203: PetscErrorCode EPSCHASEGetDegree(EPS eps,PetscInt *deg,PetscBool *opt)
204: {
205: PetscFunctionBegin;
207: PetscUseMethod(eps,"EPSCHASEGetDegree_C",(EPS,PetscInt*,PetscBool*),(eps,deg,opt));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: static PetscErrorCode EPSSetFromOptions_ChASE(EPS eps,PetscOptionItems *PetscOptionsObject)
212: {
213: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
214: PetscBool flg1,flg2;
215: PetscInt deg;
216: PetscBool opt;
218: PetscFunctionBegin;
219: PetscOptionsHeadBegin(PetscOptionsObject,"EPS ChASE Options");
221: deg = ctx->deg;
222: PetscCall(PetscOptionsInt("-eps_chase_degree","Initial degree of Chebyshev polynomial filter","EPSCHASESetDegree",deg,°,&flg1));
223: opt = ctx->opt;
224: PetscCall(PetscOptionsBool("-eps_chase_degree_opt","Internal optimization of polynomial degree","EPSCHASESetDegree",opt,&opt,&flg2));
225: if (flg1 || flg2) PetscCall(EPSCHASESetDegree(eps,deg,opt));
227: PetscOptionsHeadEnd();
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode EPSDestroy_ChASE(EPS eps)
232: {
233: PetscFunctionBegin;
234: PetscCall(PetscFree(eps->data));
235: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",NULL));
236: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",NULL));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode EPSReset_ChASE(EPS eps)
241: {
242: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
244: PetscFunctionBegin;
245: PetscCall(MatDestroy(&ctx->As));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode EPSView_ChASE(EPS eps,PetscViewer viewer)
250: {
251: EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
252: PetscBool isascii;
254: PetscFunctionBegin;
255: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
256: if (isascii) {
257: PetscCall(PetscViewerASCIIPrintf(viewer," initial degree of Chebyshev polynomial filter: %" PetscInt_FMT "\n",ctx->deg));
258: if (ctx->opt) PetscCall(PetscViewerASCIIPrintf(viewer," internal optimization of polynomial degree enabled\n"));
259: }
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: SLEPC_EXTERN PetscErrorCode EPSCreate_ChASE(EPS eps)
264: {
265: EPS_ChASE *ctx;
267: PetscFunctionBegin;
268: PetscCall(PetscNew(&ctx));
269: eps->data = (void*)ctx;
270: ctx->deg = PETSC_DEFAULT;
271: ctx->opt = PETSC_TRUE;
273: eps->categ = EPS_CATEGORY_OTHER;
275: eps->ops->solve = EPSSolve_ChASE;
276: eps->ops->setup = EPSSetUp_ChASE;
277: eps->ops->setupsort = EPSSetUpSort_Basic;
278: eps->ops->setfromoptions = EPSSetFromOptions_ChASE;
279: eps->ops->destroy = EPSDestroy_ChASE;
280: eps->ops->reset = EPSReset_ChASE;
281: eps->ops->view = EPSView_ChASE;
282: eps->ops->backtransform = EPSBackTransform_Default;
283: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
285: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",EPSCHASESetDegree_ChASE));
286: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",EPSCHASEGetDegree_ChASE));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }