GCC Code Coverage Report

Directory: ./
File: src/eps/impls/external/chase/chase.c
Date: 2025-11-19 04:19:03
Exec Total Coverage
Lines: 125 136 91.9%
Functions: 10 11 90.9%
Branches: 63 186 33.9%
Line Branch Exec Source
1 /*
2 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 SLEPc - Scalable Library for Eigenvalue Problem Computations
4 Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5
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 */
13
14 #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 #include <petsc/private/petscscalapack.h>
16
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
42
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*);
46
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;
53
54 6 static PetscErrorCode EPSSetUp_ChASE(EPS eps)
55 {
56 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
57 6 Mat A;
58 6 PetscBool isshift;
59 6 PetscScalar shift;
60
61 6 PetscFunctionBegin;
62
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 EPSCheckStandard(eps);
63
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 EPSCheckHermitian(eps);
64
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 EPSCheckNotStructured(eps);
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
66
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
67
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 if (eps->nev==0) eps->nev = 1;
68
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 6 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
71
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 6 if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
72
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 if (!eps->which) eps->which = EPS_SMALLEST_REAL;
73
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports computation of leftmost eigenvalues");
74
4/14
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 6 EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION);
75
3/12
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 6 EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
76
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(EPSAllocateSolution(eps,0));
77
78 6 ctx->nex = eps->ncv-eps->nev;
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 if (ctx->deg==PETSC_DEFAULT) ctx->deg = 10;
80
81 /* convert matrices */
82
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatDestroy(&ctx->As));
83
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(STGetMatrix(eps->st,0,&A));
84
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(STGetShift(eps->st,&shift));
86
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 if (shift != 0.0) PetscCall(MatShift(ctx->As,-shift));
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
90 6 static PetscErrorCode EPSSolve_ChASE(EPS eps)
91 {
92 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
93 6 Mat As = ctx->As,Vs,V;
94 6 Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)As->data,*v;
95 6 MPI_Comm comm;
96 6 PetscReal *w = eps->errest; /* used to store real eigenvalues */
97 6 PetscInt i,init=0,flg=1;
98
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 char grid_major='R',mode='X',opt=ctx->opt?'S':'X',qr='C';
99
100 6 PetscFunctionBegin;
101
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(BVGetMat(eps->V,&V));
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatConvert(V,MATSCALAPACK,MAT_INITIAL_MATRIX,&Vs));
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(BVRestoreMat(eps->V,&V));
105 6 v = (Mat_ScaLAPACK*)Vs->data;
106
107 /* initialization */
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectGetComm((PetscObject)As,&comm));
109 6 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheck(init==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Problem initializing ChASE");
111
112 /* solve */
113 6 CHASEchase(&ctx->deg,&eps->tol,&mode,&opt,&qr);
114 6 CHASEchase_finalize(&flg);
115
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheck(flg==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Problem solving with ChASE");
116
117
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 132 for (i=0;i<eps->ncv;i++) {
118 126 eps->eigr[i] = eps->errest[i];
119 126 eps->errest[i] = eps->tol;
120 }
121
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(BVGetMat(eps->V,&V));
123
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatConvert(Vs,MATDENSE,MAT_REUSE_MATRIX,&V));
124
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(BVRestoreMat(eps->V,&V));
125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatDestroy(&Vs));
126
127 6 eps->nconv = eps->nev;
128 6 eps->its = 1;
129 6 eps->reason = EPS_CONVERGED_TOL;
130 6 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 6 static PetscErrorCode EPSCHASESetDegree_ChASE(EPS eps,PetscInt deg,PetscBool opt)
134 {
135 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
136
137 6 PetscFunctionBegin;
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 if (deg==PETSC_DEFAULT || deg==PETSC_DECIDE) {
139 ✗ ctx->deg = PETSC_DEFAULT;
140 ✗ eps->state = EPS_STATE_INITIAL;
141 } else {
142
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheck(deg>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"The degree must be >0");
143 6 ctx->deg = deg;
144 }
145 6 ctx->opt = opt;
146 6 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
149 /*@
150 EPSCHASESetDegree - Sets the degree of the Chebyshev polynomial filter in the ChASE solver.
151
152 Logically Collective
153
154 Input Parameters:
155 + eps - the linear eigensolver context
156 . deg - initial degree of Chebyshev polynomial filter
157 - opt - internal optimization of polynomial degree
158
159 Options Database Keys:
160 + -eps_chase_degree \<deg\> - set the initial degree
161 - -eps_chase_degree_opt - toggle the optimization
162
163 Note:
164 See the documentation of ChASE {cite:p}`Win19` for details.
165
166 Level: advanced
167
168 .seealso: [](ch:eps), `EPSCHASE`, `EPSCHASEGetDegree()`
169 @*/
170 6 PetscErrorCode EPSCHASESetDegree(EPS eps,PetscInt deg,PetscBool opt)
171 {
172 6 PetscFunctionBegin;
173 6 PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
174 6 PetscValidLogicalCollectiveInt(eps,deg,2);
175 6 PetscValidLogicalCollectiveBool(eps,opt,3);
176
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 6 PetscTryMethod(eps,"EPSCHASESetDegree_C",(EPS,PetscInt,PetscBool),(eps,deg,opt));
177 6 PetscFunctionReturn(PETSC_SUCCESS);
178 }
179
180 6 static PetscErrorCode EPSCHASEGetDegree_ChASE(EPS eps,PetscInt *deg,PetscBool *opt)
181 {
182 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
183
184 6 PetscFunctionBegin;
185
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 6 if (deg) *deg = ctx->deg;
186
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 6 if (opt) *opt = ctx->opt;
187 6 PetscFunctionReturn(PETSC_SUCCESS);
188 }
189
190 /*@
191 EPSCHASEGetDegree - Gets the degree of the Chebyshev polynomial filter used in the ChASE solver.
192
193 Not Collective
194
195 Input Parameter:
196 . eps - the linear eigensolver context
197
198 Output Parameters:
199 + deg - initial degree of Chebyshev polynomial filter
200 - opt - internal optimization of polynomial degree
201
202 Level: advanced
203
204 .seealso: [](ch:eps), `EPSCHASE`, `EPSCHASESetDegree()`
205 @*/
206 6 PetscErrorCode EPSCHASEGetDegree(EPS eps,PetscInt *deg,PetscBool *opt)
207 {
208 6 PetscFunctionBegin;
209 6 PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
210
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 6 PetscUseMethod(eps,"EPSCHASEGetDegree_C",(EPS,PetscInt*,PetscBool*),(eps,deg,opt));
211 6 PetscFunctionReturn(PETSC_SUCCESS);
212 }
213
214 6 static PetscErrorCode EPSSetFromOptions_ChASE(EPS eps,PetscOptionItems PetscOptionsObject)
215 {
216 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
217 6 PetscBool flg1,flg2;
218 6 PetscInt deg;
219 6 PetscBool opt;
220
221 6 PetscFunctionBegin;
222
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 6 PetscOptionsHeadBegin(PetscOptionsObject,"EPS ChASE Options");
223
224 6 deg = ctx->deg;
225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscOptionsInt("-eps_chase_degree","Initial degree of Chebyshev polynomial filter","EPSCHASESetDegree",deg,&deg,&flg1));
226 6 opt = ctx->opt;
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscOptionsBool("-eps_chase_degree_opt","Internal optimization of polynomial degree","EPSCHASESetDegree",opt,&opt,&flg2));
228
2/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 6 if (flg1 || flg2) PetscCall(EPSCHASESetDegree(eps,deg,opt));
229
230 6 PetscOptionsHeadEnd();
231 PetscFunctionReturn(PETSC_SUCCESS);
232 }
233
234 6 static PetscErrorCode EPSDestroy_ChASE(EPS eps)
235 {
236 6 PetscFunctionBegin;
237
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 6 PetscCall(PetscFree(eps->data));
238
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",NULL));
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",NULL));
240 PetscFunctionReturn(PETSC_SUCCESS);
241 }
242
243 6 static PetscErrorCode EPSReset_ChASE(EPS eps)
244 {
245 6 EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
246
247 6 PetscFunctionBegin;
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(MatDestroy(&ctx->As));
249 PetscFunctionReturn(PETSC_SUCCESS);
250 }
251
252 ✗ static PetscErrorCode EPSView_ChASE(EPS eps,PetscViewer viewer)
253 {
254 ✗ EPS_ChASE *ctx = (EPS_ChASE*)eps->data;
255 ✗ PetscBool isascii;
256
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 }
265
266 /*MC
267 EPSCHASE - EPSCHASE = "chase" - A wrapper to CHASE {cite:p}`Win19`.
268
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.
273
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.
277
278 Level: beginner
279
280 .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
281 M*/
282 6 SLEPC_EXTERN PetscErrorCode EPSCreate_ChASE(EPS eps)
283 {
284 6 EPS_ChASE *ctx;
285
286 6 PetscFunctionBegin;
287
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscNew(&ctx));
288 6 eps->data = (void*)ctx;
289 6 ctx->deg = PETSC_DEFAULT;
290 6 ctx->opt = PETSC_TRUE;
291
292 6 eps->categ = EPS_CATEGORY_OTHER;
293
294 6 eps->ops->solve = EPSSolve_ChASE;
295 6 eps->ops->setup = EPSSetUp_ChASE;
296 6 eps->ops->setupsort = EPSSetUpSort_Basic;
297 6 eps->ops->setfromoptions = EPSSetFromOptions_ChASE;
298 6 eps->ops->destroy = EPSDestroy_ChASE;
299 6 eps->ops->reset = EPSReset_ChASE;
300 6 eps->ops->view = EPSView_ChASE;
301 6 eps->ops->backtransform = EPSBackTransform_Default;
302 6 eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
303
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASESetDegree_C",EPSCHASESetDegree_ChASE));
305
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 6 PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSCHASEGetDegree_C",EPSCHASEGetDegree_ChASE));
306 PetscFunctionReturn(PETSC_SUCCESS);
307 }
308