My Project
SmallDenseMatrixUtils.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_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED
22 #define OPM_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED
23 
24 #include <dune/common/dynmatrix.hh>
25 
26 namespace Opm
27 {
28 namespace detail
29 {
31  template< class TA, class TB, class TC, class PositiveSign >
32  static inline void multMatrixImpl( const TA &A, // n x m
33  const TB &B, // n x p
34  TC &ret, // m x p
35  const PositiveSign )
36  {
37  using size_type = typename TA :: size_type;
38  using K = typename TA :: field_type;
39  assert( A.N() == B.N() );
40  assert( A.M() == ret.N() );
41  assert( B.M() == ret.M() );
42 
43  const size_type n = A.N();
44  const size_type m = ret.N();
45  const size_type p = B.M();
46  for( size_type i = 0; i < m; ++i )
47  {
48  for( size_type j = 0; j < p; ++j )
49  {
50  K sum = 0;
51  for( size_type k = 0; k < n; ++k )
52  {
53  sum += A[ i ][ k ] * B[ k ][ j ];
54  }
55  // set value depending on given sign
56  ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
57  }
58  }
59  }
60 
64  template< class TA, class TB, class TC, class PositiveSign >
65  static inline void multMatrixTransposedImpl ( const TA &A, // n x m
66  const TB &B, // n x p
67  TC &ret, // m x p
68  const PositiveSign )
69  {
70  using size_type = typename TA :: size_type;
71  using K = typename TA :: field_type;
72  assert( A.N() == B.N() );
73  assert( A.M() == ret.N() );
74  assert( B.M() == ret.M() );
75 
76  const size_type n = A.N();
77  const size_type m = ret.N();
78  const size_type p = B.M();
79  for( size_type i = 0; i < m; ++i )
80  {
81  for( size_type j = 0; j < p; ++j )
82  {
83  K sum = 0;
84  for( size_type k = 0; k < n; ++k )
85  {
86  sum += A[ k ][ i ] * B[ k ][ j ];
87  }
88  // set value depending on given sign
89  ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
90  }
91  }
92  }
93 
95  template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
96  static inline void multMatrixTransposed(const DenseMatrixA& A,
97  const DenseMatrixB& B,
98  DenseMatrixC& ret)
99  {
100  multMatrixTransposedImpl( A, B, ret, std::true_type() );
101  }
102 
104  template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
105  static inline void negativeMultMatrixTransposed(const DenseMatrixA& A,
106  const DenseMatrixB& B,
107  DenseMatrixC& ret)
108  {
109  multMatrixTransposedImpl( A, B, ret, std::false_type() );
110  }
111 
113  template< class K>
114  static inline void multMatrix(const Dune::DynamicMatrix<K>& A,
115  const Dune::DynamicMatrix<K>& B,
116  Dune::DynamicMatrix<K>& ret )
117  {
118  using size_type = typename Dune::DynamicMatrix<K> :: size_type;
119 
120  const size_type m = A.rows();
121  const size_type n = A.cols();
122 
123  assert(n == B.rows() );
124 
125  const size_type p = B.cols();
126 
127  ret.resize(m, p);
128 
129  for( size_type i = 0; i < m; ++i )
130  {
131  for( size_type j = 0; j < p; ++j )
132  {
133  ret[ i ][ j ] = K( 0 );
134  for( size_type k = 0; k < n; ++k )
135  ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
136  }
137  }
138  }
139 
140 } // namespace detail
141 } // namespace Opm
142 
143 #endif
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27