Matrix.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: Matrix.hpp
//
// Created: Tue Jun 30 11:25:46 2009
//
// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
//            Atgeirr F Rasmussen <atgeirr@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_MATRIX_HEADER
#define OPENRS_MATRIX_HEADER

#include <algorithm>
#include <ostream>
#include <vector>
#include <boost/bind.hpp>

#include <dune/common/fvector.hh>
#include <opm/common/ErrorMacros.hpp>

#include <opm/porsol/common/fortran.hpp>
#include <opm/porsol/common/blas_lapack.hpp>

namespace Opm {

    // ----------------------------------------------------------------------
    // FullMatrix storage policies.
    //

    template<typename T>
    class OwnData {
    public:
        T&       operator[](int i)       { return data_[i]; }
        const T& operator[](int i) const { return data_[i]; }


        int size() const { return data_.size(); }

        T*       data()       { return &data_[0]; }
        const T* data() const { return &data_[0]; }

    protected:
        OwnData(int sz, const T* data)
        {
            if (data) {
                data_.assign(data, data + sz);
            } else {
                data_.resize(sz);
            }
        }

    private:
        std::vector<T> data_;
    };

    template<typename T>
    class SharedData {
    public:
        T&       operator[](int i)       { return data_[i]; }
        const T& operator[](int i) const { return data_[i]; }

        int size() const { return sz_; }

        T*       data()       { return data_; }
        const T* data() const { return data_; }

    protected:
        SharedData(int sz, T* data)
            : sz_(sz), data_(data)
        {
            assert ((sz == 0) == (data == 0));
        }

    private:
        int sz_;
        T*  data_;
    };

    template<typename T>
    class ImmutableSharedData {
    public:
        const T& operator[](int i) const { return data_[i]; }

        int size() const { return sz_; }

        const T* data() const { return data_; }

    protected:
        ImmutableSharedData(int sz, const T* data)
            : sz_(sz), data_(data)
        {
            assert (data_ != 0);
        }

    private:
        int sz_;
        const T*  data_;
    };





    // ----------------------------------------------------------------------
    // FullMatrix ordering policies.
    //

    class OrderingBase {
    public:
        int numRows() const { return rows_; }

        int numCols() const { return cols_; }

    protected:
        OrderingBase()
            : rows_(0), cols_(0)
        {}

        OrderingBase(int rows, int cols)
            : rows_(rows), cols_(cols)
        {}

    private:
        int rows_, cols_;
    };


    class COrdering : public OrderingBase {
    public:
        int leadingDimension() const { return numCols(); }

        int idx(int row, int col) const
        {
            assert ((0 <= row) && (row < numRows()));
            assert ((0 <= col) && (col < numCols()));

            return row*numCols() + col;
        }

    protected:
        COrdering()
            : OrderingBase()
        {}

        COrdering(int rows, int cols)
            : OrderingBase(rows, cols)
        {}
    };


    class FortranOrdering : public OrderingBase {
    public:
        int leadingDimension() const { return numRows(); }

        int idx(int row, int col) const
        {
            assert ((0 <= row) && (row < numRows()));
            assert ((0 <= col) && (col < numCols()));

            return row + col*numRows();
        }

    protected:
        FortranOrdering()
            : OrderingBase()
        {}

        FortranOrdering(int rows, int cols)
            : OrderingBase(rows, cols)
        {}
    };




    // ----------------------------------------------------------------------
    // Class FullMatrix.
    //


    template<typename                 T,
             template<typename> class StoragePolicy,
             class                    OrderingPolicy>
    class FullMatrix : private StoragePolicy<T>,
                       private OrderingPolicy
    {
    public:
        FullMatrix()
            : StoragePolicy<T>(0, 0),
              OrderingPolicy()
        {}

        template <typename DataPointer>
        FullMatrix(int rows, int cols, DataPointer data)
            : StoragePolicy<T>(rows * cols, data),
              OrderingPolicy(rows, cols)
        {}

        template <template<typename> class OtherSP>
        explicit FullMatrix(const FullMatrix<T, OtherSP, OrderingPolicy>& m)
            : StoragePolicy<T>(m.numRows()*m.numCols(), m.data()),
              OrderingPolicy(m.numRows(), m.numCols())
        {
        }

        template <template<typename> class OtherSP, class OtherOP>
        FullMatrix& operator=(const FullMatrix<T, OtherSP, OtherOP>& m)
        {
            assert(numRows() == m.numRows());
            assert(numCols() == m.numCols());
            for (int r = 0; r < numRows(); ++r) {
                for (int c = 0; c < numCols(); ++c) {
                    this->operator()(r, c) = m(r,c);
                }
            }
            return *this;
        }

        template <template<typename> class OtherSP, class OtherOP>
        bool operator==(const FullMatrix<T, OtherSP, OtherOP>& m)
        {
            if (numRows() != m.numRows() || numCols() != m.numCols()) {
                return false;
            }
            for (int r = 0; r < numRows(); ++r) {
                for (int c = 0; c < numCols(); ++c) {
                    if (this->operator()(r,c) != m(r,c)) {
                        return false;
                    }
                }
            }
            return true;
        }

        template <template<typename> class OtherSP>
        void operator+= (const FullMatrix<T, OtherSP, OrderingPolicy>& m)
        {
            assert(numRows() == m.numRows() && numCols() == m.numCols());
            std::transform(data(), data() + this->size(),
                           m.data(), data(), std::plus<T>());
        }

        void operator*= (const T& scalar)
        {
            std::transform(data(), data() + this->size(),
                           data(), boost::bind(std::multiplies<T>(), _1, scalar));
        }

        typedef T value_type;

        using StoragePolicy<T>::data;
        using OrderingPolicy::numRows;
        using OrderingPolicy::numCols;
        using OrderingPolicy::leadingDimension;

        value_type&       operator()(int row, int col)
        {
            return this->operator[](this->idx(row, col));
        }

        const value_type& operator()(int row, int col) const
        {
            return this->operator[](this->idx(row, col));
        }
    };



    // ----------------------------------------------------------------------
    // FullMatrix operations.
    //


    // Convenience typedefs.

    typedef FullMatrix<double, OwnData,             COrdering>        OwnCMatrix;
    typedef FullMatrix<double, SharedData,          COrdering>        SharedCMatrix;
    typedef const FullMatrix<double, ImmutableSharedData, COrdering>  ImmutableCMatrix;


    typedef FullMatrix<double, OwnData,             FortranOrdering>       OwnFortranMatrix;
    typedef FullMatrix<double, SharedData,          FortranOrdering>       SharedFortranMatrix;
    typedef const FullMatrix<double, ImmutableSharedData, FortranOrdering> ImmutableFortranMatrix;


    template<class Matrix>
    void zero(Matrix& A)
    {
        std::fill_n(A.data(), A.numRows() * A.numCols(),
                    typename Matrix::value_type(0.0));
    }

    template<class Matrix>
    void eye(Matrix& A)
    {
        zero(A);
        for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
            A(i, i) = 1.0;
        }
    }

    template<class Matrix>
    typename Matrix::value_type
    trace(const Matrix& A)
    {
        typename Matrix::value_type ret(0);

        for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
            ret += A(i,i);
        }
        return ret;
    }


    template<class Matrix, int rows>
    Dune::FieldVector<typename Matrix::value_type, rows>
    prod(const Matrix& A, const Dune::FieldVector<typename Matrix::value_type,rows>& x)
    {
        const int cols = rows;
        assert (A.numRows() == rows);
        assert (A.numCols() == cols);

        Dune::FieldVector<typename Matrix::value_type, rows> res(0.0);
        for (int c = 0; c < cols; ++c) {
            for (int r = 0; r < rows; ++r) {
                res[r] += A(r, c)*x[c];
            }
        }
        return res;
    }

    template<class Matrix1, class Matrix2, class MutableMatrix>
    void prod(const Matrix1& A, const Matrix2& B, MutableMatrix& C)
    {
        int result_rows = A.numRows();
        int result_cols = B.numCols();
        int inner_dim = A.numCols();
        assert (inner_dim == B.numRows());
        assert(C.numRows() == result_rows);
        assert(C.numCols() == result_cols);

        for (int c = 0; c < result_cols; ++c) {
            for (int r = 0; r < result_rows; ++r) {
                C(r,c) = 0.0;
                for (int i = 0; i < inner_dim; ++i) {
                    C(r,c) += A(r,i)*B(i,c);
                }
            }
        }
    }


    template<typename T, template<typename> class StoragePolicy>
    int orthogonalizeColumns(FullMatrix<T,StoragePolicy,FortranOrdering>& A)
    {
        typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;

        static std::vector<value_type> tau;
        static std::vector<value_type> work;

        if (int(tau .size()) <      A.numCols()) tau .resize(     A.numCols());
        if (int(work.size()) < 64 * A.numRows()) work.resize(64 * A.numRows());  // 64 from ILAENV

        int info = 0;

        // Generate factorization.  Matrix Q is implicitly represented.
        Opm::BLAS_LAPACK::GEQRF(A.numRows(), A.numCols()               ,
                                 A.data()   , A.leadingDimension()      ,
                                 &tau[0]    , &work[0], int(work.size()),
                                 info);

        if (info == 0)
        {
            // QR factorization successfully generated--extract
            // explict representation of orthogonal matrix Q.
            Opm::BLAS_LAPACK::ORGQR(A.numRows(), A.numCols(), A.numCols()  ,
                                     A.data()   , A.leadingDimension()      ,
                                     &tau[0]    , &work[0], int(work.size()),
                                     info);
        }

        return info;
    }


    template<typename T, template<typename> class StoragePolicy, class OrderingPolicy>
    int invert(FullMatrix<T,StoragePolicy,OrderingPolicy>& A)
    {
        typedef typename FullMatrix<T,StoragePolicy,OrderingPolicy>::value_type value_type;

        assert (A.numRows() == A.numCols());

        std::vector<int> ipiv(A.numRows());
        int info = 0;

        // Correct both for COrdering and FortranOrdering (inv(A)' == inv(A')).
        Opm::BLAS_LAPACK::GETRF(A.numRows(), A.numCols(), A.data(),
                                 A.leadingDimension(), &ipiv[0], info);

        if (info == 0) {
            std::vector<value_type> work(A.numRows());

            Opm::BLAS_LAPACK::GETRI(A.numRows(), A.data(), A.leadingDimension(),
                                     &ipiv[0], &work[0], int(work.size()), info);
        }

        return info;
    }


    template<typename T, template<typename> class StoragePolicy>
    void symmetricUpdate(const T&                                           a1,
                         const FullMatrix<T,StoragePolicy,FortranOrdering>& A ,
                         const T&                                           a2,
                         FullMatrix<T,StoragePolicy,FortranOrdering>&       C )
    {
        Opm::BLAS_LAPACK::SYRK("Upper"     , "No transpose"      ,
                                C.numRows() , A.numCols()         ,
                                a1, A.data(), A.leadingDimension(),
                                a2, C.data(), C.leadingDimension());

        // Account for SYRK (in this case) only updating the upper
        // leading n-by-n submatrix.
        //
        for (int j = 0; j < C.numCols(); ++j) {
            for (int i = j+1; i < C.numRows(); ++i) {
                C(i,j) = C(j,i);
            }
        }
    }


    // B <- A*B*A'
    // Assumes T is an arithmetic (floating point) type, and that A==A'.
    template<typename T, template<typename> class StoragePolicy>
    void symmetricUpdate(const FullMatrix<T,StoragePolicy,FortranOrdering>& A,
                         FullMatrix<T,StoragePolicy,FortranOrdering>&       B)
    {
        typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;

        // B <- A*B
        Opm::BLAS_LAPACK::TRMM("Left" , "Upper", "No transpose", "Non-unit",
                                B.numRows(), B.numCols(), value_type(1.0),
                                A.data(), A.leadingDimension(),
                                B.data(), B.leadingDimension());

        // B <- B*A (== A * B_orig * A)
        Opm::BLAS_LAPACK::TRMM("Right", "Upper", "No transpose", "Non-unit",
                                B.numRows(), B.numCols(), value_type(1.0),
                                A.data(), A.leadingDimension(),
                                B.data(), B.leadingDimension());

        // Account for TRMM (in this case) only updating the upper
        // leading n-by-n submatrix.
        //
        for (int j = 0; j < B.numCols(); ++j) {
            for (int i = j+1; i < B.numRows(); ++i) {
                B(i,j) = B(j,i);
            }
        }
    }


    template<typename                 T ,
             template<typename> class SP>
    void vecMulAdd_N(const T&                                a1,
                     const FullMatrix<T,SP,FortranOrdering>& A ,
                     const std::vector<T>&                   x ,
                     const T&                                a2,
                     std::vector<T>&                         y)
    {
        assert(A.numRows() == y.size());
        assert(A.numCols() == x.size());

        Opm::BLAS_LAPACK::GEMV("No Transpose",
                                A.numRows(), A.numCols(),
                                a1, A.data(), A.leadingDimension(),
                                &x[0], 1, a2, &y[0], 1);
    }


    template<typename                 T ,
             template<typename> class SP>
    void vecMulAdd_N(const T&                                a1,
                     const FullMatrix<T,SP,FortranOrdering>& A ,
                     const T*                                x ,
                     const T&                                a2,
                     T*                                      y)
    {
        Opm::BLAS_LAPACK::GEMV("No Transpose",
                                A.numRows(), A.numCols(),
                                a1, A.data(), A.leadingDimension(),
                                x, 1, a2, y, 1);
    }


    template<typename                 T ,
             template<typename> class SP>
    void vecMulAdd_T(const T&                                a1,
                     const FullMatrix<T,SP,FortranOrdering>& A ,
                     const std::vector<T>&                   x ,
                     const T&                                a2,
                     std::vector<T>&                         y)
    {
        assert (A.numCols() == y.size());
        assert (A.numRows() == x.size());

        Opm::BLAS_LAPACK::GEMV("Transpose",
                                A.numRows(), A.numCols(),
                                a1, A.data(), A.leadingDimension(),
                                &x[0], 1, a2, &y[0], 1);
    }


    template<typename                 T ,
             template<typename> class SP>
    void vecMulAdd_T(const T&                                a1,
                     const FullMatrix<T,SP,FortranOrdering>& A ,
                     const T*                                x ,
                     const T&                                a2,
                     T*                                      y)
    {
        Opm::BLAS_LAPACK::GEMV("Transpose",
                                A.numRows(), A.numCols(),
                                a1, A.data(), A.leadingDimension(),
                                x, 1, a2, y, 1);
    }


    template<typename                 T ,
             template<typename> class SP>
    void vecMulAdd_N(const T&                          a1,
                     const FullMatrix<T,SP,COrdering>& A ,
                     const T*                          x ,
                     const T&                          a2,
                     T*                                y)
    {
        typedef FullMatrix<T, ImmutableSharedData, FortranOrdering> FMAT;

        const FMAT At(A.numCols(), A.numRows(), A.data());

        vecMulAdd_T(a1, At, x, a2, y);
    }


    template<typename                 T  ,
             template<typename> class SP1,
             template<typename> class SP2,
             template<typename> class SP3>
    void matMulAdd_NN(const T&                                 a1,
                      const FullMatrix<T,SP1,FortranOrdering>& A ,
                      const FullMatrix<T,SP2,FortranOrdering>& B ,
                      const T&                                 a2,
                      FullMatrix<T,SP3,FortranOrdering>&       C)
    {
        assert(A.numRows() == C.numRows());
        assert(A.numCols() == B.numRows());
        assert(B.numCols() == C.numCols());

        int m = A.numRows();  // Number of *rows* in A
        int n = B.numCols();  // Number of *cols* in B
        int k = A.numCols();  // Number of *cols* in A (== numer of *rows* in B)

        Opm::BLAS_LAPACK::GEMM("No Transpose", "No Transpose", m, n, k,
                                a1, A.data(), A.leadingDimension(),
                                    B.data(), B.leadingDimension(),
                                a2, C.data(), C.leadingDimension());
    }


    template<typename                 T  ,
             template<typename> class SP1,
             template<typename> class SP2,
             template<typename> class SP3>
    void matMulAdd_NT(const T&                                 a1,
                      const FullMatrix<T,SP1,FortranOrdering>& A ,
                      const FullMatrix<T,SP2,FortranOrdering>& B ,
                      const T&                                 a2,
                      FullMatrix<T,SP3,FortranOrdering>&       C)
    {
        assert(A.numRows() == C.numRows());
        assert(B.numRows() == C.numCols());
        assert(A.numCols() == B.numCols());

        int m = A.numRows();  // Number of *rows* in A
        int n = B.numRows();  // Number of *cols* in B'
        int k = A.numCols();  // Number of *cols* in A (== numer of *rows* in B')

        Opm::BLAS_LAPACK::GEMM("No Transpose", "Transpose", m, n, k,
                                a1, A.data(), A.leadingDimension(),
                                    B.data(), B.leadingDimension(),
                                a2, C.data(), C.leadingDimension());
    }


    template<typename                 T  ,
             template<typename> class SP1,
             template<typename> class SP2,
             template<typename> class SP3>
    void matMulAdd_TN(const T&                                 a1,
                      const FullMatrix<T,SP1,FortranOrdering>& A ,
                      const FullMatrix<T,SP2,FortranOrdering>& B ,
                      const T&                                 a2,
                      FullMatrix<T,SP3,FortranOrdering>&       C)
    {
        assert (A.numCols() == C.numRows());
        assert (A.numRows() == B.numRows());
        assert (B.numCols() == C.numCols());

        int m = A.numCols();  // Number of *rows* in A'
        int n = B.numCols();  // Number of *cols* in B
        int k = A.numRows();  // Number of *cols* in A' (== numer of *rows* in B)

        Opm::BLAS_LAPACK::GEMM("Transpose", "No Transpose", m, n, k,
                                a1, A.data(), A.leadingDimension(),
                                    B.data(), B.leadingDimension(),
                                a2, C.data(), C.leadingDimension());
    }


    template<typename                 T  ,
             template<typename> class SP1,
             template<typename> class SP2,
             template<typename> class SP3>
    void matMulAdd_NN(const T&                                 a1,
                      const FullMatrix<T,SP1,FortranOrdering>& A ,
                      const FullMatrix<T,SP2,COrdering>&       B ,
                      const T&                                 a2,
                      FullMatrix<T,SP3,FortranOrdering>&       C)
    {
        typedef typename FullMatrix<T,SP2,COrdering>::value_type           value_type;
        typedef FullMatrix<value_type,ImmutableSharedData,FortranOrdering> FMat;

        const FMat Bt(B.numCols(), B.numRows(), B.data());

        matMulAdd_NT(a1, A, Bt, a2, C);
    }


    template<class charT, class traits,
             typename T, template<typename> class SP, class OP>
    std::basic_ostream<charT,traits>&
    operator<<(std::basic_ostream<charT,traits>& os,
               const FullMatrix<T,SP,OP>& A)
    {
        for (int i = 0; i < A.numRows(); ++i) {
            for (int j = 0; j < A.numCols(); ++j)
                os << A(i,j) << ' ';
            os << '\n';
        }

        return os;
    }
} // namespace Opm
#endif // OPENRS_MATRIX_HEADER