GCC Code Coverage Report

Directory: ./
File: src/svd/impls/external/scalapack/svdscalap.c
Date: 2025-12-10 04:20:18
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 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 info,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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 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 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 30 PetscCheckScaLapackInfo("gesvd",info);
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
80
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 PetscCall(PetscMalloc1(lwork,&work));
81 /* call computational routine */
82
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 30 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 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 30 PetscCheckScaLapackInfo("gesvd",info);
84
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 30 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 4 times.
 40 PetscCall(PetscFPTrapPop());
98
99
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatHermitianTranspose(QT,MAT_INITIAL_MATRIX,&Q));
100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&QT));
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVGetMat(svd->U,&U));
102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVGetMat(svd->V,&V));
103
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 40 if (M>=N) {
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 36 PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 36 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
106 } else {
107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
109 }
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVRestoreMat(svd->U,&U));
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(BVRestoreMat(svd->V,&V));
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&Z));
113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(MatDestroy(&Q));
114
115 40 svd->nconv = svd->ncv;
116 40 svd->its = 1;
117 40 svd->reason = SVD_CONVERGED_TOL;
118 40 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
121 40 static PetscErrorCode SVDDestroy_ScaLAPACK(SVD svd)
122 {
123 40 PetscFunctionBegin;
124
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));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 40 static PetscErrorCode SVDReset_ScaLAPACK(SVD svd)
129 {
130 40 SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
131
132 40 PetscFunctionBegin;
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 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 40 SLEPC_EXTERN PetscErrorCode SVDCreate_ScaLAPACK(SVD svd)
156 {
157 40 SVD_ScaLAPACK *ctx;
158
159 40 PetscFunctionBegin;
160
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 40 PetscCall(PetscNew(&ctx));
161 40 svd->data = (void*)ctx;
162
163 40 svd->ops->solve = SVDSolve_ScaLAPACK;
164 40 svd->ops->setup = SVDSetUp_ScaLAPACK;
165 40 svd->ops->destroy = SVDDestroy_ScaLAPACK;
166 40 svd->ops->reset = SVDReset_ScaLAPACK;
167 40 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169