BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_gesv.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,
4 * through 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_GESV_HPP
23#define BELFEM_FN_GESV_HPP
24
25#include "lapacktools.hpp"
26#include "cl_Vector.hpp"
27
28namespace belfem
29{
30 namespace lapack
31 {
32//------------------------------------------------------------------------------
33#ifdef __cplusplus
34 extern "C"
35 {
36#endif
37//------------------------------------------------------------------------------
38
39 void
40 sgesv_( int_t * n,
41 int_t * nrhs,
42 float * a,
43 int_t * lda,
44 int_t * ipiv,
45 float * b,
46 int_t * ldb,
47 int_t * info );
48
49//------------------------------------------------------------------------------
50
51 void
52 dgesv_( int_t * n,
53 int_t * nrhs,
54 double * a,
55 int_t * lda,
56 int_t * ipiv,
57 double * b,
58 int_t * ldb,
59 int_t * info );
60
61//------------------------------------------------------------------------------
62
63 void
64 cgesv_( int_t * n,
65 int_t * nrhs,
66 cplx_float_t * a,
67 int_t * lda,
68 int_t * ipiv,
69 cplx_float_t * b,
70 int_t * ldb,
71 int_t * info );
72
73//------------------------------------------------------------------------------
74
75 void
76 zgesv_( int_t * n,
77 int_t * nrhs,
78 cplx_double_t * a,
79 int_t * lda,
80 int_t * ipiv,
81 cplx_double_t * b,
82 int_t * ldb,
83 int_t * info );
84
85//------------------------------------------------------------------------------
86#ifdef __cplusplus
87 }
88#endif
89//------------------------------------------------------------------------------
90
91 template< typename T >
92 void
93 gesv( const int_t * n,
94 const int_t * nrhs,
95 T * a,
96 const int_t * lda,
97 int_t * ipiv,
98 T * b,
99 const int_t * ldb,
100 int_t * info )
101 {
102 static_assert( dependent_false< T >,
103 "gesv not implemented for selected data type" );
104 }
105
106// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107
108 template<>
109 inline void
110 gesv( const int_t * n,
111 const int_t * nrhs,
112 float * a,
113 const int_t * lda,
114 int_t * ipiv,
115 float * b,
116 const int_t * ldb,
117 int_t * info )
118 {
119 sgesv_(
120 const_cast< int_t * > ( n ),
121 const_cast< int_t * > ( nrhs ),
122 a,
123 const_cast< int_t * > ( lda ),
124 ipiv,
125 b,
126 const_cast< int_t * > ( ldb ),
127 info );
128 }
129
130// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131
132 template<>
133 inline void
134 gesv( const int_t * n,
135 const int_t * nrhs,
136 double * a,
137 const int_t * lda,
138 int_t * ipiv,
139 double * b,
140 const int_t * ldb,
141 int_t * info )
142 {
143 dgesv_(
144 const_cast< int_t * > ( n ),
145 const_cast< int_t * > ( nrhs ),
146 a,
147 const_cast< int_t * > ( lda ),
148 ipiv,
149 b,
150 const_cast< int_t * > ( ldb ),
151 info );
152 }
153
154// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155
156 template<>
157 inline void
158 gesv( const int_t * n,
159 const int_t * nrhs,
160 std::complex< float > * a,
161 const int_t * lda,
162 int_t * ipiv,
163 std::complex< float > * b,
164 const int_t * ldb,
165 int_t * info )
166 {
167 cgesv_(
168 const_cast< int_t * > ( n ),
169 const_cast< int_t * > ( nrhs ),
170 reinterpret_cast< cplx_float_t * >( a ),
171 const_cast< int_t * > ( lda ),
172 ipiv,
173 reinterpret_cast< cplx_float_t * >( b ),
174 const_cast< int_t * > ( ldb ),
175 info );
176 }
177
178// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179
180 template<>
181 inline void
182 gesv( const int_t * n,
183 const int_t * nrhs,
184 std::complex< double > * a,
185 const int_t * lda,
186 int_t * ipiv,
187 std::complex< double > * b,
188 const int_t * ldb,
189 int_t * info )
190 {
191 zgesv_(
192 const_cast< int_t * > ( n ),
193 const_cast< int_t * > ( nrhs ),
194 reinterpret_cast< cplx_double_t * >( a ),
195 const_cast< int_t * > ( lda ),
196 ipiv,
197 reinterpret_cast< cplx_double_t * >( b ),
198 const_cast< int_t * > ( ldb ) ,
199 info );
200 }
201
202//------------------------------------------------------------------------------
203 }
204//------------------------------------------------------------------------------
205
220 template< typename T >
221 int_t
222 gesv( Matrix< T > & A, Vector< T > & B, Vector< int_t > & Pivot,
223 const bool AbortOnError = true )
224 {
225 BELFEM_ASSERT( A.n_rows() == A.n_cols(),
226 "Matrix A must be square ( is %lu x %lu )",
227 ( long unsigned int ) A.n_rows(),
228 ( long unsigned int ) A.n_cols() );
229
230 BELFEM_ASSERT( B.length() == A.n_rows(),
231 "Length of right hand side does not match ( %lu vs %lu )",
232 ( long unsigned int ) B.length(),
233 ( long unsigned int ) A.n_rows() );
234
235 BELFEM_ASSERT( Pivot.length() >= A.n_rows(),
236 "Pivot vector is too short ( %lu, need %lu )",
237 ( long unsigned int ) Pivot.length(),
238 ( long unsigned int ) A.n_rows() );
239
240 // size of matrix
241 int_t n = ( int_t ) A.n_rows();
242
243 int_t nrhs = 1;
244
245 int_t lda = lapack::leading_dimension( A );
246 int_t ldb = lapack::leading_dimension( B );
247
248 // error code
249 int_t info = 0;
250
251 // call lapack
252 lapack::gesv( &n,
253 &nrhs,
254 A.data(),
255 &lda,
256 Pivot.data(),
257 B.data(),
258 &ldb,
259 &info );
260
261 BELFEM_ERROR( info == 0 || ! AbortOnError,
262 "LAPACK gesv has thrown an error: %i", ( int ) info );
263
264 return info ;
265 }
266
267//------------------------------------------------------------------------------
268
281 template< typename T >
282 int_t
283 gesv( Matrix< T > & A, Matrix< T > & B, Vector< int_t > & Pivot, const bool AbortOnError = true )
284 {
285 BELFEM_ASSERT( A.n_rows() == A.n_cols(),
286 "Matrix A must be square ( is %lu x %lu )",
287 ( long unsigned int ) A.n_rows(),
288 ( long unsigned int ) A.n_cols() );
289
290 BELFEM_ASSERT( B.n_rows() == A.n_rows(),
291 "Number of rows of right hand side does not match ( %lu vs %lu )",
292 ( long unsigned int ) B.n_rows(),
293 ( long unsigned int ) A.n_rows() );
294
295 BELFEM_ASSERT( Pivot.length() >= A.n_rows(),
296 "Pivot vector is too short ( %lu, need %lu )",
297 ( long unsigned int ) Pivot.length(),
298 ( long unsigned int ) A.n_rows() );
299
300 // size of matrix
301 int_t n = ( int_t ) A.n_rows();
302
303 int_t nrhs = ( int_t ) B.n_cols();
304
305 int_t lda = lapack::leading_dimension( A );
306 int_t ldb = lapack::leading_dimension( B );
307
308 // error code
309 int_t info = 0;
310
311 // call lapack
312 lapack::gesv( &n,
313 & nrhs,
314 A.data(),
315 & lda,
316 Pivot.data(),
317 B.data(),
318 & ldb,
319 & info );
320
321 BELFEM_ERROR( info == 0 || ! AbortOnError,
322 "LAPACK gesv has thrown an error: %i", ( int ) info );
323
324 return info ;
325 }
326
327//------------------------------------------------------------------------------
328}
329#endif //BELFEM_FN_GESV_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::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::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::dgesv_
void dgesv_(int_t *n, int_t *nrhs, double *a, int_t *lda, int_t *ipiv, double *b, int_t *ldb, int_t *info)
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::sgesv_
void sgesv_(int_t *n, int_t *nrhs, float *a, int_t *lda, int_t *ipiv, float *b, int_t *ldb, int_t *info)
belfem::lapack::cgesv_
void cgesv_(int_t *n, int_t *nrhs, cplx_float_t *a, int_t *lda, int_t *ipiv, cplx_float_t *b, int_t *ldb, int_t *info)
belfem::lapack::gesv
void gesv(const int_t *n, const int_t *nrhs, T *a, const int_t *lda, int_t *ipiv, T *b, const int_t *ldb, int_t *info)
Definition fn_gesv.hpp:93
belfem::lapack::cplx_double_t
double cplx_double_t
Definition lapacktools.hpp:91
belfem::lapack::zgesv_
void zgesv_(int_t *n, int_t *nrhs, cplx_double_t *a, int_t *lda, int_t *ipiv, cplx_double_t *b, int_t *ldb, int_t *info)
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::gesv
int_t gesv(Matrix< T > &A, Vector< T > &B, Vector< int_t > &Pivot, const bool AbortOnError=true)
solve the square linear system A * x = b via LAPACK ?gesv ( LU factorization with partial pivoting )
Definition fn_gesv.hpp:222
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51