GCC Code Coverage Report

Directory: ./
File: src/svd/impls/external/scalapack/svdscalap.c
Date: 2026-07-10 04:05:41
Exec Total Coverage
Lines: 88 88 100.0%
Functions: 5 5 100.0%
Branches: 64 146 43.8%
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 the ScaLAPACK SVD solver
12 */
13
14 #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
15 #include <slepc/private/slepcscalapack.h>
16
17 typedef struct {
18 Mat As; /* converted matrix */
19 } SVD_ScaLAPACK;
20
21 40 static PetscErrorCode SVDSetUp_ScaLAPACK(SVD svd)
22 {
23 40 SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
24 40 PetscInt M,N;
25
26 40 PetscFunctionBegin;
27
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 SVDCheckStandard(svd);
28
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 SVDCheckDefinite(svd);
29
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 40 if (svd->nsv==0) svd->nsv = 1;
30
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatGetSize(svd->A,&M,&N));
31 40 svd->ncv = N;
32
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
33
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 40 if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
34 40 svd->leftbasis = PETSC_TRUE;
35
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
36
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(SVDAllocateSolution(svd,0));
37
38 /* convert matrix */
39
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&ctx->As));
40
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatConvert(svd->OP,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
41 PetscFunctionReturn(PETSC_SUCCESS);
42 }
43
44 40 static PetscErrorCode SVDSolve_ScaLAPACK(SVD svd)
45 {
46 40 SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
47 40 Mat A = ctx->As,Z,Q,QT,U,V;
48 40 Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)A->data,*q,*z;
49 40 PetscScalar *work,minlwork;
50 40 PetscBLASInt lwork=-1,one=1;
51 40 PetscInt M,N,m,n,mn;
52 #if defined(PETSC_USE_COMPLEX)
53 10 PetscBLASInt lrwork;
54 10 PetscReal *rwork,dummy;
55 #endif
56
57 40 PetscFunctionBegin;
58
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatGetSize(A,&M,&N));
59
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatGetLocalSize(A,&m,&n));
60
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 40 mn = (M>=N)? n: m;
61
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Z));
62
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatSetSizes(Z,m,mn,PETSC_DECIDE,PETSC_DECIDE));
63
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatSetType(Z,MATSCALAPACK));
64
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatAssemblyBegin(Z,MAT_FINAL_ASSEMBLY));
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatAssemblyEnd(Z,MAT_FINAL_ASSEMBLY));
66 40 z = (Mat_ScaLAPACK*)Z->data;
67
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&QT));
68
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatSetSizes(QT,mn,n,PETSC_DECIDE,PETSC_DECIDE));
69
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatSetType(QT,MATSCALAPACK));
70
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatAssemblyBegin(QT,MAT_FINAL_ASSEMBLY));
71
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatAssemblyEnd(QT,MAT_FINAL_ASSEMBLY));
72 40 q = (Mat_ScaLAPACK*)QT->data;
73
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
75 #if !defined(PETSC_USE_COMPLEX)
76 /* allocate workspace */
77
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.
 30 PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,&minlwork,&lwork,&info));
78
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 PetscCall(PetscMalloc1(lwork,&work));
80 /* call computational routine */
81
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.
 30 PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,work,&lwork,&info));
82
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 30 PetscCall(PetscFree(work));
83 #else
84 /* allocate workspace */
85
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.
 10 PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,&minlwork,&lwork,&dummy,&info));
86
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 10 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork),&lwork));
87 10 lrwork = 1+4*PetscMax(a->M,a->N);
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 10 PetscCall(PetscMalloc2(lwork,&work,lrwork,&rwork));
89 /* call computational routine */
90
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.
 10 PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,work,&lwork,rwork,&info));
91
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 10 PetscCall(PetscFree2(work,rwork));
92 #endif
93
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(PetscFPTrapPop());
94
95
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatHermitianTranspose(QT,MAT_INITIAL_MATRIX,&Q));
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&QT));
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVGetMat(svd->U,&U));
98
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVGetMat(svd->V,&V));
99
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 40 if (M>=N) {
100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 36 PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 36 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
102 } else {
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
105 }
106
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVRestoreMat(svd->U,&U));
107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVRestoreMat(svd->V,&V));
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&Z));
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&Q));
110
111 40 svd->nconv = svd->ncv;
112 40 svd->its = 1;
113 40 svd->reason = SVD_CONVERGED_TOL;
114 40 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
117 40 static PetscErrorCode SVDDestroy_ScaLAPACK(SVD svd)
118 {
119 40 PetscFunctionBegin;
120
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 40 PetscCall(PetscFree(svd->data));
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
124 40 static PetscErrorCode SVDReset_ScaLAPACK(SVD svd)
125 {
126 40 SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
127
128 40 PetscFunctionBegin;
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&ctx->As));
130 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 /*MC
134 SVDSCALAPACK - SVDSCALAPACK = "scalapack" - A wrapper to the ScaLAPACK
135 singular value solver {cite:p}`Bla97`.
136
137 Notes:
138 Only available for standard SVD problems, using subroutine `pdgesvd` and
139 the analogs for other precisions.
140
141 This is a direct singular value solver, that is, the full decomposition
142 is computed. The computation involves redistributing the matrices from PETSc
143 storage to ScaLAPACK distribution, and vice versa (this is done automatically
144 by SLEPc). Alternatively, the user may create the problem matrices
145 already with type `MATSCALAPACK`.
146
147 Level: beginner
148
149 .seealso: [](ch:svd), `SVD`, `SVDType`, `SVDSetType()`
150 M*/
151 40 SLEPC_EXTERN PetscErrorCode SVDCreate_ScaLAPACK(SVD svd)
152 {
153 40 SVD_ScaLAPACK *ctx;
154
155 40 PetscFunctionBegin;
156
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(PetscNew(&ctx));
157 40 svd->data = (void*)ctx;
158
159 40 svd->ops->solve = SVDSolve_ScaLAPACK;
160 40 svd->ops->setup = SVDSetUp_ScaLAPACK;
161 40 svd->ops->destroy = SVDDestroy_ScaLAPACK;
162 40 svd->ops->reset = SVDReset_ScaLAPACK;
163 40 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165