BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_getri.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_GETRI_HPP
23#define BELFEM_FN_GETRI_HPP
24
25#include "assert.hpp"
26#include "lapacktools.hpp"
27
28//------------------------------------------------------------------------------
29namespace belfem
30{
31 namespace lapack
32 {
33//------------------------------------------------------------------------------
34#ifdef __cplusplus
35 extern "C"
36 {
37#endif
38//------------------------------------------------------------------------------
39
40 void
41 sgetri_( int_t * n,
42 float * a,
43 int_t * lda,
44 int_t * ipiv,
45 float * work,
46 int_t * lwork,
47 int_t * info );
48
49//------------------------------------------------------------------------------
50
51 void
52 dgetri_( int_t * n,
53 double * a,
54 int_t * lda,
55 int_t * ipiv,
56 double * work,
57 int_t * lwork,
58 int_t * info );
59//------------------------------------------------------------------------------
60
61 void
62 cgetri_( int_t * n,
63 cplx_float_t * a,
64 int_t * lda,
65 int_t * ipiv,
66 cplx_float_t * work,
67 int_t * lwork,
68 int_t * info );
69
70//------------------------------------------------------------------------------
71
72 void
73 zgetri_( int_t * n,
74 cplx_double_t * a,
75 int_t * lda,
76 int_t * ipiv,
77 cplx_double_t * work,
78 int_t * lwork,
79 int_t * info );
80
81
82//------------------------------------------------------------------------------
83
84#ifdef __cplusplus
85}
86#endif
87
88//------------------------------------------------------------------------------
89
90 template< typename T >
91 void
92 getri( const int_t * n,
93 T * a,
94 const int_t * lda,
95 int_t * ipiv,
96 T * work,
97 const int_t * lwork,
98 int_t * info )
99 {
100 static_assert( dependent_false< T >,
101 "getri not implemented for selected data type" );
102 }
103
104// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105
106 template<>
107 inline void
108 getri( const int_t * n,
109 float * a,
110 const int_t * lda,
111 int_t * ipiv,
112 float * work,
113 const int_t * lwork,
114 int_t * info )
115 {
116 sgetri_( const_cast< int_t * >( n ),
117 a,
118 const_cast< int_t * >( lda ),
119 ipiv,
120 work,
121 const_cast< int_t * >( lwork ),
122 info );
123 }
124
125// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126
127 template<>
128 inline void
129 getri( const int_t * n,
130 double * a,
131 const int_t * lda,
132 int_t * ipiv,
133 double * work,
134 const int_t * lwork,
135 int_t * info )
136 {
137 dgetri_( const_cast< int_t * >( n ),
138 a,
139 const_cast< int_t * >( lda ),
140 ipiv,
141 work,
142 const_cast< int_t * >( lwork ),
143 info );
144 }
145
146// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147
148 template<>
149 inline void
150 getri( const int_t * n,
151 std::complex< float > * a,
152 const int_t * lda,
153 int_t * ipiv,
154 std::complex< float > * work,
155 const int_t * lwork,
156 int_t * info )
157 {
158 cgetri_( const_cast< int_t * >( n ),
159 reinterpret_cast< cplx_float_t * >( a ),
160 const_cast< int_t * >( lda ),
161 ipiv,
162 reinterpret_cast< cplx_float_t * >( work ),
163 const_cast< int_t * >( lwork ),
164 info );
165 }
166
167// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168
169 template<>
170 inline void
171 getri( const int_t * n,
172 std::complex< double > * a,
173 const int_t * lda,
174 int_t * ipiv,
175 std::complex< double > * work,
176 const int_t * lwork,
177 int_t * info )
178 {
179 zgetri_( const_cast< int_t * >( n ),
180 reinterpret_cast< cplx_double_t * >( a ),
181 const_cast< int_t * >( lda ),
182 ipiv,
183 reinterpret_cast< cplx_double_t * >( work ),
184 const_cast< int_t * >( lwork ),
185 info );
186 }
187
188//------------------------------------------------------------------------------
189 } /* end namespace lapack */
190//------------------------------------------------------------------------------
191
205 template< typename T >
206 int_t
207 getri( Matrix< T > & A , Vector< int_t > & Pivot, Vector< T > & Work, const bool AbortOnError = true )
208 {
209 BELFEM_ASSERT( A.n_rows() == A.n_cols(),
210 "Matrix A must be square ( is %lu x %lu )",
211 ( long unsigned int ) A.n_rows(),
212 ( long unsigned int ) A.n_cols() );
213
214 BELFEM_ASSERT( Pivot.length() >= A.n_rows(),
215 "Pivot vector is too short ( %lu, need %lu )",
216 ( long unsigned int ) Pivot.length(),
217 ( long unsigned int ) A.n_rows() );
218
219 int_t n = ( int_t ) A.n_cols() ;
220 int_t lda = lapack::leading_dimension( A );
221 int_t info = 0 ;
222
223 int_t lwork = std::max< int_t >( 1, n );
224
225 if ( static_cast< int_t >( Work.length() ) < lwork )
226 {
227 // ask lapack for the optimal size, this leaves A untouched
228 Work.set_size( 1 );
229
230 int_t query = -1 ;
231
232 lapack::getri( &n, A.data(), &lda, Pivot.data(), Work.data(), &query, &info );
233
234 BELFEM_ERROR( info == 0 || ! AbortOnError,
235 "LAPACK getri workspace query has thrown an error: %i", ( int ) info );
236
237 if ( info != 0 ) return info ;
238
239 lwork = lapack::work_size( Work( 0 ) );
240
241 Work.set_size( lwork );
242 }
243 else
244 {
245 // use the full buffer the caller has provided
246 lwork = ( int_t ) Work.length();
247 }
248
249 lapack::getri( &n, A.data(), &lda, Pivot.data(), Work.data(), &lwork, &info );
250
251 BELFEM_ERROR( info == 0 || ! AbortOnError,
252 "LAPACK getri has thrown an error: %i", ( int ) info );
253
254 return info ;
255 }
256} /* end namespace belfem */
257#endif //BELFEM_FN_GETRI_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::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::sgetri_
void sgetri_(int_t *n, float *a, int_t *lda, int_t *ipiv, float *work, int_t *lwork, int_t *info)
belfem::lapack::zgetri_
void zgetri_(int_t *n, cplx_double_t *a, int_t *lda, int_t *ipiv, cplx_double_t *work, int_t *lwork, int_t *info)
belfem::lapack::dependent_false
constexpr bool dependent_false
Definition lapacktools.hpp:49
belfem::lapack::getri
void getri(const int_t *n, T *a, const int_t *lda, int_t *ipiv, T *work, const int_t *lwork, int_t *info)
Definition fn_getri.hpp:92
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::cgetri_
void cgetri_(int_t *n, cplx_float_t *a, int_t *lda, int_t *ipiv, cplx_float_t *work, int_t *lwork, int_t *info)
belfem::lapack::dgetri_
void dgetri_(int_t *n, double *a, int_t *lda, int_t *ipiv, double *work, int_t *lwork, int_t *info)
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::cplx_double_t
double cplx_double_t
Definition lapacktools.hpp:91
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::getri
int_t getri(Matrix< T > &A, Vector< int_t > &Pivot, Vector< T > &Work, const bool AbortOnError=true)
invert a square matrix in place via LAPACK ?getri, using the LU factorization computed by getrf()
Definition fn_getri.hpp:207
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51