My Project
MatrixBlock.hpp
1 /*
2  Copyright 2016 IRIS AS
3  Copyright 2019 NORCE
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_MATRIX_BLOCK_HEADER_INCLUDED
22 #define OPM_MATRIX_BLOCK_HEADER_INCLUDED
23 
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>
30 
31 namespace Dune
32 {
33 namespace FMatrixHelp {
35 template <typename K>
36 static inline K invertMatrix(const FieldMatrix<K,4,4>& matrix, FieldMatrix<K,4,4>& inverse)
37 {
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];
44 
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];
51 
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];
58 
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];
65 
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];
72 
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];
79 
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];
86 
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];
93 
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];
100 
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];
107 
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];
114 
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];
121 
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];
128 
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];
135 
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];
142 
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];
149 
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];
152 
153  // return identity for singular or nearly singular matrices.
154  if (std::abs(det) < 1e-40) {
155  for (int i = 0; i < 4; ++i){
156  inverse[i][i] = 1.0;
157  }
158  return 1.0;
159  }
160  K inv_det = 1.0 / det;
161  inverse *= inv_det;
162 
163  return det;
164 }
165 
166 template <typename K>
167 static inline K invertMatrix(const DynamicMatrix<K>& matrix, DynamicMatrix<K>& inverse)
168 {
169  // this function is only for 4 X 4 matrix
170  assert (matrix.rows() == 4);
171 
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];
178 
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];
185 
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];
192 
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];
199 
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];
206 
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];
213 
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];
220 
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];
227 
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];
234 
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];
241 
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];
248 
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];
255 
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];
262 
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];
269 
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];
276 
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];
283 
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];
286 
287  // return identity for singular or nearly singular matrices.
288  if (std::abs(det) < 1e-40) {
289  for (int i = 0; i < 4; ++i){
290  inverse[i][i] = 1.0;
291  }
292  return 1.0;
293  }
294  K inv_det = 1.0 / det;
295  inverse *= inv_det;
296 
297  return det;
298 }
299 } // end FMatrixHelp
300 
301 namespace ISTLUtility {
302 
304 template <typename K>
305 static inline void invertMatrix(FieldMatrix<K,1,1>& matrix)
306 {
307  FieldMatrix<K,1,1> A ( matrix );
308  FMatrixHelp::invertMatrix(A, matrix );
309 }
310 
312 template <typename K>
313 static inline void invertMatrix(FieldMatrix<K,2,2>& matrix)
314 {
315  FieldMatrix<K,2,2> A ( matrix );
316  FMatrixHelp::invertMatrix(A, matrix );
317 }
318 
320 template <typename K>
321 static inline void invertMatrix(FieldMatrix<K,3,3>& matrix)
322 {
323  FieldMatrix<K,3,3> A ( matrix );
324  FMatrixHelp::invertMatrix(A, matrix );
325 }
326 
328 template <typename K>
329 static inline void invertMatrix(FieldMatrix<K,4,4>& matrix)
330 {
331  FieldMatrix<K,4,4> A ( matrix );
332  FMatrixHelp::invertMatrix(A, matrix );
333 }
334 
336 template <typename K, int n>
337 static inline void invertMatrix(FieldMatrix<K,n,n>& matrix)
338 {
339 #if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
340  Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
341 #endif
342  matrix.invert();
343 }
344 
346 template <typename K>
347 static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
348 {
349  // for 4 X 4 matrix, using the invertMatrix() function above
350  // it is for temporary usage, mainly to reduce the huge burden of testing
351  // what algorithm should be used to invert 4 X 4 matrix will be handled
352  // as a seperate issue
353  if (matrix.rows() == 4) {
354  Dune::DynamicMatrix<K> A = matrix;
355  FMatrixHelp::invertMatrix(A, matrix);
356  return;
357  }
358 
359 #if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
360  Dune::FMatrixPrecision<K>::set_singular_limit(1.e-30);
361 #endif
362  matrix.invert();
363 
364 }
365 
366 } // end ISTLUtility
367 
368 template <class Scalar, int n, int m>
369 class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
370 {
371 public:
372  typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
373 
374  using BaseType :: operator= ;
375  using BaseType :: rows;
376  using BaseType :: cols;
377  explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
378  void invert()
379  {
380  ISTLUtility::invertMatrix( *this );
381  }
382  const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
383  BaseType& asBase() { return static_cast< BaseType& > (*this); }
384 };
385 
386 template<class K, int n, int m>
387 void
388 print_row(std::ostream& s, const MatrixBlock<K,n,m>& A,
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,
392  int precision)
393 {
394  print_row(s, A.asBase(), I, J, therow, width, precision);
395 }
396 
397 template<class K, int n, int m>
398 K& firstmatrixelement(MatrixBlock<K,n,m>& A)
399 {
400  return firstmatrixelement( A.asBase() );
401 }
402 
403 
404 
405 template<typename Scalar, int n, int m>
406 struct MatrixDimension< MatrixBlock< Scalar, n, m > >
407 : public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
408 {
409 };
410 
411 
412 #if HAVE_UMFPACK
413 
417 template<typename T, typename A, int n, int m>
418 class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
419  : public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
420 {
421  typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
422  typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
423 
424 public:
425  typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
426 
427  UMFPack(const RealMatrix& matrix, int verbose, bool)
428  : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
429  {}
430 };
431 #endif
432 
433 #if HAVE_SUPERLU
434 
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> >
441 {
442  typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
443  typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
444 
445 public:
446  typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
447 
448  SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
449  : Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
450  {}
451 };
452 #endif
453 
454 
455 } // end namespace Dune
456 
457 namespace Opm
458 {
459 namespace Detail
460 {
462  template< class TA, class TB, class TC, class PositiveSign >
463  static inline void multMatrixImpl( const TA &A, // n x m
464  const TB &B, // n x p
465  TC &ret, // m x p
466  const PositiveSign )
467  {
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() );
473 
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 )
478  {
479  for( size_type j = 0; j < p; ++j )
480  {
481  K sum = 0;
482  for( size_type k = 0; k < n; ++k )
483  {
484  sum += A[ i ][ k ] * B[ k ][ j ];
485  }
486  // set value depending on given sign
487  ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
488  }
489  }
490  }
491 
495  template< class TA, class TB, class TC, class PositiveSign >
496  static inline void multMatrixTransposedImpl ( const TA &A, // n x m
497  const TB &B, // n x p
498  TC &ret, // m x p
499  const PositiveSign )
500  {
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() );
506 
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 )
511  {
512  for( size_type j = 0; j < p; ++j )
513  {
514  K sum = 0;
515  for( size_type k = 0; k < n; ++k )
516  {
517  sum += A[ k ][ i ] * B[ k ][ j ];
518  }
519  // set value depending on given sign
520  ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
521  }
522  }
523  }
524 
526  template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
527  static inline void multMatrixTransposed(const DenseMatrixA& A,
528  const DenseMatrixB& B,
529  DenseMatrixC& ret)
530  {
531  multMatrixTransposedImpl( A, B, ret, std::true_type() );
532  }
533 
535  template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
536  static inline void negativeMultMatrixTransposed(const DenseMatrixA& A,
537  const DenseMatrixB& B,
538  DenseMatrixC& ret)
539  {
540  multMatrixTransposedImpl( A, B, ret, std::false_type() );
541  }
542 
544  template< class K>
545  static inline void multMatrix(const Dune::DynamicMatrix<K>& A,
546  const Dune::DynamicMatrix<K>& B,
547  Dune::DynamicMatrix<K>& ret )
548  {
549  typedef typename Dune::DynamicMatrix<K> :: size_type size_type;
550 
551  const size_type m = A.rows();
552  const size_type n = A.cols();
553 
554  assert(n == B.rows() );
555 
556  const size_type p = B.cols();
557 
558  ret.resize(m, p);
559 
560  for( size_type i = 0; i < m; ++i )
561  {
562  for( size_type j = 0; j < p; ++j )
563  {
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 ];
567  }
568  }
569  }
570 
573  template <int block_size>
574  struct Inverter
575  {
576  template <typename K>
577  void operator()(const K *matrix [[maybe_unused]], K *inverse [[maybe_unused]])
578  {
579  throw std::logic_error("Not implemented");
580  }
581  };
582 
584  template <>
585  struct Inverter<4>
586  {
587  template <typename K>
588  void operator()(const K *matrix, K *inverse)
589  {
590  // based on Dune::FMatrixHelp::invertMatrix
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];
597 
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];
604 
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];
611 
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];
618 
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];
625 
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];
632 
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];
639 
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];
646 
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];
653 
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];
660 
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];
667 
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];
674 
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];
681 
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];
688 
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];
695 
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];
702 
703  K det = matrix[0] * inverse[0] + matrix[1] * inverse[4] +
704  matrix[2] * inverse[8] + matrix[3] * inverse[12];
705 
706  // return identity for singular or nearly singular matrices.
707  if (std::abs(det) < 1e-40) {
708  for (int i = 0; i < 4; ++i){
709  inverse[4*i + i] = 1.0;
710  }
711  }
712  K inv_det = 1.0 / det;
713 
714  for (unsigned int i = 0; i < 4 * 4; ++i) {
715  inverse[i] *= inv_det;
716  }
717  }
718  };
719 
721  template <>
722  struct Inverter<3>
723  {
724  template <typename K>
725  void operator()(const K *matrix, K *inverse)
726  {
727  // code generated by maple, copied from Dune::DenseMatrix
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];
734 
735  K det = (t4 * matrix[8] - t6 * matrix[7] - t8 * matrix[8] +
736  t10 * matrix[7] + t12 * matrix[5] - t14 * matrix[4]);
737  K t17 = 1.0 / det;
738 
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;
748  }
749  };
750 
752  template <>
753  struct Inverter<2>
754  {
755  template <typename K>
756  void operator()(const K *matrix, K *inverse)
757  {
758  // code based on Dune::DenseMatrix
759  K detinv = matrix[0] * matrix[3] - matrix[1] * matrix[2];
760  detinv = 1 / detinv;
761  inverse[0] = matrix[3] * detinv;
762  inverse[1] = -matrix[1] * detinv;
763  inverse[2] = -matrix[2] * detinv;
764  inverse[3] = matrix[0] * detinv;
765  }
766  };
767 
769  template <>
770  struct Inverter<1>
771  {
772  template <typename K>
773  void operator()(const K *matrix, K *inverse)
774  {
775  inverse[0] = 1.0 / matrix[0];
776  }
777  };
778 
779 } // namespace Detail
780 } // namespace Opm
781 
782 #endif
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