SmallDenseMatrixUtils.hpp

Go to the documentation of this file.

/*
  Copyright 2016 IRIS AS
  Copyright 2019 NORCE
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED
#define OPM_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED
 
#include <dune/common/dynmatrix.hh>
 
namespace Opm
{
namespace detail
{
    template< class TA, class TB, class TC, class PositiveSign >
    static inline void multMatrixImpl( const TA &A, // n x m
                                       const TB &B, // n x p
                                       TC &ret,     // m x p
                                       const PositiveSign )
    {
        using size_type = typename TA :: size_type;
        using K = typename TA :: field_type;
        assert( A.N() == B.N() );
        assert( A.M() == ret.N() );
        assert( B.M() == ret.M() );
 
        const size_type n = A.N();
        const size_type m = ret.N();
        const size_type p = B.M();
        for( size_type i = 0; i < m; ++i )
        {
            for( size_type j = 0; j < p; ++j )
            {
                K sum = 0;
                for( size_type k = 0; k < n; ++k )
                {
                    sum += A[ i ][ k ] * B[ k ][ j ];
                }
                // set value depending on given sign
                ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
            }
        }
    }
 
    template< class TA, class TB, class TC, class PositiveSign >
    static inline void multMatrixTransposedImpl ( const TA &A, // n x m
                                                  const TB &B, // n x p
                                                  TC &ret,     // m x p
                                                  const PositiveSign )
    {
        using size_type = typename TA :: size_type;
        using K = typename TA :: field_type;
        assert( A.N() == B.N() );
        assert( A.M() == ret.N() );
        assert( B.M() == ret.M() );
 
        const size_type n = A.N();
        const size_type m = ret.N();
        const size_type p = B.M();
        for( size_type i = 0; i < m; ++i )
        {
            for( size_type j = 0; j < p; ++j )
            {
                K sum = 0;
                for( size_type k = 0; k < n; ++k )
                {
                    sum += A[ k ][ i ] * B[ k ][ j ];
                }
                // set value depending on given sign
                ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
            }
        }
    }
 
    template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
    static inline void multMatrixTransposed(const DenseMatrixA& A,
                                            const DenseMatrixB& B,
                                            DenseMatrixC& ret)
    {
        multMatrixTransposedImpl( A, B, ret, std::true_type() );
    }
 
    template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
    static inline void negativeMultMatrixTransposed(const DenseMatrixA& A,
                                                    const DenseMatrixB& B,
                                                    DenseMatrixC& ret)
    {
        multMatrixTransposedImpl( A, B, ret, std::false_type() );
    }
 
    template< class K>
    static inline
    Dune::DynamicMatrix<K>
    multMatrix(const Dune::DynamicMatrix<K>& A,
               const Dune::DynamicMatrix<K>& B)
    {
        using size_type = typename Dune::DynamicMatrix<K> :: size_type;
 
        const size_type m = A.rows();
        const size_type n = A.cols();
 
        assert(n == B.rows() );
 
        const size_type p = B.cols();
 
        Dune::DynamicMatrix<K> ret(m, p);
 
        for( size_type i = 0; i < m; ++i )
        {
            for( size_type j = 0; j < p; ++j )
            {
                ret[ i ][ j ] = K( 0 );
                for( size_type k = 0; k < n; ++k )
                    ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
            }
        }
 
        return ret;
    }
 
} // namespace detail
} // namespace Opm
 
#endif