BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_invert_symmetric.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
12#ifndef BELFEM_INVERT_SYMMETRIC_HPP
13#define BELFEM_INVERT_SYMMETRIC_HPP
14
15#include "typedefs.hpp"
16#include "assert.hpp"
17#include "fn_getrf.hpp"
18#include "fn_getri.hpp"
19#include "fn_gemm.hpp"
20
21
22#include "cl_Tensor.hpp"
23#include "cl_Vector.hpp"
24
25
26
27namespace belfem
28{
29//----------------------------------------------------------------------------
30
31 namespace tensor
32 {
33 template < typename T >
34 void
35 invert_symmetric(
36 Tensor< T > & aTensor,
37 Vector< T > & aWork,
38 Vector< int_t > & aPivot )
39 {
40 BELFEM_ASSERT( aTensor.is_3333(), "The tensor must be a 3x3x3x3 tensor" );
41
42 BELFEM_ASSERT( aWork.length() >= 72,
43 "The work vector must have a length of at least 72." );
44
45 BELFEM_ASSERT( aPivot.length() >= 36,
46 "The pivot vector must have a length of at least 36" );
47
48 // pointer to tensor data
49 T * tData = aTensor.data();
50
51 // pointer to matrix data
52 T * tWorkA = aWork.data();
53
54 // pointer to Indentity data
55 T * tWorkB = aWork.data() + 36 ;
56
57 // pointer to Indentity data, can use tensor here
58 T * tWorkI = aTensor.data();
59
60 // copy data into WorkA vector
61 tWorkA[ 0 ] = tData[ 0 ];
62 tWorkA[ 1 ] = tData[ 1 ];
63 tWorkA[ 2 ] = tData[ 2 ];
64 tWorkA[ 3 ] = tData[ 4 ];
65 tWorkA[ 4 ] = tData[ 5 ];
66 tWorkA[ 5 ] = tData[ 8 ];
67 tWorkA[ 6 ] = tData[ 9 ] + tData[ 9 ];
68 tWorkA[ 7 ] = tData[ 10 ] + tData[ 10 ];
69 tWorkA[ 8 ] = tData[ 11 ] + tData[ 11 ];
70 tWorkA[ 9 ] = tData[ 13 ] + tData[ 13 ];
71 tWorkA[ 10 ] = tData[ 14 ] + tData[ 14 ];
72 tWorkA[ 11 ] = tData[ 17 ] + tData[ 17 ];
73 tWorkA[ 12 ] = tData[ 18 ] + tData[ 18 ];
74 tWorkA[ 13 ] = tData[ 19 ] + tData[ 19 ];
75 tWorkA[ 14 ] = tData[ 20 ] + tData[ 20 ];
76 tWorkA[ 15 ] = tData[ 22 ] + tData[ 22 ];
77 tWorkA[ 16 ] = tData[ 23 ] + tData[ 23 ];
78 tWorkA[ 17 ] = tData[ 26 ] + tData[ 26 ];
79 tWorkA[ 18 ] = tData[ 36 ];
80 tWorkA[ 19 ] = tData[ 37 ];
81 tWorkA[ 20 ] = tData[ 38 ];
82 tWorkA[ 21 ] = tData[ 40 ];
83 tWorkA[ 22 ] = tData[ 41 ];
84 tWorkA[ 23 ] = tData[ 44 ];
85 tWorkA[ 24 ] = tData[ 45 ] + tData[ 45 ];
86 tWorkA[ 25 ] = tData[ 46 ] + tData[ 46 ];
87 tWorkA[ 26 ] = tData[ 47 ] + tData[ 47 ];
88 tWorkA[ 27 ] = tData[ 49 ] + tData[ 49 ];
89 tWorkA[ 28 ] = tData[ 50 ] + tData[ 50 ];
90 tWorkA[ 29 ] = tData[ 53 ] + tData[ 53 ];
91 tWorkA[ 30 ] = tData[ 72 ];
92 tWorkA[ 31 ] = tData[ 73 ];
93 tWorkA[ 32 ] = tData[ 74 ];
94 tWorkA[ 33 ] = tData[ 76 ];
95 tWorkA[ 34 ] = tData[ 77 ];
96 tWorkA[ 35 ] = tData[ 80 ];
97
98
99 // call lapack for matrix-matrix multiplication
100 int_t info = 0 ;
101 char trans = 'N';
102 int_t n = 6 ;
103 int_t lwork = 36 ;
104 T alpha = 1.0 ;
105 T beta = 0.0 ;
106 // create an LU factorization for A
107 lapack::getrf( &n, &n, tWorkA, &n, aPivot.data(), &info );
108
109 BELFEM_ERROR( info == 0, "LAPACK getrf has thrown error %i", ( int ) info );
110
111 // create the inverse of A using the LU factorization
112 lapack::getri( &n, tWorkA, &n, aPivot.data(), tWorkB, &lwork, &info );
113
114 BELFEM_ERROR( info == 0, "LAPACK getri has thrown error %i", ( int ) info );
115
116 // initialize identity matrix
117 tWorkI[ 0 ] = 1.0 ;
118 tWorkI[ 1 ] = 0.0 ;
119 tWorkI[ 2 ] = 0.0 ;
120 tWorkI[ 3 ] = 0.0 ;
121 tWorkI[ 4 ] = 0.0 ;
122 tWorkI[ 5 ] = 0.0 ;
123
124 tWorkI[ 6 ] = 0.0 ;
125 tWorkI[ 7 ] = 0.5 ;
126 tWorkI[ 8 ] = 0.0 ;
127 tWorkI[ 9 ] = 0.0 ;
128 tWorkI[ 10 ] = 0.0 ;
129 tWorkI[ 11 ] = 0.0 ;
130
131 tWorkI[ 12 ] = 0.0 ;
132 tWorkI[ 13 ] = 0.0 ;
133 tWorkI[ 14 ] = 0.5 ;
134 tWorkI[ 15 ] = 0.0 ;
135 tWorkI[ 16 ] = 0.0 ;
136 tWorkI[ 17 ] = 0.0 ;
137
138 tWorkI[ 18 ] = 0.0 ;
139 tWorkI[ 19 ] = 0.0 ;
140 tWorkI[ 20 ] = 0.0 ;
141 tWorkI[ 21 ] = 1.0 ;
142 tWorkI[ 22 ] = 0.0 ;
143 tWorkI[ 23 ] = 0.0 ;
144
145 tWorkI[ 24 ] = 0.0 ;
146 tWorkI[ 25 ] = 0.0 ;
147 tWorkI[ 26 ] = 0.0 ;
148 tWorkI[ 27 ] = 0.0 ;
149 tWorkI[ 28 ] = 0.5 ;
150 tWorkI[ 29 ] = 0.0 ;
151
152 tWorkI[ 30 ] = 0.0 ;
153 tWorkI[ 31 ] = 0.0 ;
154 tWorkI[ 32 ] = 0.0 ;
155 tWorkI[ 33 ] = 0.0 ;
156 tWorkI[ 34 ] = 0.0 ;
157 tWorkI[ 35 ] = 1.0 ;
158
159
160 lapack::gemm( &trans, &trans,
161 &n, &n, &n,
162 &alpha,
163 tWorkA,
164 &n,
165 tWorkI,
166 &n,
167 &beta,
168 tWorkB,
169 &n );
170
171 // write data back
172 tData[ 0 ] = tWorkB[ 0 ];
173 tData[ 1 ] = tWorkB[ 1 ];
174 tData[ 2 ] = tWorkB[ 2 ];
175 tData[ 3 ] = tWorkB[ 1 ];
176 tData[ 4 ] = tWorkB[ 3 ];
177 tData[ 5 ] = tWorkB[ 4 ];
178 tData[ 6 ] = tWorkB[ 2 ];
179 tData[ 7 ] = tWorkB[ 4 ];
180 tData[ 8 ] = tWorkB[ 5 ];
181 tData[ 9 ] = tWorkB[ 6 ];
182 tData[ 10 ] = tWorkB[ 7 ];
183 tData[ 11 ] = tWorkB[ 8 ];
184 tData[ 12 ] = tWorkB[ 7 ];
185 tData[ 13 ] = tWorkB[ 9 ];
186 tData[ 14 ] = tWorkB[ 10 ];
187 tData[ 15 ] = tWorkB[ 8 ];
188 tData[ 16 ] = tWorkB[ 10 ];
189 tData[ 17 ] = tWorkB[ 11 ];
190 tData[ 18 ] = tWorkB[ 12 ];
191 tData[ 19 ] = tWorkB[ 13 ];
192 tData[ 20 ] = tWorkB[ 14 ];
193 tData[ 21 ] = tWorkB[ 13 ];
194 tData[ 22 ] = tWorkB[ 15 ];
195 tData[ 23 ] = tWorkB[ 16 ];
196 tData[ 24 ] = tWorkB[ 14 ];
197 tData[ 25 ] = tWorkB[ 16 ];
198 tData[ 26 ] = tWorkB[ 17 ];
199 tData[ 27 ] = tWorkB[ 6 ];
200 tData[ 28 ] = tWorkB[ 7 ];
201 tData[ 29 ] = tWorkB[ 8 ];
202 tData[ 30 ] = tWorkB[ 7 ];
203 tData[ 31 ] = tWorkB[ 9 ];
204 tData[ 32 ] = tWorkB[ 10 ];
205 tData[ 33 ] = tWorkB[ 8 ];
206 tData[ 34 ] = tWorkB[ 10 ];
207 tData[ 35 ] = tWorkB[ 11 ];
208 tData[ 36 ] = tWorkB[ 18 ];
209 tData[ 37 ] = tWorkB[ 19 ];
210 tData[ 38 ] = tWorkB[ 20 ];
211 tData[ 39 ] = tWorkB[ 19 ];
212 tData[ 40 ] = tWorkB[ 21 ];
213 tData[ 41 ] = tWorkB[ 22 ];
214 tData[ 42 ] = tWorkB[ 20 ];
215 tData[ 43 ] = tWorkB[ 22 ];
216 tData[ 44 ] = tWorkB[ 23 ];
217 tData[ 45 ] = tWorkB[ 24 ];
218 tData[ 46 ] = tWorkB[ 25 ];
219 tData[ 47 ] = tWorkB[ 26 ];
220 tData[ 48 ] = tWorkB[ 25 ];
221 tData[ 49 ] = tWorkB[ 27 ];
222 tData[ 50 ] = tWorkB[ 28 ];
223 tData[ 51 ] = tWorkB[ 26 ];
224 tData[ 52 ] = tWorkB[ 28 ];
225 tData[ 53 ] = tWorkB[ 29 ];
226 tData[ 54 ] = tWorkB[ 12 ];
227 tData[ 55 ] = tWorkB[ 13 ];
228 tData[ 56 ] = tWorkB[ 14 ];
229 tData[ 57 ] = tWorkB[ 13 ];
230 tData[ 58 ] = tWorkB[ 15 ];
231 tData[ 59 ] = tWorkB[ 16 ];
232 tData[ 60 ] = tWorkB[ 14 ];
233 tData[ 61 ] = tWorkB[ 16 ];
234 tData[ 62 ] = tWorkB[ 17 ];
235 tData[ 63 ] = tWorkB[ 24 ];
236 tData[ 64 ] = tWorkB[ 25 ];
237 tData[ 65 ] = tWorkB[ 26 ];
238 tData[ 66 ] = tWorkB[ 25 ];
239 tData[ 67 ] = tWorkB[ 27 ];
240 tData[ 68 ] = tWorkB[ 28 ];
241 tData[ 69 ] = tWorkB[ 26 ];
242 tData[ 70 ] = tWorkB[ 28 ];
243 tData[ 71 ] = tWorkB[ 29 ];
244 tData[ 72 ] = tWorkB[ 30 ];
245 tData[ 73 ] = tWorkB[ 31 ];
246 tData[ 74 ] = tWorkB[ 32 ];
247 tData[ 75 ] = tWorkB[ 31 ];
248 tData[ 76 ] = tWorkB[ 33 ];
249 tData[ 77 ] = tWorkB[ 34 ];
250 tData[ 78 ] = tWorkB[ 32 ];
251 tData[ 79 ] = tWorkB[ 34 ];
252 tData[ 80 ] = tWorkB[ 35 ];
253
254 }
255 }
256//----------------------------------------------------------------------------
257
258}
259
260#endif //BELFEM_INVERT_SYMMETRIC_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Tensor.hpp
cl_Vector.hpp
belfem::Tensor
Third- or fourth-order tensor container; the constitutive helpers (contraction, rotation,...
Definition cl_Tensor.hpp:43
belfem::Tensor::data
T * data()
expose the underlying raw pointer
Definition cl_Tensor.hpp:211
belfem::Tensor::is_3333
bool is_3333() const
returns true if this is a 3x3x3x3 tensor
Definition cl_Tensor.hpp:319
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
fn_gemm.hpp
Matrix-matrix product C := alpha*op(A)*op(B) + beta*C (BLAS ?gemm).
fn_getrf.hpp
LU factorization with partial pivoting, A = P*L*U (LAPACK ?getrf).
fn_getri.hpp
Inverts a square matrix in place from its LU factors (LAPACK ?getri).
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::gemm
void gemm(const char *transa, const char *transb, const int_t *m, const int_t *n, const int_t *k, const T *alpha, const T *a, const int_t *lda, const T *b, const int_t *ldb, const T *beta, T *c, const int_t *ldc)
Definition fn_gemm.hpp:124
belfem::lapack::getrf
void getrf(const int_t *m, const int_t *n, T *a, const int_t *lda, int_t *ipiv, int_t *info)
Definition fn_getrf.hpp:89
belfem::tensor
Definition fn_TR_contract42_arma.hpp:21
belfem::tensor::invert_symmetric
void invert_symmetric(Tensor< T > &aTensor, Vector< T > &aWork, Vector< int_t > &aPivot)
Definition fn_invert_symmetric.hpp:35
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::trans
auto trans(Matrix< T > &aMatrix) -> decltype(trans(aMatrix.matrix_data()))
Definition fn_trans.hpp:74
belfem::MaterialProperty::alpha
@ alpha
Definition cl_Material.hpp:161
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51