GCC Code Coverage Report

Directory: ./
File: src/eps/impls/external/scalapack/scalapack.c
Date: 2026-07-06 04:10:30
Exec Total Coverage
Lines: 106 108 98.1%
Functions: 5 5 100.0%
Branches: 92 232 39.7%
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 ScaLAPACK.
12 */
13
14 #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 #include <slepc/private/slepcscalapack.h>
16
17 typedef struct {
18 Mat As,Bs; /* converted matrices */
19 } EPS_ScaLAPACK;
20
21 44 static PetscErrorCode EPSSetUp_ScaLAPACK(EPS eps)
22 {
23 44 EPS_ScaLAPACK *ctx = (EPS_ScaLAPACK*)eps->data;
24 44 Mat A,B;
25 44 PetscInt nmat;
26 44 PetscBool isshift;
27 44 PetscScalar shift;
28
29 44 PetscFunctionBegin;
30
4/10
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 44 EPSCheckHermitianDefinite(eps);
31
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 44 EPSCheckNotStructured(eps);
32
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
33
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 44 PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
34
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 if (eps->nev==0) eps->nev = 1;
35 44 eps->ncv = eps->n;
36
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 44 if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
37
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 44 if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
38
3/4
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 44 if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
39
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 44 PetscCheck(eps->which!=EPS_ALL || eps->inta==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support interval computation");
40
4/14
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 44 EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION);
41
5/12
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 44 EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
42
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(EPSAllocateSolution(eps,0));
43
44 /* convert matrices */
45
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatDestroy(&ctx->As));
46
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatDestroy(&ctx->Bs));
47
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(STGetNumMatrices(eps->st,&nmat));
48
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(STGetMatrix(eps->st,0,&A));
49
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
50
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 44 if (nmat>1) {
51
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 24 PetscCall(STGetMatrix(eps->st,1,&B));
52
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 24 PetscCall(MatConvert(B,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->Bs));
53 }
54
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(STGetShift(eps->st,&shift));
55
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 if (shift != 0.0) {
56 ✗ if (nmat>1) PetscCall(MatAXPY(ctx->As,-shift,ctx->Bs,SAME_NONZERO_PATTERN));
57 ✗ else PetscCall(MatShift(ctx->As,-shift));
58 }
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
62 44 static PetscErrorCode EPSSolve_ScaLAPACK(EPS eps)
63 {
64 44 EPS_ScaLAPACK *ctx = (EPS_ScaLAPACK*)eps->data;
65 44 Mat A = ctx->As,B = ctx->Bs,Q,V;
66 44 Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)A->data,*b,*q;
67 44 PetscReal rdummy=0.0,abstol=0.0,*gap=NULL,orfac=-1.0,*w = eps->errest; /* used to store real eigenvalues */
68 44 PetscScalar *work,minlwork[3];
69 44 PetscBLASInt i,m,idummy=0,lwork=-1,liwork=-1,minliwork,*iwork,*ifail=NULL,*iclustr=NULL,one=1;
70 #if defined(PETSC_USE_COMPLEX)
71 11 PetscReal *rwork,minlrwork[3];
72 11 PetscBLASInt lrwork=-1;
73 #endif
74
75 44 PetscFunctionBegin;
76
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
77
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
78 44 q = (Mat_ScaLAPACK*)Q->data;
79
80
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 44 if (B) {
81
82 24 b = (Mat_ScaLAPACK*)B->data;
83
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 24 PetscCall(PetscMalloc3(a->grid->nprow*a->grid->npcol,&gap,a->N,&ifail,2*a->grid->nprow*a->grid->npcol,&iclustr));
84 #if !defined(PETSC_USE_COMPLEX)
85 /* allocate workspace */
86
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 18 PetscCallScaLAPACKInfo("sygvx",SCALAPACKsygvx_(&one,"V","A","L",&a->N,a->loc,&one,&one,a->desc,b->loc,&one,&one,b->desc,&rdummy,&rdummy,&idummy,&idummy,&abstol,&m,&idummy,w,&orfac,q->loc,&one,&one,q->desc,minlwork,&lwork,&minliwork,&liwork,ifail,iclustr,gap,&info));
87
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscBLASIntCast((PetscInt)minlwork[0],&lwork));
88 18 liwork = minliwork;
89 /* call computational routine */
90
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscMalloc2(lwork,&work,liwork,&iwork));
91
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 18 PetscCallScaLAPACKInfo("sygvx",SCALAPACKsygvx_(&one,"V","A","L",&a->N,a->loc,&one,&one,a->desc,b->loc,&one,&one,b->desc,&rdummy,&rdummy,&idummy,&idummy,&abstol,&m,&idummy,w,&orfac,q->loc,&one,&one,q->desc,work,&lwork,iwork,&liwork,ifail,iclustr,gap,&info));
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscFree2(work,iwork));
93 #else
94 /* allocate workspace */
95
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 6 PetscCallScaLAPACKInfo("sygvx",SCALAPACKsygvx_(&one,"V","A","L",&a->N,a->loc,&one,&one,a->desc,b->loc,&one,&one,b->desc,&rdummy,&rdummy,&idummy,&idummy,&abstol,&m,&idummy,w,&orfac,q->loc,&one,&one,q->desc,minlwork,&lwork,minlrwork,&lrwork,&minliwork,&liwork,ifail,iclustr,gap,&info));
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork[0]),&lwork));
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscBLASIntCast((PetscInt)minlrwork[0],&lrwork));
98 6 lrwork += a->N*a->N;
99 6 liwork = minliwork;
100 /* call computational routine */
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscMalloc3(lwork,&work,lrwork,&rwork,liwork,&iwork));
102
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 6 PetscCallScaLAPACKInfo("sygvx",SCALAPACKsygvx_(&one,"V","A","L",&a->N,a->loc,&one,&one,a->desc,b->loc,&one,&one,b->desc,&rdummy,&rdummy,&idummy,&idummy,&abstol,&m,&idummy,w,&orfac,q->loc,&one,&one,q->desc,work,&lwork,rwork,&lrwork,iwork,&liwork,ifail,iclustr,gap,&info));
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscFree3(work,rwork,iwork));
104 #endif
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 24 PetscCall(PetscFree3(gap,ifail,iclustr));
106
107 } else {
108
109 #if !defined(PETSC_USE_COMPLEX)
110 /* allocate workspace */
111
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 15 PetscCallScaLAPACKInfo("syev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,minlwork,&lwork,&info));
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCall(PetscBLASIntCast((PetscInt)minlwork[0],&lwork));
113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCall(PetscMalloc1(lwork,&work));
114 /* call computational routine */
115
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 15 PetscCallScaLAPACKInfo("syev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,work,&lwork,&info));
116
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 15 PetscCall(PetscFree(work));
117 #else
118 /* allocate workspace */
119
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 5 PetscCallScaLAPACKInfo("syev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,minlwork,&lwork,minlrwork,&lrwork,&info));
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork[0]),&lwork));
121 5 lrwork = 4*a->N; /* PetscCall(PetscBLASIntCast((PetscInt)minlrwork[0],&lrwork)); */
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscMalloc2(lwork,&work,lrwork,&rwork));
123 /* call computational routine */
124
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 5 PetscCallScaLAPACKInfo("syev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,work,&lwork,rwork,&lrwork,&info));
125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscFree2(work,rwork));
126 #endif
127
128 }
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(PetscFPTrapPop());
130
131
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8420 for (i=0;i<eps->ncv;i++) {
132 8376 eps->eigr[i] = eps->errest[i];
133 8376 eps->errest[i] = PETSC_MACHINE_EPSILON;
134 }
135
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(BVGetMat(eps->V,&V));
137
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(BVRestoreMat(eps->V,&V));
139
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatDestroy(&Q));
140
141 44 eps->nconv = eps->ncv;
142 44 eps->its = 1;
143 44 eps->reason = EPS_CONVERGED_TOL;
144 44 PetscFunctionReturn(PETSC_SUCCESS);
145 }
146
147 32 static PetscErrorCode EPSDestroy_ScaLAPACK(EPS eps)
148 {
149 32 PetscFunctionBegin;
150
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 32 PetscCall(PetscFree(eps->data));
151 PetscFunctionReturn(PETSC_SUCCESS);
152 }
153
154 32 static PetscErrorCode EPSReset_ScaLAPACK(EPS eps)
155 {
156 32 EPS_ScaLAPACK *ctx = (EPS_ScaLAPACK*)eps->data;
157
158 32 PetscFunctionBegin;
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(MatDestroy(&ctx->As));
160
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(MatDestroy(&ctx->Bs));
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
164 /*MC
165 EPSSCALAPACK - EPSSCALAPACK = "scalapack" - A wrapper to ScaLAPACK
166 eigensolvers {cite:p}`Bla97`.
167
168 Notes:
169 Only available for Hermitian problems, using subroutines `pdsyev` and
170 `pdsygvx` and the analogs for other precisions.
171
172 This is a direct eigensolver, that is, the full spectrum is computed.
173 The computation involves redistributing the matrices from PETSc storage
174 to ScaLAPACK distribution, and vice versa (this is done automatically
175 by SLEPc). Alternatively, the user may create the problem matrices
176 already with type `MATSCALAPACK`.
177
178 Level: beginner
179
180 .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
181 M*/
182 32 SLEPC_EXTERN PetscErrorCode EPSCreate_ScaLAPACK(EPS eps)
183 {
184 32 EPS_ScaLAPACK *ctx;
185
186 32 PetscFunctionBegin;
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(PetscNew(&ctx));
188 32 eps->data = (void*)ctx;
189
190 32 eps->categ = EPS_CATEGORY_OTHER;
191
192 32 eps->ops->solve = EPSSolve_ScaLAPACK;
193 32 eps->ops->setup = EPSSetUp_ScaLAPACK;
194 32 eps->ops->setupsort = EPSSetUpSort_Basic;
195 32 eps->ops->destroy = EPSDestroy_ScaLAPACK;
196 32 eps->ops->reset = EPSReset_ScaLAPACK;
197 32 eps->ops->backtransform = EPSBackTransform_Default;
198 32 eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
199 32 PetscFunctionReturn(PETSC_SUCCESS);
200 }
201