BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_gels.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_GELS_HPP
23#define BELFEM_FN_GELS_HPP
24
25#include <algorithm>
26#include <complex>
27
28#include "assert.hpp"
29#include "cl_Vector.hpp"
30#include "cl_Matrix.hpp"
31#include "lapacktools.hpp"
32
33//------------------------------------------------------------------------------
34namespace belfem
35{
36 namespace lapack
37 {
38//------------------------------------------------------------------------------
39#ifdef __cplusplus
40 extern "C"
41 {
42#endif
43//------------------------------------------------------------------------------
44
45 void
46 sgels_( char * trans,
47 int_t * m,
48 int_t * n,
49 int_t * nrhs,
50 float * a,
51 int_t * lda,
52 float * b,
53 int_t * ldb,
54 float * work,
55 int_t * lwork,
56 int_t * info,
57 fortran_charlen_t lt );
58
59//------------------------------------------------------------------------------
60
61 void
62 dgels_( char * trans,
63 int_t * m,
64 int_t * n,
65 int_t * nrhs,
66 double * a,
67 int_t * lda,
68 double * b,
69 int_t * ldb,
70 double * work,
71 int_t * lwork,
72 int_t * info,
73 fortran_charlen_t lt );
74
75//------------------------------------------------------------------------------
76
77 void
78 cgels_( char * trans,
79 int_t * m,
80 int_t * n,
81 int_t * nrhs,
82 cplx_float_t * a,
83 int_t * lda,
84 cplx_float_t * b,
85 int_t * ldb,
86 cplx_float_t * work,
87 int_t * lwork,
88 int_t * info,
89 fortran_charlen_t lt );
90
91//------------------------------------------------------------------------------
92
93 void
94 zgels_( char * trans,
95 int_t * m,
96 int_t * n,
97 int_t * nrhs,
98 cplx_double_t * a,
99 int_t * lda,
100 cplx_double_t * b,
101 int_t * ldb,
102 cplx_double_t * work,
103 int_t * lwork,
104 int_t * info,
105 fortran_charlen_t lt );
106
107//------------------------------------------------------------------------------
108
109#ifdef __cplusplus
110 }
111#endif
112
113//------------------------------------------------------------------------------
114
115 template< typename T >
116 void
117 gels( const char * trans,
118 const int_t * m,
119 const int_t * n,
120 const int_t * nrhs,
121 T * a,
122 const int_t * lda,
123 T * b,
124 const int_t * ldb,
125 T * work,
126 const int_t * lwork,
127 int_t * info )
128 {
129 static_assert( dependent_false< T >,
130 "gels not implemented for selected data type" );
131 }
132
133// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134
135 template<>
136 inline void
137 gels( const char * trans,
138 const int_t * m,
139 const int_t * n,
140 const int_t * nrhs,
141 float * a,
142 const int_t * lda,
143 float * b,
144 const int_t * ldb,
145 float * work,
146 const int_t * lwork,
147 int_t * info )
148 {
149 sgels_(
150 const_cast< char * >( trans ),
151 const_cast< int_t * >( m ),
152 const_cast< int_t * >( n ),
153 const_cast< int_t * >( nrhs ),
154 a,
155 const_cast< int_t * >( lda ),
156 b,
157 const_cast< int_t * >( ldb ),
158 work,
159 const_cast< int_t * >( lwork ),
160 info,
161 1 );
162 }
163
164// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165
166 template<>
167 inline void
168 gels( const char * trans,
169 const int_t * m,
170 const int_t * n,
171 const int_t * nrhs,
172 double * a,
173 const int_t * lda,
174 double * b,
175 const int_t * ldb,
176 double * work,
177 const int_t * lwork,
178 int_t * info )
179 {
180 dgels_(
181 const_cast< char * >( trans ),
182 const_cast< int_t * >( m ),
183 const_cast< int_t * >( n ),
184 const_cast< int_t * >( nrhs ),
185 a,
186 const_cast< int_t * >( lda ),
187 b,
188 const_cast< int_t * >( ldb ),
189 work,
190 const_cast< int_t * >( lwork ),
191 info,
192 1 );
193 }
194
195// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196
197 template<>
198 inline void
199 gels( const char * trans,
200 const int_t * m,
201 const int_t * n,
202 const int_t * nrhs,
203 std::complex< float > * a,
204 const int_t * lda,
205 std::complex< float > * b,
206 const int_t * ldb,
207 std::complex< float > * work,
208 const int_t * lwork,
209 int_t * info )
210 {
211 cgels_(
212 const_cast< char * >( trans ),
213 const_cast< int_t * >( m ),
214 const_cast< int_t * >( n ),
215 const_cast< int_t * >( nrhs ),
216 reinterpret_cast< cplx_float_t * > ( a ),
217 const_cast< int_t * >( lda ),
218 reinterpret_cast< cplx_float_t * > ( b ),
219 const_cast< int_t * >( ldb ),
220 reinterpret_cast< cplx_float_t * > ( work ),
221 const_cast< int_t * >( lwork ),
222 info,
223 1 );
224 }
225
226// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227
228 template<>
229 inline void
230 gels( const char * trans,
231 const int_t * m,
232 const int_t * n,
233 const int_t * nrhs,
234 std::complex< double > * a,
235 const int_t * lda,
236 std::complex< double > * b,
237 const int_t * ldb,
238 std::complex< double > * work,
239 const int_t * lwork,
240 int_t * info )
241 {
242 zgels_(
243 const_cast< char * >( trans ),
244 const_cast< int_t * >( m ),
245 const_cast< int_t * >( n ),
246 const_cast< int_t * >( nrhs ),
247 reinterpret_cast< cplx_double_t * > ( a ),
248 const_cast< int_t * >( lda ),
249 reinterpret_cast< cplx_double_t * > ( b ),
250 const_cast< int_t * >( ldb ),
251 reinterpret_cast< cplx_double_t * > ( work ),
252 const_cast< int_t * >( lwork ),
253 info,
254 1 );
255 }
256
257//------------------------------------------------------------------------------
258 }
259//------------------------------------------------------------------------------
260
276 template< typename T >
277 int_t
278 gels( Matrix< T > & A, Vector< T > & B, Vector< T > & Work, const bool AbortOnError = true )
279 {
280 int_t m = A.n_rows();
281 int_t n = A.n_cols();
282
283 // B carries the right hand side ( m ) as well as the solution ( n )
284 BELFEM_ASSERT( static_cast< int_t >( B.length() ) >= std::max( m, n ),
285 "Length of rhs vector is %i, but must be at least %i.",
286 ( int ) B.length(), ( int ) std::max( m, n ) );
287
288 char trans = 'N';
289
290 int_t nrhs = 1 ;
291 int_t lda = lapack::leading_dimension( A );
292
293 // leading dimension of B, taken from the actual allocation
294 int_t ldb = lapack::leading_dimension( B );
295
296 int_t mn = std::min( m, n );
297 int_t lwork = std::max< int_t >( 1, mn + std::max( mn, nrhs ) );
298 int_t info = 0 ;
299
300 // check length of work array
301 if ( static_cast< int_t >( Work.length() ) < lwork )
302 {
303 // ask lapack for the optimal size, this leaves A and B untouched
304 Work.set_size( 1 );
305
306 int_t query = -1 ;
307
308 lapack::gels( &trans, &m, &n, &nrhs, A.data(), &lda, B.data(), &ldb, Work.data(), &query, &info );
309
310 BELFEM_ERROR( info == 0 || ! AbortOnError,
311 "LAPACK gels workspace query has thrown an error: %i", ( int ) info );
312
313 if ( info != 0 ) return info ;
314
315 lwork = lapack::work_size( Work( 0 ) );
316
317 Work.set_size( lwork );
318 }
319 else
320 {
321 // use the full buffer the caller has provided
322 lwork = ( int_t ) Work.length();
323 }
324
325 lapack::gels(
326 &trans,
327 &m,
328 &n,
329 &nrhs,
330 A.data(),
331 &lda,
332 B.data(),
333 &ldb,
334 Work.data(),
335 &lwork, &info );
336
337 // info > 0 : a diagonal entry of the triangular factor is zero, so A
338 // does not have full rank and no solution was computed
339 BELFEM_ERROR( info == 0 || ! AbortOnError,
340 "LAPACK gels has thrown an error: %i", ( int ) info );
341
342 return info;
343 }
344
345//------------------------------------------------------------------------------
346
359 template< typename T >
360 int_t
361 gels( Matrix< T > & A, Matrix< T > & B, Vector< T > & Work, const bool AbortOnError = true )
362 {
363 int_t m = A.n_rows();
364 int_t n = A.n_cols();
365
366 // B carries the right hand sides ( m ) as well as the solutions ( n )
367 BELFEM_ASSERT( static_cast< int_t >( B.n_rows() ) >= std::max( m, n ),
368 "Number of rows of rhs matrix is %i, but must be at least %i.",
369 ( int ) B.n_rows(), ( int ) std::max( m, n ) );
370
371 char trans = 'N';
372
373 int_t nrhs = B.n_cols() ;
374
375 int_t lda = lapack::leading_dimension( A );
376 int_t ldb = lapack::leading_dimension( B );
377
378 int_t mn = std::min( m, n );
379 int_t lwork = std::max< int_t >( 1, mn + std::max( mn, nrhs ) );
380 int_t info = 0 ;
381
382 // check length of work array
383 if ( static_cast< int_t >( Work.length() ) < lwork )
384 {
385 // ask lapack for the optimal size, this leaves A and B untouched
386 Work.set_size( 1 );
387
388 int_t query = -1 ;
389
390 lapack::gels( &trans, &m, &n, &nrhs, A.data(), &lda, B.data(), &ldb, Work.data(), &query, &info );
391
392 BELFEM_ERROR( info == 0 || ! AbortOnError,
393 "LAPACK gels workspace query has thrown an error: %i", ( int ) info );
394
395 if ( info != 0 ) return info ;
396
397 lwork = lapack::work_size( Work( 0 ) );
398
399 Work.set_size( lwork );
400 }
401 else
402 {
403 // use the full buffer the caller has provided
404 lwork = ( int_t ) Work.length();
405 }
406
407 lapack::gels( &trans, &m, &n, &nrhs, A.data(), &lda, B.data(), &ldb, Work.data(), &lwork, &info );
408
409 // info > 0 : a diagonal entry of the triangular factor is zero, so A
410 // does not have full rank and no solution was computed
411 BELFEM_ERROR( info == 0 || ! AbortOnError,
412 "LAPACK gels has thrown an error: %i", ( int ) info );
413
414 return info;
415 }
416
417//------------------------------------------------------------------------------
418}
419#endif //BELFEM_FN_GELS_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Matrix.hpp
cl_Vector.hpp
belfem::Matrix
Dense column-major matrix.
Definition cl_BZ_Matrix.hpp:28
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
belfem::Vector::data
T * data()
expose the underlying raw pointer ( writable version )
Definition cl_AR_Vector.hpp:182
belfem::Vector::set_size
void set_size(const size_t aNumRows)
change the size of the vector
Definition cl_AR_Vector.hpp:237
belfem::Vector::length
size_t length() const
get the length of the vector
Definition cl_AR_Vector.hpp:257
lapacktools.hpp
Shared helpers for the LAPACK wrappers, such as leading_dimension().
belfem::lapack
Definition fn_gees.hpp:33
belfem::lapack::dgels_
void dgels_(char *trans, int_t *m, int_t *n, int_t *nrhs, double *a, int_t *lda, double *b, int_t *ldb, double *work, int_t *lwork, int_t *info, fortran_charlen_t lt)
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::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::sgels_
void sgels_(char *trans, int_t *m, int_t *n, int_t *nrhs, float *a, int_t *lda, float *b, int_t *ldb, float *work, int_t *lwork, int_t *info, fortran_charlen_t lt)
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::cgels_
void cgels_(char *trans, int_t *m, int_t *n, int_t *nrhs, cplx_float_t *a, int_t *lda, cplx_float_t *b, int_t *ldb, cplx_float_t *work, int_t *lwork, int_t *info, fortran_charlen_t lt)
belfem::lapack::fortran_charlen_t
size_t fortran_charlen_t
Definition lapacktools.hpp:100
belfem::lapack::zgels_
void zgels_(char *trans, int_t *m, int_t *n, int_t *nrhs, cplx_double_t *a, int_t *lda, cplx_double_t *b, int_t *ldb, cplx_double_t *work, int_t *lwork, int_t *info, fortran_charlen_t lt)
belfem::lapack::cplx_double_t
double cplx_double_t
Definition lapacktools.hpp:91
belfem::lapack::gels
void gels(const char *trans, const int_t *m, const int_t *n, const int_t *nrhs, T *a, const int_t *lda, T *b, const int_t *ldb, T *work, const int_t *lwork, int_t *info)
Definition fn_gels.hpp:117
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::gels
int_t gels(Matrix< T > &A, Vector< T > &B, Vector< T > &Work, const bool AbortOnError=true)
solve the least squares problem min || A * x - b || via LAPACK ?gels ( QR or LQ factorization,...
Definition fn_gels.hpp:278
belfem::trans
auto trans(Matrix< T > &aMatrix) -> decltype(trans(aMatrix.matrix_data()))
Definition fn_trans.hpp:74
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51