MatrixInverse.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: MatrixInverse.hpp<2>
//
// Created: Wed Sep  3 14:49:08 2008
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Bj�rn Spjelkavik    <bsp@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.

  This file is part of The Open Reservoir Simulator Project (OpenRS).

  OpenRS 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.

  OpenRS 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 OpenRS.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPENRS_MATRIXINVERSE_HEADER
#define OPENRS_MATRIXINVERSE_HEADER

#include <type_traits>

namespace Opm {


    template <typename M>
    M inverse2x2(const M& m)
    {
        // Because then the divisions below would compile but not be correct, we must guard
        // against integral types.
        typedef typename M::value_type T;
        static_assert(!std::is_integral<T>::value, "");
        assert(m.numRows() == 2 && m.numCols() == 2);

        T det = m(0,0)*m(1,1) - m(0,1)*m(1,0);
        M mi(2, 2, (double*)0);
        mi(0,0) = m(1,1);
        mi(1,0) = -m(1,0);
        mi(0,1) = -m(0,1);
        mi(1,1) = m(0,0);
        mi /= det;
        return mi;
    }

    template <typename M>
    M matprod(const M& m1, const M& m2)
    {
        assert(m1.numCols() == m2.numRows());
        int num_contracting = m1.numCols();
        M m(m1.numRows(), m2.numCols(), (double*)0);
        for (int r = 0; r < m1.numRows(); ++r) {
        for (int c = 0; c < m2.numCols(); ++c) {
            m(r, c) = 0.0;
            for (int kk = 0; kk < num_contracting; ++kk) {
            m(r, c) += m1(r, kk)*m2(kk, c);
            }
        }
        }
        return m;
    }

    template <typename M>
    M inverse3x3(const M& m)
    {
        // Because then the divisions below would compile but not be correct, we must guard
        // against integral types.
        typedef typename M::value_type T;
        static_assert(!std::is_integral<T>::value, "");
        assert(m.numRows() == 3 && m.numCols() == 3);
//      double det = m(0,0)*(m(1,1)*m(2,2)-m(1,2)*m(2,1))
//      - m(0,1)*(m(1,0)*m(2,2)-m(1,2)*m(2,0))
//      + m(0,2)*(m(1,0)*m(2,1)-m(1,1)*m(2,0));

        T a = m(0,0);
        T b = m(0,1);
        T c = m(0,2);
        T d = m(1,0);
        T e = m(1,1);
        T f = m(1,2);
        T g = m(2,0);
        T h = m(2,1);
        T i = m(2,2);
        T t1 = (e-f*h/i);
        T t2 = (c*h/i-b);
        T t3 = (f*g/i-d);
        T t4 = (a-c*g/i);
        T x =  t4*t1-t2*t3;

        M mi(3, 3, (double*)0);
        mi(0,0) = t1/x;
        mi(0,1) = t2/x;
        mi(0,2) = -(c*t1+f*t2)/(i*x);
        mi(1,0) = t3/x;
        mi(1,1) = t4/x;
        mi(1,2) = -(c*t3+f*t4)/(i*x);
        mi(2,0) = -(g*t1+h*t3)/(i*x);
        mi(2,1) = -(g*t2+h*t4)/(i*x);
        mi(2,2) =  1/i+1/(i*i*x)*(c*(g*t1+h*t3)+f*(g*t2+h*t4));
        return mi;
    }


} // namespace Opm

#endif // OPENRS_MATRIXINVERSE_HEADER