BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_gesvd.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California, through
4 * Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
21
22#ifndef BELFEM_FN_GESVD_HPP
23#define BELFEM_FN_GESVD_HPP
24
25#include "assert.hpp"
26#include "lapacktools.hpp"
27#include "cl_Vector.hpp"
28
29namespace belfem
30{
31 namespace lapack
32 {
33#ifdef __cplusplus
34 extern "C"
35 {
36#endif
37// -----------------------------------------------------------------------------
38
39 // two gesvd particularities: the singular values s ( and the
40 // complex-only scratch rwork ) stay REAL valued in all flavors,
41 // and the two job chars carry two hidden Fortran length arguments
42
43 void
44 sgesvd_(
45 char * jobu,
46 char * jobvt,
47 int_t * m,
48 int_t * n,
49 float * a,
50 int_t * lda,
51 float * s,
52 float * u,
53 int_t * ldu,
54 float * vt,
55 int_t * ldvt,
56 float * work,
57 int_t * lwork,
58 int_t * info,
59 fortran_charlen_t lu,
60 fortran_charlen_t lvt );
61
62// -----------------------------------------------------------------------------
63
64 void
65 dgesvd_(
66 char * jobu,
67 char * jobvt,
68 int_t * m,
69 int_t * n,
70 double * a,
71 int_t * lda,
72 double * s,
73 double * u,
74 int_t * ldu,
75 double * vt,
76 int_t * ldvt,
77 double * work,
78 int_t * lwork,
79 int_t * info,
80 fortran_charlen_t lu,
81 fortran_charlen_t lvt );
82
83// -----------------------------------------------------------------------------
84
85 void
86 cgesvd_(
87 char * jobu,
88 char * jobvt,
89 int_t * m,
90 int_t * n,
91 cplx_float_t * a,
92 int_t * lda,
93 float * s,
94 cplx_float_t * u,
95 int_t * ldu,
96 cplx_float_t * vt,
97 int_t * ldvt,
98 cplx_float_t * work,
99 int_t * lwork,
100 float * rwork,
101 int_t * info,
102 fortran_charlen_t lu,
103 fortran_charlen_t lvt );
104
105// -----------------------------------------------------------------------------
106
107 void
108 zgesvd_(
109 char * jobu,
110 char * jobvt,
111 int_t * m,
112 int_t * n,
113 cplx_double_t * a,
114 int_t * lda,
115 double * s,
116 cplx_double_t * u,
117 int_t * ldu,
118 cplx_double_t * vt,
119 int_t * ldvt,
120 cplx_double_t * work,
121 int_t * lwork,
122 double * rwork,
123 int_t * info,
124 fortran_charlen_t lu,
125 fortran_charlen_t lvt );
126
127#ifdef __cplusplus
128 }
129#endif
130// -----------------------------------------------------------------------------
131
132 // unified dispatch: rwork is referenced by the complex flavors only,
133 // the real specializations ignore it
134
135 template< typename T >
136 void gesvd(
137 const char * jobu,
138 const char * jobvt,
139 const int_t * m,
140 const int_t * n,
141 T * a,
142 const int_t * lda,
143 real_t< T > * s,
144 T * u,
145 const int_t * ldu,
146 T * vt,
147 const int_t * ldvt,
148 T * work,
149 const int_t * lwork,
150 real_t< T > * rwork,
151 int_t * info )
152 {
153 static_assert( dependent_false< T >,
154 "gesvd not implemented for selected data type" );
155 }
156
157// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158
159 template<>
160 inline void gesvd(
161 const char * jobu,
162 const char * jobvt,
163 const int_t * m,
164 const int_t * n,
165 float * a,
166 const int_t * lda,
167 float * s,
168 float * u,
169 const int_t * ldu,
170 float * vt,
171 const int_t * ldvt,
172 float * work,
173 const int_t * lwork,
174 float * /* rwork */,
175 int_t * info )
176 {
177 sgesvd_(
178 const_cast< char * >( jobu ),
179 const_cast< char * >( jobvt ),
180 const_cast< int_t * >( m ),
181 const_cast< int_t * >( n ),
182 a,
183 const_cast< int_t * >( lda ),
184 s,
185 u,
186 const_cast< int_t * >( ldu ),
187 vt,
188 const_cast< int_t * >( ldvt ),
189 work,
190 const_cast< int_t * >( lwork ),
191 info,
192 1, 1 );
193 }
194
195// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196
197 template<>
198 inline void gesvd(
199 const char * jobu,
200 const char * jobvt,
201 const int_t * m,
202 const int_t * n,
203 double * a,
204 const int_t * lda,
205 double * s,
206 double * u,
207 const int_t * ldu,
208 double * vt,
209 const int_t * ldvt,
210 double * work,
211 const int_t * lwork,
212 double * /* rwork */,
213 int_t * info )
214 {
215 dgesvd_(
216 const_cast< char * >( jobu ),
217 const_cast< char * >( jobvt ),
218 const_cast< int_t * >( m ),
219 const_cast< int_t * >( n ),
220 a,
221 const_cast< int_t * >( lda ),
222 s,
223 u,
224 const_cast< int_t * >( ldu ),
225 vt,
226 const_cast< int_t * >( ldvt ),
227 work,
228 const_cast< int_t * >( lwork ),
229 info,
230 1, 1 );
231 }
232
233// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234
235 template<>
236 inline void gesvd(
237 const char * jobu,
238 const char * jobvt,
239 const int_t * m,
240 const int_t * n,
241 std::complex< float > * a,
242 const int_t * lda,
243 float * s,
244 std::complex< float > * u,
245 const int_t * ldu,
246 std::complex< float > * vt,
247 const int_t * ldvt,
248 std::complex< float > * work,
249 const int_t * lwork,
250 float * rwork,
251 int_t * info )
252 {
253 cgesvd_(
254 const_cast< char * >( jobu ),
255 const_cast< char * >( jobvt ),
256 const_cast< int_t * >( m ),
257 const_cast< int_t * >( n ),
258 reinterpret_cast< cplx_float_t * >( a ),
259 const_cast< int_t * >( lda ),
260 s,
261 reinterpret_cast< cplx_float_t * >( u ),
262 const_cast< int_t * >( ldu ),
263 reinterpret_cast< cplx_float_t * >( vt ),
264 const_cast< int_t * >( ldvt ),
265 reinterpret_cast< cplx_float_t * >( work ),
266 const_cast< int_t * >( lwork ),
267 rwork,
268 info,
269 1, 1 );
270 }
271
272// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273
274 template<>
275 inline void gesvd(
276 const char * jobu,
277 const char * jobvt,
278 const int_t * m,
279 const int_t * n,
280 std::complex< double > * a,
281 const int_t * lda,
282 double * s,
283 std::complex< double > * u,
284 const int_t * ldu,
285 std::complex< double > * vt,
286 const int_t * ldvt,
287 std::complex< double > * work,
288 const int_t * lwork,
289 double * rwork,
290 int_t * info )
291 {
292 zgesvd_(
293 const_cast< char * >( jobu ),
294 const_cast< char * >( jobvt ),
295 const_cast< int_t * >( m ),
296 const_cast< int_t * >( n ),
297 reinterpret_cast< cplx_double_t * >( a ),
298 const_cast< int_t * >( lda ),
299 s,
300 reinterpret_cast< cplx_double_t * >( u ),
301 const_cast< int_t * >( ldu ),
302 reinterpret_cast< cplx_double_t * >( vt ),
303 const_cast< int_t * >( ldvt ),
304 reinterpret_cast< cplx_double_t * >( work ),
305 const_cast< int_t * >( lwork ),
306 rwork,
307 info,
308 1, 1 );
309 }
310
311// -----------------------------------------------------------------------------
312 } /* end namespace lapack */
313
338 template< typename T >
339 int_t
340 gesvd(
341 Matrix< T > & A,
342 Vector< lapack::real_t< T > > & S,
343 Matrix< T > & U,
344 Matrix< T > & VT,
345 Vector< lapack::real_t< T > > & Work,
346 const char jobu = 'A',
347 const char jobvt = 'A',
348 const bool AbortOnError = true )
349 {
350 BELFEM_ASSERT( jobu == 'A' || jobu == 'S' || jobu == 'N',
351 "unsupported jobu flag '%c'", jobu );
352 BELFEM_ASSERT( jobvt == 'A' || jobvt == 'S' || jobvt == 'N',
353 "unsupported jobvt flag '%c'", jobvt );
354
355 // reals per LAPACK work entry
356 constexpr int_t tRealsPerT =
357 std::is_same< T, lapack::real_t< T > >::value ? 1 : 2 ;
358
359 int_t m = ( int_t ) A.n_rows();
360 int_t n = ( int_t ) A.n_cols();
361 int_t mn = std::min( m, n );
362
363 int_t lda = lapack::leading_dimension( A );
364
365 S.set_size( mn );
366
367 // U is m x m for 'A' and m x mn for 'S'; lapack requires
368 // ldu >= 1 even when U is not referenced ( jobu = 'N' )
369 int_t ldu = 1 ;
370 if ( jobu != 'N' )
371 {
372 U.set_size( m, jobu == 'A' ? m : mn );
373 ldu = lapack::leading_dimension( U );
374 }
375
376 // VT is n x n for 'A' and mn x n for 'S'
377 int_t ldvt = 1 ;
378 if ( jobvt != 'N' )
379 {
380 VT.set_size( jobvt == 'A' ? n : mn, n );
381 ldvt = lapack::leading_dimension( VT );
382 }
383
384 // real scratch in the tail, referenced by the complex flavors only
385 int_t tRWorkSize = tRealsPerT == 1 ? 0 : 5 * mn ;
386
387 // minimum work sizes: real flavors need
388 // max( 3*mn + max( m, n ), 5*mn ), complex ones 2*mn + max( m, n )
389 int_t lwork = std::max< int_t >( 1, tRealsPerT == 1 ?
390 std::max( 3 * mn + std::max( m, n ), 5 * mn ) :
391 2 * mn + std::max( m, n ) );
392
393 int_t info = 0 ;
394
395 // required buffer: lwork entries of T in the head plus the rwork tail
396 if ( static_cast< int_t >( Work.length() ) < tRealsPerT * lwork + tRWorkSize )
397 {
398 // ask lapack for the optimal size ( one T entry plus the
399 // tail ); this leaves A untouched. Grow-only: never shrink a
400 // buffer that already covers the query call
401 if ( static_cast< int_t >( Work.length() ) < tRealsPerT + tRWorkSize )
402 {
403 Work.set_size( tRealsPerT + tRWorkSize );
404 }
405
406 int_t query = -1 ;
407
408 lapack::gesvd( &jobu, &jobvt, &m, &n, A.data(), &lda, S.data(),
409 U.data(), &ldu, VT.data(), &ldvt,
410 reinterpret_cast< T * >( Work.data() ), &query,
411 Work.data() + tRealsPerT, &info );
412
413 BELFEM_ERROR( info == 0 || ! AbortOnError,
414 "LAPACK gesvd workspace query has thrown an error: %i", ( int ) info );
415
416 if ( info != 0 ) return info ;
417
418 // the optimal size sits in the real part of the first entry.
419 // Keep the reference minimum as floor: a vendor query may
420 // undercut it ( MKL: 7 vs 15 for a 4x3 'A','A' ), and a buffer
421 // sized below the length test above would re-enter this branch
422 // on every call instead of being reused
423 lwork = std::max( lwork, lapack::work_size( Work( 0 ) ) );
424
425 Work.set_size( tRealsPerT * lwork + tRWorkSize );
426 }
427 else
428 {
429 // use the full buffer the caller has provided
430 lwork = ( ( int_t ) Work.length() - tRWorkSize ) / tRealsPerT ;
431 }
432
433 // work segment in the head, rwork in the tail
434 lapack::gesvd( &jobu, &jobvt, &m, &n, A.data(), &lda, S.data(),
435 U.data(), &ldu, VT.data(), &ldvt,
436 reinterpret_cast< T * >( Work.data() ), &lwork,
437 Work.data() + tRealsPerT * lwork, &info );
438
439 // info > 0 : the QR iteration did not converge
440 BELFEM_ERROR( info == 0 || ! AbortOnError,
441 "LAPACK gesvd has thrown an error: %i", ( int ) info );
442
443 return info;
444 }
445
446//------------------------------------------------------------------------------
447}
448#endif //BELFEM_FN_GESVD_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Vector.hpp
belfem::Matrix
Dense column-major matrix.
Definition cl_BZ_Matrix.hpp:28
belfem::Matrix::set_size
void set_size(const size_t aNumRows, const size_t aNumCols)
Definition cl_AR_Matrix.hpp:186
belfem::Matrix::n_rows
size_t n_rows() const
Definition cl_AR_Matrix.hpp:205
belfem::Matrix::n_cols
size_t n_cols() const
Definition cl_AR_Matrix.hpp:213
belfem::Matrix::data
T * data()
Definition cl_AR_Matrix.hpp:135
belfem::Vector
Column vector.
Definition cl_BZ_Vector.hpp:41
lapacktools.hpp
Shared helpers for the LAPACK wrappers, such as leading_dimension().
belfem::lapack
Definition fn_gees.hpp:33
belfem::lapack::dependent_false
constexpr bool dependent_false
Definition lapacktools.hpp:49
belfem::lapack::cplx_float_t
float cplx_float_t
Definition lapacktools.hpp:90
belfem::lapack::gesvd
void gesvd(const char *jobu, const char *jobvt, const int_t *m, const int_t *n, T *a, const int_t *lda, real_t< T > *s, T *u, const int_t *ldu, T *vt, const int_t *ldvt, T *work, const int_t *lwork, real_t< T > *rwork, int_t *info)
Definition fn_gesvd.hpp:136
belfem::lapack::leading_dimension
int_t leading_dimension(const Vector< T > &A)
logical length of a vector operand, as passed to LAPACK as LDB; vector storage is contiguous under bo...
Definition lapacktools.hpp:143
belfem::lapack::dgesvd_
void dgesvd_(char *jobu, char *jobvt, int_t *m, int_t *n, double *a, int_t *lda, double *s, double *u, int_t *ldu, double *vt, int_t *ldvt, double *work, int_t *lwork, int_t *info, fortran_charlen_t lu, fortran_charlen_t lvt)
belfem::lapack::zgesvd_
void zgesvd_(char *jobu, char *jobvt, int_t *m, int_t *n, cplx_double_t *a, int_t *lda, double *s, cplx_double_t *u, int_t *ldu, cplx_double_t *vt, int_t *ldvt, cplx_double_t *work, int_t *lwork, double *rwork, int_t *info, fortran_charlen_t lu, fortran_charlen_t lvt)
belfem::lapack::work_size
int_t work_size(const float &aValue)
lapack reports the optimal work size in the first entry of the work array, which stays real valued fo...
Definition lapacktools.hpp:191
belfem::lapack::fortran_charlen_t
size_t fortran_charlen_t
Definition lapacktools.hpp:100
belfem::lapack::cplx_double_t
double cplx_double_t
Definition lapacktools.hpp:91
belfem::lapack::cgesvd_
void cgesvd_(char *jobu, char *jobvt, int_t *m, int_t *n, cplx_float_t *a, int_t *lda, float *s, cplx_float_t *u, int_t *ldu, cplx_float_t *vt, int_t *ldvt, cplx_float_t *work, int_t *lwork, float *rwork, int_t *info, fortran_charlen_t lu, fortran_charlen_t lvt)
belfem::lapack::real_t
typename real_type< T >::type real_t
Definition lapacktools.hpp:63
belfem::lapack::sgesvd_
void sgesvd_(char *jobu, char *jobvt, int_t *m, int_t *n, float *a, int_t *lda, float *s, float *u, int_t *ldu, float *vt, int_t *ldvt, float *work, int_t *lwork, int_t *info, fortran_charlen_t lu, fortran_charlen_t lvt)
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::gesvd
int_t gesvd(Matrix< T > &A, Vector< lapack::real_t< T > > &S, Matrix< T > &U, Matrix< T > &VT, Vector< lapack::real_t< T > > &Work, const char jobu='A', const char jobvt='A', const bool AbortOnError=true)
singular value decomposition A = U * diag( S ) * VT via LAPACK ?gesvd
Definition fn_gesvd.hpp:340
belfem::value
std::pair< real, unit > value
Definition typedefs.hpp:74
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51