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