21 #ifndef OPM_MATRIX_BLOCK_HEADER_INCLUDED
22 #define OPM_MATRIX_BLOCK_HEADER_INCLUDED
24 #include <dune/common/fmatrix.hh>
25 #include <dune/common/fvector.hh>
26 #include <dune/common/version.hh>
27 #include <dune/istl/matrixutils.hh>
28 #include <dune/istl/umfpack.hh>
29 #include <dune/istl/superlu.hh>
33 namespace FMatrixHelp {
36 static inline K invertMatrix(
const FieldMatrix<K,4,4>& matrix, FieldMatrix<K,4,4>& inverse)
38 inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
39 matrix[1][1] * matrix[2][3] * matrix[3][2] -
40 matrix[2][1] * matrix[1][2] * matrix[3][3] +
41 matrix[2][1] * matrix[1][3] * matrix[3][2] +
42 matrix[3][1] * matrix[1][2] * matrix[2][3] -
43 matrix[3][1] * matrix[1][3] * matrix[2][2];
45 inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
46 matrix[1][0] * matrix[2][3] * matrix[3][2] +
47 matrix[2][0] * matrix[1][2] * matrix[3][3] -
48 matrix[2][0] * matrix[1][3] * matrix[3][2] -
49 matrix[3][0] * matrix[1][2] * matrix[2][3] +
50 matrix[3][0] * matrix[1][3] * matrix[2][2];
52 inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
53 matrix[1][0] * matrix[2][3] * matrix[3][1] -
54 matrix[2][0] * matrix[1][1] * matrix[3][3] +
55 matrix[2][0] * matrix[1][3] * matrix[3][1] +
56 matrix[3][0] * matrix[1][1] * matrix[2][3] -
57 matrix[3][0] * matrix[1][3] * matrix[2][1];
59 inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
60 matrix[1][0] * matrix[2][2] * matrix[3][1] +
61 matrix[2][0] * matrix[1][1] * matrix[3][2] -
62 matrix[2][0] * matrix[1][2] * matrix[3][1] -
63 matrix[3][0] * matrix[1][1] * matrix[2][2] +
64 matrix[3][0] * matrix[1][2] * matrix[2][1];
66 inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
67 matrix[0][1] * matrix[2][3] * matrix[3][2] +
68 matrix[2][1] * matrix[0][2] * matrix[3][3] -
69 matrix[2][1] * matrix[0][3] * matrix[3][2] -
70 matrix[3][1] * matrix[0][2] * matrix[2][3] +
71 matrix[3][1] * matrix[0][3] * matrix[2][2];
73 inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
74 matrix[0][0] * matrix[2][3] * matrix[3][2] -
75 matrix[2][0] * matrix[0][2] * matrix[3][3] +
76 matrix[2][0] * matrix[0][3] * matrix[3][2] +
77 matrix[3][0] * matrix[0][2] * matrix[2][3] -
78 matrix[3][0] * matrix[0][3] * matrix[2][2];
80 inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
81 matrix[0][0] * matrix[2][3] * matrix[3][1] +
82 matrix[2][0] * matrix[0][1] * matrix[3][3] -
83 matrix[2][0] * matrix[0][3] * matrix[3][1] -
84 matrix[3][0] * matrix[0][1] * matrix[2][3] +
85 matrix[3][0] * matrix[0][3] * matrix[2][1];
87 inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
88 matrix[0][0] * matrix[2][2] * matrix[3][1] -
89 matrix[2][0] * matrix[0][1] * matrix[3][2] +
90 matrix[2][0] * matrix[0][2] * matrix[3][1] +
91 matrix[3][0] * matrix[0][1] * matrix[2][2] -
92 matrix[3][0] * matrix[0][2] * matrix[2][1];
94 inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
95 matrix[0][1] * matrix[1][3] * matrix[3][2] -
96 matrix[1][1] * matrix[0][2] * matrix[3][3] +
97 matrix[1][1] * matrix[0][3] * matrix[3][2] +
98 matrix[3][1] * matrix[0][2] * matrix[1][3] -
99 matrix[3][1] * matrix[0][3] * matrix[1][2];
101 inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
102 matrix[0][0] * matrix[1][3] * matrix[3][2] +
103 matrix[1][0] * matrix[0][2] * matrix[3][3] -
104 matrix[1][0] * matrix[0][3] * matrix[3][2] -
105 matrix[3][0] * matrix[0][2] * matrix[1][3] +
106 matrix[3][0] * matrix[0][3] * matrix[1][2];
108 inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
109 matrix[0][0] * matrix[1][3] * matrix[3][1] -
110 matrix[1][0] * matrix[0][1] * matrix[3][3] +
111 matrix[1][0] * matrix[0][3] * matrix[3][1] +
112 matrix[3][0] * matrix[0][1] * matrix[1][3] -
113 matrix[3][0] * matrix[0][3] * matrix[1][1];
115 inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
116 matrix[0][0] * matrix[1][2] * matrix[3][1] +
117 matrix[1][0] * matrix[0][1] * matrix[3][2] -
118 matrix[1][0] * matrix[0][2] * matrix[3][1] -
119 matrix[3][0] * matrix[0][1] * matrix[1][2] +
120 matrix[3][0] * matrix[0][2] * matrix[1][1];
122 inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
123 matrix[0][1] * matrix[1][3] * matrix[2][2] +
124 matrix[1][1] * matrix[0][2] * matrix[2][3] -
125 matrix[1][1] * matrix[0][3] * matrix[2][2] -
126 matrix[2][1] * matrix[0][2] * matrix[1][3] +
127 matrix[2][1] * matrix[0][3] * matrix[1][2];
129 inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
130 matrix[0][0] * matrix[1][3] * matrix[2][2] -
131 matrix[1][0] * matrix[0][2] * matrix[2][3] +
132 matrix[1][0] * matrix[0][3] * matrix[2][2] +
133 matrix[2][0] * matrix[0][2] * matrix[1][3] -
134 matrix[2][0] * matrix[0][3] * matrix[1][2];
136 inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
137 matrix[0][0] * matrix[1][3] * matrix[2][1] +
138 matrix[1][0] * matrix[0][1] * matrix[2][3] -
139 matrix[1][0] * matrix[0][3] * matrix[2][1] -
140 matrix[2][0] * matrix[0][1] * matrix[1][3] +
141 matrix[2][0] * matrix[0][3] * matrix[1][1];
143 inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
144 matrix[0][0] * matrix[1][2] * matrix[2][1] -
145 matrix[1][0] * matrix[0][1] * matrix[2][2] +
146 matrix[1][0] * matrix[0][2] * matrix[2][1] +
147 matrix[2][0] * matrix[0][1] * matrix[1][2] -
148 matrix[2][0] * matrix[0][2] * matrix[1][1];
150 K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] +
151 matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
154 if (std::abs(det) < 1e-40) {
155 for (
int i = 0; i < 4; ++i){
160 K inv_det = 1.0 / det;
166 template <
typename K>
167 static inline K invertMatrix(
const DynamicMatrix<K>& matrix, DynamicMatrix<K>& inverse)
170 assert (matrix.rows() == 4);
172 inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
173 matrix[1][1] * matrix[2][3] * matrix[3][2] -
174 matrix[2][1] * matrix[1][2] * matrix[3][3] +
175 matrix[2][1] * matrix[1][3] * matrix[3][2] +
176 matrix[3][1] * matrix[1][2] * matrix[2][3] -
177 matrix[3][1] * matrix[1][3] * matrix[2][2];
179 inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
180 matrix[1][0] * matrix[2][3] * matrix[3][2] +
181 matrix[2][0] * matrix[1][2] * matrix[3][3] -
182 matrix[2][0] * matrix[1][3] * matrix[3][2] -
183 matrix[3][0] * matrix[1][2] * matrix[2][3] +
184 matrix[3][0] * matrix[1][3] * matrix[2][2];
186 inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
187 matrix[1][0] * matrix[2][3] * matrix[3][1] -
188 matrix[2][0] * matrix[1][1] * matrix[3][3] +
189 matrix[2][0] * matrix[1][3] * matrix[3][1] +
190 matrix[3][0] * matrix[1][1] * matrix[2][3] -
191 matrix[3][0] * matrix[1][3] * matrix[2][1];
193 inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
194 matrix[1][0] * matrix[2][2] * matrix[3][1] +
195 matrix[2][0] * matrix[1][1] * matrix[3][2] -
196 matrix[2][0] * matrix[1][2] * matrix[3][1] -
197 matrix[3][0] * matrix[1][1] * matrix[2][2] +
198 matrix[3][0] * matrix[1][2] * matrix[2][1];
200 inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
201 matrix[0][1] * matrix[2][3] * matrix[3][2] +
202 matrix[2][1] * matrix[0][2] * matrix[3][3] -
203 matrix[2][1] * matrix[0][3] * matrix[3][2] -
204 matrix[3][1] * matrix[0][2] * matrix[2][3] +
205 matrix[3][1] * matrix[0][3] * matrix[2][2];
207 inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
208 matrix[0][0] * matrix[2][3] * matrix[3][2] -
209 matrix[2][0] * matrix[0][2] * matrix[3][3] +
210 matrix[2][0] * matrix[0][3] * matrix[3][2] +
211 matrix[3][0] * matrix[0][2] * matrix[2][3] -
212 matrix[3][0] * matrix[0][3] * matrix[2][2];
214 inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
215 matrix[0][0] * matrix[2][3] * matrix[3][1] +
216 matrix[2][0] * matrix[0][1] * matrix[3][3] -
217 matrix[2][0] * matrix[0][3] * matrix[3][1] -
218 matrix[3][0] * matrix[0][1] * matrix[2][3] +
219 matrix[3][0] * matrix[0][3] * matrix[2][1];
221 inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
222 matrix[0][0] * matrix[2][2] * matrix[3][1] -
223 matrix[2][0] * matrix[0][1] * matrix[3][2] +
224 matrix[2][0] * matrix[0][2] * matrix[3][1] +
225 matrix[3][0] * matrix[0][1] * matrix[2][2] -
226 matrix[3][0] * matrix[0][2] * matrix[2][1];
228 inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
229 matrix[0][1] * matrix[1][3] * matrix[3][2] -
230 matrix[1][1] * matrix[0][2] * matrix[3][3] +
231 matrix[1][1] * matrix[0][3] * matrix[3][2] +
232 matrix[3][1] * matrix[0][2] * matrix[1][3] -
233 matrix[3][1] * matrix[0][3] * matrix[1][2];
235 inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
236 matrix[0][0] * matrix[1][3] * matrix[3][2] +
237 matrix[1][0] * matrix[0][2] * matrix[3][3] -
238 matrix[1][0] * matrix[0][3] * matrix[3][2] -
239 matrix[3][0] * matrix[0][2] * matrix[1][3] +
240 matrix[3][0] * matrix[0][3] * matrix[1][2];
242 inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
243 matrix[0][0] * matrix[1][3] * matrix[3][1] -
244 matrix[1][0] * matrix[0][1] * matrix[3][3] +
245 matrix[1][0] * matrix[0][3] * matrix[3][1] +
246 matrix[3][0] * matrix[0][1] * matrix[1][3] -
247 matrix[3][0] * matrix[0][3] * matrix[1][1];
249 inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
250 matrix[0][0] * matrix[1][2] * matrix[3][1] +
251 matrix[1][0] * matrix[0][1] * matrix[3][2] -
252 matrix[1][0] * matrix[0][2] * matrix[3][1] -
253 matrix[3][0] * matrix[0][1] * matrix[1][2] +
254 matrix[3][0] * matrix[0][2] * matrix[1][1];
256 inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
257 matrix[0][1] * matrix[1][3] * matrix[2][2] +
258 matrix[1][1] * matrix[0][2] * matrix[2][3] -
259 matrix[1][1] * matrix[0][3] * matrix[2][2] -
260 matrix[2][1] * matrix[0][2] * matrix[1][3] +
261 matrix[2][1] * matrix[0][3] * matrix[1][2];
263 inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
264 matrix[0][0] * matrix[1][3] * matrix[2][2] -
265 matrix[1][0] * matrix[0][2] * matrix[2][3] +
266 matrix[1][0] * matrix[0][3] * matrix[2][2] +
267 matrix[2][0] * matrix[0][2] * matrix[1][3] -
268 matrix[2][0] * matrix[0][3] * matrix[1][2];
270 inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
271 matrix[0][0] * matrix[1][3] * matrix[2][1] +
272 matrix[1][0] * matrix[0][1] * matrix[2][3] -
273 matrix[1][0] * matrix[0][3] * matrix[2][1] -
274 matrix[2][0] * matrix[0][1] * matrix[1][3] +
275 matrix[2][0] * matrix[0][3] * matrix[1][1];
277 inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
278 matrix[0][0] * matrix[1][2] * matrix[2][1] -
279 matrix[1][0] * matrix[0][1] * matrix[2][2] +
280 matrix[1][0] * matrix[0][2] * matrix[2][1] +
281 matrix[2][0] * matrix[0][1] * matrix[1][2] -
282 matrix[2][0] * matrix[0][2] * matrix[1][1];
284 K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] +
285 matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
288 if (std::abs(det) < 1e-40) {
289 for (
int i = 0; i < 4; ++i){
294 K inv_det = 1.0 / det;
301 namespace ISTLUtility {
304 template <
typename K>
305 static inline void invertMatrix(FieldMatrix<K,1,1>& matrix)
307 FieldMatrix<K,1,1> A ( matrix );
308 FMatrixHelp::invertMatrix(A, matrix );
312 template <
typename K>
313 static inline void invertMatrix(FieldMatrix<K,2,2>& matrix)
315 FieldMatrix<K,2,2> A ( matrix );
316 FMatrixHelp::invertMatrix(A, matrix );
320 template <
typename K>
321 static inline void invertMatrix(FieldMatrix<K,3,3>& matrix)
323 FieldMatrix<K,3,3> A ( matrix );
324 FMatrixHelp::invertMatrix(A, matrix );
328 template <
typename K>
329 static inline void invertMatrix(FieldMatrix<K,4,4>& matrix)
331 FieldMatrix<K,4,4> A ( matrix );
332 FMatrixHelp::invertMatrix(A, matrix );
336 template <
typename K,
int n>
337 static inline void invertMatrix(FieldMatrix<K,n,n>& matrix)
339 #if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
340 Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
346 template <
typename K>
347 static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
353 if (matrix.rows() == 4) {
354 Dune::DynamicMatrix<K> A = matrix;
355 FMatrixHelp::invertMatrix(A, matrix);
359 #if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
360 Dune::FMatrixPrecision<K>::set_singular_limit(1.e-30);
368 template <
class Scalar,
int n,
int m>
372 typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
374 using BaseType :: operator= ;
375 using BaseType :: rows;
376 using BaseType :: cols;
377 explicit MatrixBlock(
const Scalar scalar = 0 ) : BaseType( scalar ) {}
380 ISTLUtility::invertMatrix( *
this );
382 const BaseType& asBase()
const {
return static_cast< const BaseType&
> (*this); }
383 BaseType& asBase() {
return static_cast< BaseType&
> (*this); }
386 template<
class K,
int n,
int m>
389 typename FieldMatrix<K,n,m>::size_type I,
390 typename FieldMatrix<K,n,m>::size_type J,
391 typename FieldMatrix<K,n,m>::size_type therow,
int width,
394 print_row(s, A.asBase(), I, J, therow, width, precision);
397 template<
class K,
int n,
int m>
398 K& firstmatrixelement(MatrixBlock<K,n,m>& A)
400 return firstmatrixelement( A.asBase() );
405 template<
typename Scalar,
int n,
int m>
407 :
public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
417 template<
typename T,
typename A,
int n,
int m>
419 :
public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
421 typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
422 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
425 typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
427 UMFPack(
const RealMatrix& matrix,
int verbose,
bool)
428 : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
438 template<
typename T,
typename A,
int n,
int m>
439 class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
440 :
public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
442 typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
443 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
446 typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
448 SuperLU(
const RealMatrix& matrix,
int verbose,
bool reuse=
true)
449 : Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
462 template<
class TA,
class TB,
class TC,
class PositiveSign >
463 static inline void multMatrixImpl(
const TA &A,
468 typedef typename TA :: size_type size_type;
469 typedef typename TA :: field_type K;
470 assert( A.N() == B.N() );
471 assert( A.M() == ret.N() );
472 assert( B.M() == ret.M() );
474 const size_type n = A.N();
475 const size_type m = ret.N();
476 const size_type p = B.M();
477 for( size_type i = 0; i < m; ++i )
479 for( size_type j = 0; j < p; ++j )
482 for( size_type k = 0; k < n; ++k )
484 sum += A[ i ][ k ] * B[ k ][ j ];
487 ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
495 template<
class TA,
class TB,
class TC,
class PositiveSign >
496 static inline void multMatrixTransposedImpl (
const TA &A,
501 typedef typename TA :: size_type size_type;
502 typedef typename TA :: field_type K;
503 assert( A.N() == B.N() );
504 assert( A.M() == ret.N() );
505 assert( B.M() == ret.M() );
507 const size_type n = A.N();
508 const size_type m = ret.N();
509 const size_type p = B.M();
510 for( size_type i = 0; i < m; ++i )
512 for( size_type j = 0; j < p; ++j )
515 for( size_type k = 0; k < n; ++k )
517 sum += A[ k ][ i ] * B[ k ][ j ];
520 ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
526 template <
class DenseMatrixA,
class DenseMatrixB,
class DenseMatrixC>
527 static inline void multMatrixTransposed(
const DenseMatrixA& A,
528 const DenseMatrixB& B,
531 multMatrixTransposedImpl( A, B, ret, std::true_type() );
535 template <
class DenseMatrixA,
class DenseMatrixB,
class DenseMatrixC>
536 static inline void negativeMultMatrixTransposed(
const DenseMatrixA& A,
537 const DenseMatrixB& B,
540 multMatrixTransposedImpl( A, B, ret, std::false_type() );
545 static inline void multMatrix(
const Dune::DynamicMatrix<K>& A,
546 const Dune::DynamicMatrix<K>& B,
547 Dune::DynamicMatrix<K>& ret )
549 typedef typename Dune::DynamicMatrix<K> :: size_type size_type;
551 const size_type m = A.rows();
552 const size_type n = A.cols();
554 assert(n == B.rows() );
556 const size_type p = B.cols();
560 for( size_type i = 0; i < m; ++i )
562 for( size_type j = 0; j < p; ++j )
564 ret[ i ][ j ] = K( 0 );
565 for( size_type k = 0; k < n; ++k )
566 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
573 template <
int block_size>
576 template <
typename K>
577 void operator()(
const K *matrix [[maybe_unused]], K *inverse [[maybe_unused]])
579 throw std::logic_error(
"Not implemented");
587 template <
typename K>
588 void operator()(
const K *matrix, K *inverse)
591 inverse[0] = matrix[5] * matrix[10] * matrix[15] -
592 matrix[5] * matrix[11] * matrix[14] -
593 matrix[9] * matrix[6] * matrix[15] +
594 matrix[9] * matrix[7] * matrix[14] +
595 matrix[13] * matrix[6] * matrix[11] -
596 matrix[13] * matrix[7] * matrix[10];
598 inverse[4] = -matrix[4] * matrix[10] * matrix[15] +
599 matrix[4] * matrix[11] * matrix[14] +
600 matrix[8] * matrix[6] * matrix[15] -
601 matrix[8] * matrix[7] * matrix[14] -
602 matrix[12] * matrix[6] * matrix[11] +
603 matrix[12] * matrix[7] * matrix[10];
605 inverse[8] = matrix[4] * matrix[9] * matrix[15] -
606 matrix[4] * matrix[11] * matrix[13] -
607 matrix[8] * matrix[5] * matrix[15] +
608 matrix[8] * matrix[7] * matrix[13] +
609 matrix[12] * matrix[5] * matrix[11] -
610 matrix[12] * matrix[7] * matrix[9];
612 inverse[12] = -matrix[4] * matrix[9] * matrix[14] +
613 matrix[4] * matrix[10] * matrix[13] +
614 matrix[8] * matrix[5] * matrix[14] -
615 matrix[8] * matrix[6] * matrix[13] -
616 matrix[12] * matrix[5] * matrix[10] +
617 matrix[12] * matrix[6] * matrix[9];
619 inverse[1]= -matrix[1] * matrix[10] * matrix[15] +
620 matrix[1] * matrix[11] * matrix[14] +
621 matrix[9] * matrix[2] * matrix[15] -
622 matrix[9] * matrix[3] * matrix[14] -
623 matrix[13] * matrix[2] * matrix[11] +
624 matrix[13] * matrix[3] * matrix[10];
626 inverse[5] = matrix[0] * matrix[10] * matrix[15] -
627 matrix[0] * matrix[11] * matrix[14] -
628 matrix[8] * matrix[2] * matrix[15] +
629 matrix[8] * matrix[3] * matrix[14] +
630 matrix[12] * matrix[2] * matrix[11] -
631 matrix[12] * matrix[3] * matrix[10];
633 inverse[9] = -matrix[0] * matrix[9] * matrix[15] +
634 matrix[0] * matrix[11] * matrix[13] +
635 matrix[8] * matrix[1] * matrix[15] -
636 matrix[8] * matrix[3] * matrix[13] -
637 matrix[12] * matrix[1] * matrix[11] +
638 matrix[12] * matrix[3] * matrix[9];
640 inverse[13] = matrix[0] * matrix[9] * matrix[14] -
641 matrix[0] * matrix[10] * matrix[13] -
642 matrix[8] * matrix[1] * matrix[14] +
643 matrix[8] * matrix[2] * matrix[13] +
644 matrix[12] * matrix[1] * matrix[10] -
645 matrix[12] * matrix[2] * matrix[9];
647 inverse[2] = matrix[1] * matrix[6] * matrix[15] -
648 matrix[1] * matrix[7] * matrix[14] -
649 matrix[5] * matrix[2] * matrix[15] +
650 matrix[5] * matrix[3] * matrix[14] +
651 matrix[13] * matrix[2] * matrix[7] -
652 matrix[13] * matrix[3] * matrix[6];
654 inverse[6] = -matrix[0] * matrix[6] * matrix[15] +
655 matrix[0] * matrix[7] * matrix[14] +
656 matrix[4] * matrix[2] * matrix[15] -
657 matrix[4] * matrix[3] * matrix[14] -
658 matrix[12] * matrix[2] * matrix[7] +
659 matrix[12] * matrix[3] * matrix[6];
661 inverse[10] = matrix[0] * matrix[5] * matrix[15] -
662 matrix[0] * matrix[7] * matrix[13] -
663 matrix[4] * matrix[1] * matrix[15] +
664 matrix[4] * matrix[3] * matrix[13] +
665 matrix[12] * matrix[1] * matrix[7] -
666 matrix[12] * matrix[3] * matrix[5];
668 inverse[14] = -matrix[0] * matrix[5] * matrix[14] +
669 matrix[0] * matrix[6] * matrix[13] +
670 matrix[4] * matrix[1] * matrix[14] -
671 matrix[4] * matrix[2] * matrix[13] -
672 matrix[12] * matrix[1] * matrix[6] +
673 matrix[12] * matrix[2] * matrix[5];
675 inverse[3] = -matrix[1] * matrix[6] * matrix[11] +
676 matrix[1] * matrix[7] * matrix[10] +
677 matrix[5] * matrix[2] * matrix[11] -
678 matrix[5] * matrix[3] * matrix[10] -
679 matrix[9] * matrix[2] * matrix[7] +
680 matrix[9] * matrix[3] * matrix[6];
682 inverse[7] = matrix[0] * matrix[6] * matrix[11] -
683 matrix[0] * matrix[7] * matrix[10] -
684 matrix[4] * matrix[2] * matrix[11] +
685 matrix[4] * matrix[3] * matrix[10] +
686 matrix[8] * matrix[2] * matrix[7] -
687 matrix[8] * matrix[3] * matrix[6];
689 inverse[11] = -matrix[0] * matrix[5] * matrix[11] +
690 matrix[0] * matrix[7] * matrix[9] +
691 matrix[4] * matrix[1] * matrix[11] -
692 matrix[4] * matrix[3] * matrix[9] -
693 matrix[8] * matrix[1] * matrix[7] +
694 matrix[8] * matrix[3] * matrix[5];
696 inverse[15] = matrix[0] * matrix[5] * matrix[10] -
697 matrix[0] * matrix[6] * matrix[9] -
698 matrix[4] * matrix[1] * matrix[10] +
699 matrix[4] * matrix[2] * matrix[9] +
700 matrix[8] * matrix[1] * matrix[6] -
701 matrix[8] * matrix[2] * matrix[5];
703 K det = matrix[0] * inverse[0] + matrix[1] * inverse[4] +
704 matrix[2] * inverse[8] + matrix[3] * inverse[12];
707 if (std::abs(det) < 1e-40) {
708 for (
int i = 0; i < 4; ++i){
709 inverse[4*i + i] = 1.0;
712 K inv_det = 1.0 / det;
714 for (
unsigned int i = 0; i < 4 * 4; ++i) {
715 inverse[i] *= inv_det;
724 template <
typename K>
725 void operator()(
const K *matrix, K *inverse)
728 K t4 = matrix[0] * matrix[4];
729 K t6 = matrix[0] * matrix[5];
730 K t8 = matrix[1] * matrix[3];
731 K t10 = matrix[2] * matrix[3];
732 K t12 = matrix[1] * matrix[6];
733 K t14 = matrix[2] * matrix[6];
735 K det = (t4 * matrix[8] - t6 * matrix[7] - t8 * matrix[8] +
736 t10 * matrix[7] + t12 * matrix[5] - t14 * matrix[4]);
739 inverse[0] = (matrix[4] * matrix[8] - matrix[5] * matrix[7]) * t17;
740 inverse[1] = -(matrix[1] * matrix[8] - matrix[2] * matrix[7]) * t17;
741 inverse[2] = (matrix[1] * matrix[5] - matrix[2] * matrix[4]) * t17;
742 inverse[3] = -(matrix[3] * matrix[8] - matrix[5] * matrix[6]) * t17;
743 inverse[4] = (matrix[0] * matrix[8] - t14) * t17;
744 inverse[5] = -(t6 - t10) * t17;
745 inverse[6] = (matrix[3] * matrix[7] - matrix[4] * matrix[6]) * t17;
746 inverse[7] = -(matrix[0] * matrix[7] - t12) * t17;
747 inverse[8] = (t4 - t8) * t17;
755 template <
typename K>
756 void operator()(
const K *matrix, K *inverse)
759 K detinv = matrix[0] * matrix[3] - matrix[1] * matrix[2];
761 inverse[0] = matrix[3] * detinv;
762 inverse[1] = -matrix[1] * detinv;
763 inverse[2] = -matrix[2] * detinv;
764 inverse[3] = matrix[0] * detinv;
772 template <
typename K>
773 void operator()(
const K *matrix, K *inverse)
775 inverse[0] = 1.0 / matrix[0];
Definition: MatrixBlock.hpp:370
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
perform out of place matrix inversion on C-style arrays must have a specified block_size
Definition: MatrixBlock.hpp:575