GCC Code Coverage Report

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