GCC Code Coverage Report

Directory: ./
File: src/eps/impls/external/scalapack/scalapack.c
Date: 2025-11-19 04:19:03
Exec Total Coverage
Lines: 114 116 98.3%
Functions: 5 5 100.0%
Branches: 84 200 42.0%
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,info,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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCallBLAS("SCALAPACKsygvx",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/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 18 PetscCheckScaLapackInfo("sygvx",info);
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscBLASIntCast((PetscInt)minlwork[0],&lwork));
89 18 liwork = minliwork;
90 /* call computational routine */
91
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscMalloc2(lwork,&work,liwork,&iwork));
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCallBLAS("SCALAPACKsygvx",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));
93
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 18 PetscCheckScaLapackInfo("sygvx",info);
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 18 PetscCall(PetscFree2(work,iwork));
95 #else
96 /* allocate workspace */
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCallBLAS("SCALAPACKsygvx",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));
98
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheckScaLapackInfo("sygvx",info);
99
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork[0]),&lwork));
100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscBLASIntCast((PetscInt)minlrwork[0],&lrwork));
101 6 lrwork += a->N*a->N;
102 6 liwork = minliwork;
103 /* call computational routine */
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscMalloc3(lwork,&work,lrwork,&rwork,liwork,&iwork));
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCallBLAS("SCALAPACKsygvx",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));
106
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 PetscCheckScaLapackInfo("sygvx",info);
107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 6 PetscCall(PetscFree3(work,rwork,iwork));
108 #endif
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 24 PetscCall(PetscFree3(gap,ifail,iclustr));
110
111 } else {
112
113 #if !defined(PETSC_USE_COMPLEX)
114 /* allocate workspace */
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCallBLAS("SCALAPACKsyev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,minlwork,&lwork,&info));
116
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 15 PetscCheckScaLapackInfo("syev",info);
117
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCall(PetscBLASIntCast((PetscInt)minlwork[0],&lwork));
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCall(PetscMalloc1(lwork,&work));
119 /* call computational routine */
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 15 PetscCallBLAS("SCALAPACKsyev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,work,&lwork,&info));
121
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 15 PetscCheckScaLapackInfo("syev",info);
122
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 15 PetscCall(PetscFree(work));
123 #else
124 /* allocate workspace */
125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCallBLAS("SCALAPACKsyev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,minlwork,&lwork,minlrwork,&lrwork,&info));
126
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 5 PetscCheckScaLapackInfo("syev",info);
127
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork[0]),&lwork));
128 5 lrwork = 4*a->N; /* PetscCall(PetscBLASIntCast((PetscInt)minlrwork[0],&lrwork)); */
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscMalloc2(lwork,&work,lrwork,&rwork));
130 /* call computational routine */
131
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCallBLAS("SCALAPACKsyev",SCALAPACKsyev_("V","L",&a->N,a->loc,&one,&one,a->desc,w,q->loc,&one,&one,q->desc,work,&lwork,rwork,&lrwork,&info));
132
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 5 PetscCheckScaLapackInfo("syev",info);
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 5 PetscCall(PetscFree2(work,rwork));
134 #endif
135
136 }
137
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(PetscFPTrapPop());
138
139
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8420 for (i=0;i<eps->ncv;i++) {
140 8376 eps->eigr[i] = eps->errest[i];
141 8376 eps->errest[i] = PETSC_MACHINE_EPSILON;
142 }
143
144
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(BVGetMat(eps->V,&V));
145
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
146
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(BVRestoreMat(eps->V,&V));
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 44 PetscCall(MatDestroy(&Q));
148
149 44 eps->nconv = eps->ncv;
150 44 eps->its = 1;
151 44 eps->reason = EPS_CONVERGED_TOL;
152 44 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
155 32 static PetscErrorCode EPSDestroy_ScaLAPACK(EPS eps)
156 {
157 32 PetscFunctionBegin;
158
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));
159 PetscFunctionReturn(PETSC_SUCCESS);
160 }
161
162 32 static PetscErrorCode EPSReset_ScaLAPACK(EPS eps)
163 {
164 32 EPS_ScaLAPACK *ctx = (EPS_ScaLAPACK*)eps->data;
165
166 32 PetscFunctionBegin;
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(MatDestroy(&ctx->As));
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(MatDestroy(&ctx->Bs));
169 PetscFunctionReturn(PETSC_SUCCESS);
170 }
171
172 /*MC
173 EPSSCALAPACK - EPSSCALAPACK = "scalapack" - A wrapper to ScaLAPACK
174 eigensolvers {cite:p}`Bla97`.
175
176 Notes:
177 Only available for Hermitian problems, using subroutines `pdsyev` and
178 `pdsygvx` and the analogs for other precisions.
179
180 This is a direct eigensolver, that is, the full spectrum is computed.
181 The computation involves redistributing the matrices from PETSc storage
182 to ScaLAPACK distribution, and vice versa (this is done automatically
183 by SLEPc). Alternatively, the user may create the problem matrices
184 already with type `MATSCALAPACK`.
185
186 Level: beginner
187
188 .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
189 M*/
190 32 SLEPC_EXTERN PetscErrorCode EPSCreate_ScaLAPACK(EPS eps)
191 {
192 32 EPS_ScaLAPACK *ctx;
193
194 32 PetscFunctionBegin;
195
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 32 PetscCall(PetscNew(&ctx));
196 32 eps->data = (void*)ctx;
197
198 32 eps->categ = EPS_CATEGORY_OTHER;
199
200 32 eps->ops->solve = EPSSolve_ScaLAPACK;
201 32 eps->ops->setup = EPSSetUp_ScaLAPACK;
202 32 eps->ops->setupsort = EPSSetUpSort_Basic;
203 32 eps->ops->destroy = EPSDestroy_ScaLAPACK;
204 32 eps->ops->reset = EPSReset_ScaLAPACK;
205 32 eps->ops->backtransform = EPSBackTransform_Default;
206 32 eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
207 32 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209