general_matrix.h

/* Copyright (C) 2003-2005  Martin Förg
 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library 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
 * Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact:
 *   martin.o.foerg@googlemail.com
 *
 */
 
#ifndef general_matrix_h
#define general_matrix_h
 
#include "_memory.h"
#include "types.h"
#include "indices.h"
#include "matrix.h"
#include <cstdlib>
#include <stdexcept>
 
namespace fmatvec {
 
  template <class AT> class Matrix<General,Ref,Ref,AT> {
 
    public:
      static constexpr bool isVector {false};
 
      using value_type = AT;
      using shape_type = General;
 
 
      friend class Vector<Ref,AT>;
      friend class RowVector<Ref,AT>;
      friend class SquareMatrix<Ref,AT>;
      friend class Matrix<Symmetric,Ref,Ref,AT>;
      
    protected:
 
      template<bool Const>
      class Iterator {
        using Mat = Matrix<General, Ref, Ref, AT>;
        public:
          Iterator(std::conditional_t<Const, const Mat, Mat> &mat_, int row_, int col_) : mat(mat_), row(row_), col(col_) {}
          bool operator!=(const Iterator& x) const { return row!=x.row || col!=x.col; }
          Iterator& operator++() { ++row; if(row>=mat.m) { row=0; ++col; } return *this; }
          std::conditional_t<Const, const AT, AT>& operator*() { return mat.e(row,col); }
          std::conditional_t<Const, const AT, AT>& operator->() { return mat.e(row,col); }
          const AT& operator*() const { return mat.e(row,col); }
          const AT& operator->() const { return mat.e(row,col); }
        private:
          std::conditional_t<Const, const Mat, Mat> &mat;
          int row;
          int col;
      };
 
      Memory<AT> memory;
      AT *ele;
      int m{0};
      int n{0};
      int lda{0};
 
      template <class Type, class Row, class Col> inline Matrix<General,Ref,Ref,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
      inline Matrix<General,Ref,Ref,AT>& copy(const Matrix<Sparse,Ref,Ref,AT> &A);
      inline Matrix<General,Ref,Ref,AT>& copy(const Matrix<SymmetricSparse,Ref,Ref,AT> &A);
 
      const AT* elePtr(int i, int j) const {
    return ele+i+j*lda;
      }
 
      AT* elePtr(int i, int j) {
    return ele+i+j*lda;
      }
 
 
    public:
 
      explicit Matrix() : memory(), ele(nullptr) { }
 
// Works with -std=gnu++0x only
//      template<class Ini=All<AT>>
//      Matrix(int m_, int n_, Ini ini=All<AT>()) : memory(m_*n_), ele((AT*)memory.get()), m(m_), n(n_), lda(m_) {
//        init(ini);
//      }
 
      explicit Matrix(int m_, int n_, Noinit) : memory(m_*n_), ele((AT*)memory.get()), m(m_), n(n_), lda(m_) { }
      explicit Matrix(int m_, int n_, Init ini=INIT, const AT &a=AT()) : memory(m_*n_), ele((AT*)memory.get()), m(m_), n(n_), lda(m_) { init0(a); }
      explicit Matrix(int m_, int n_, Eye ini, const AT &a=1) : memory(m_*n_), ele((AT*)memory.get()), m(m_), n(n_), lda(m_) { init(ini,a); }
 
      Matrix(const Matrix<General,Ref,Ref,AT> &A) : memory(A.m*A.n), ele((AT*)memory.get()) ,m(A.m), n(A.n), lda(m) {
        copy(A);
      }
 
      template<class Row, class Col>
      Matrix(const Matrix<General,Row,Col,AT> &A) : memory(A.rows()*A.cols()), ele((AT*)memory.get()), m(A.rows()), n(A.cols()), lda(m) {
        copy(A);
      }
 
      template<class Type, class Row, class Col>
      explicit Matrix(const Matrix<Type,Row,Col,AT> &A) : memory(A.rows()*A.cols()), ele((AT*)memory.get()), m(A.rows()), n(A.cols()), lda(m) {
    copy(A);
      }
 
      explicit Matrix(int m_, int n_, AT* ele_) : memory(), ele(ele_), m(m_), n(n_), lda(m_) { 
      }
 
      Matrix(const std::string &strs);
      Matrix(const char *strs);
 
      ~Matrix() = default;
 
      using iterator = Iterator<false>;
      using const_iterator = Iterator<true>;
      iterator begin() { return iterator(*this,0,0); }
      iterator end() { return iterator(*this,0,n); }
      const_iterator begin() const { return const_iterator(*this,0,0); }
      const_iterator end() const { return const_iterator(*this,0,n); }
 
      Matrix<General,Ref,Ref,AT>& resize(int m_, int n_, Noinit) {
    m = m_; n = n_; lda = m;
    memory.resize(m*n);
    ele = (AT*)memory.get();
        return *this;
      }
 
      Matrix<General,Ref,Ref,AT>& resize(int m, int n, Init ini=INIT, const AT &a=AT()) { return resize(m,n,Noinit()).init0(a); }
 
      Matrix<General,Ref,Ref,AT>& resize(int m, int n, Eye ini, const AT &a=1) { return resize(m,n,Noinit()).init(ini,a); } 
 
      Matrix<General,Ref,Ref,AT>& operator=(const Matrix<General,Ref,Ref,AT> &A) {
        FMATVEC_ASSERT(m == A.rows(), AT);
        FMATVEC_ASSERT(n == A.cols(), AT);
        return copy(A);
      }
 
      template<class Type, class Row, class Col>
      Matrix<General,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
        FMATVEC_ASSERT(m == A.rows(), AT);
        FMATVEC_ASSERT(n == A.cols(), AT);
        return copy(A);
      }
 
      inline Matrix<General,Ref,Ref,AT>& operator&=(Matrix<General,Ref,Ref,AT> &A) {
        m=A.m; 
        n=A.n;
        memory = A.memory;
        ele = A.ele;
        lda = A.lda;
        return *this;
      }
 
      template<class Type, class Row, class Col>
        inline Matrix<General,Ref,Ref,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
          if(m!=A.rows() || n!=A.cols()) resize(A.rows(),A.cols(),NONINIT);
          return copy(A);
        }
 
      template <class AT2>
      operator Matrix<General,Ref,Ref,AT2>() const {
        Matrix<General,Ref,Ref,AT2> ret(rows(), cols());
        for(size_t r=0; r<rows(); ++r)
          for(size_t c=0; c<cols(); ++c)
            ret(r,c) = (*this)(r,c);
        return ret;
      }
 
      AT& operator()(int i, int j) {
    FMATVEC_ASSERT(i>=0, AT);
    FMATVEC_ASSERT(j>=0, AT);
    FMATVEC_ASSERT(i<m, AT);
    FMATVEC_ASSERT(j<n, AT);
 
    return e(i,j);
      }
 
      const AT& operator()(int i, int j) const {
    FMATVEC_ASSERT(i>=0, AT);
    FMATVEC_ASSERT(j>=0, AT);
    FMATVEC_ASSERT(i<m, AT);
    FMATVEC_ASSERT(j<n, AT);
 
    return e(i,j);
      }
 
      AT& e(int i, int j) {
        return ele[i+j*lda];
      }
 
      const AT& e(int i, int j) const {
    return ele[i+j*lda];
      }
 
      AT* operator()() {return ele;}
 
      const AT* operator()() const {return ele;}
 
      int rows() const {return m;}
 
      int cols() const {return n;}
 
      int ldim() const {return lda;}
 
      CBLAS_TRANSPOSE blasTrans() const {
        return CblasNoTrans;
      }
 
      CBLAS_ORDER blasOrder() const {
        return CblasColMajor;
      }
 
      inline const Matrix<General,Ref,Ref,AT> operator()(const Range<Var,Var> &I, const Range<Var,Var> &J) const;
 
      inline const SquareMatrix<Ref,AT> operator()(const Range<Var,Var> &I) const;
 
      inline const RowVector<Ref,AT> row(int i) const;
 
      inline const Vector<Ref,AT> col(int j) const;
 
      inline Matrix<General,Ref,Ref,AT>& init(const AT &val=AT());
      inline Matrix<General,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
      inline Matrix<General,Ref,Ref,AT>& init(Eye, const AT &val=1);
      inline Matrix<General,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
      inline Matrix<General,Ref,Ref,AT>& init0(const AT &val=AT());
      inline Matrix<General,Ref,Ref,AT>& init0(Init, const AT &a=AT()) { return init0(a); }
 
      explicit inline operator std::vector<std::vector<AT>>() const;
 
      explicit inline Matrix(const std::vector<std::vector<AT>> &m);
 
//      /*! \brief Cast to AT.
//       *
//       * \return The AT representation of the matrix
//       * */
//      explicit operator AT() const {
//        FMATVEC_ASSERT(m==1, AT);
//        FMATVEC_ASSERT(n==1, AT);
//        return e(0);
//      }
//
//      /*! \brief AT Constructor.
//       * Constructs and initializes a matrix with a AT object.
//       * \param x The AT the matrix will be initialized with.
//       * */
//      explicit Matrix(const AT &x) : memory(1), ele((AT*)memory.get()), m(1), n(1), lda(1) { ele[0] = x; }
 
      const Matrix<General,Ref,Ref,AT> T() const;
 
      template<class Row> inline void set(int j, const Vector<Row,AT> &x);
 
      template<class Col> inline void set(int i, const RowVector<Col,AT> &x);
 
      template<class Type, class Row, class Col> inline void set(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A);
 
      template<class Row> inline void add(int j, const Vector<Row,AT> &x);
 
      template<class Col> inline void add(int i, const RowVector<Col,AT> &x);
 
      template<class Type, class Row, class Col> inline void add(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A);
 
      template<class Type, class Row, class Col> inline void sub(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A);
 
      inline const Matrix<General,Ref,Ref,AT> operator()(const Indices &I, const Indices &J) const;
 
      template<class Type, class Row, class Col> inline void set(const Indices &I, const Indices &J, const Matrix<Type,Row,Col,AT> &A);
 
      template<class Row> inline void set(const Indices &I, int j, const Vector<Row,AT> &x);
 
      inline void ref(Matrix<General,Ref,Ref,AT> &A, const Range<Var,Var> &I, const Range<Var,Var> &J);
 
      inline void ref(Matrix<Symmetric,Ref,Ref,AT> &A, const Range<Var,Var> &I, const Range<Var,Var> &J);
 
      int nonZeroElements() const;
  };
 
  template <class AT> 
    Matrix<General,Ref,Ref,AT>::Matrix(const std::string &strs) : memory(), ele(nullptr) {
      std::istringstream iss(strs);
      iss>>*this;
 
      // check end of stream
      iss>>std::ws;
      if(!iss.eof())
        throw std::runtime_error("Input not fully read.");
    }
  template <class AT> Matrix<General,Ref,Ref,AT>::Matrix(const char *strs) : Matrix<General,Ref,Ref,AT>::Matrix(std::string(strs)) {}
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>&  Matrix<General,Ref,Ref,AT>::init0(const AT &val) {
      for(int i=0; i<m*n; i++) ele[i]=val;
      return *this;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>&  Matrix<General,Ref,Ref,AT>::init(const AT &val) {
      for(int i=0; i<rows(); i++)
        for(int j=0; j<cols(); j++)
          e(i,j) = val;
      return *this;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>&  Matrix<General,Ref,Ref,AT>::init(Eye, const AT &val) {
      for(int i=0; i<m; i++)
        for(int j=0; j<n; j++)
          e(i,j) = (i==j) ? val : 0;
      return *this;
    }
 
  template <class AT> 
    inline const Matrix<General,Ref,Ref,AT> Matrix<General,Ref,Ref,AT>::operator()(const Range<Var,Var> &I, const Range<Var,Var> &J) const {
      FMATVEC_ASSERT(I.end()<m, AT);
      FMATVEC_ASSERT(J.end()<n, AT);
      Matrix<General,Ref,Ref,AT> A(I.size(),J.size(),NONINIT);
 
      for(int i=0; i<A.rows(); i++)
        for(int j=0; j<A.cols(); j++)
          A.e(i,j) = e(I.start()+i,J.start()+j);
 
      return A;
    }
 
  template <class AT> 
    inline const SquareMatrix<Ref,AT> Matrix<General,Ref,Ref,AT>::operator()(const Range<Var,Var> &I) const {
      FMATVEC_ASSERT(I.end()<m, AT);
      SquareMatrix<Ref,AT> A(I.size(),NONINIT);
 
      for(int i=0; i<A.rows(); i++)
        for(int j=0; j<A.cols(); j++)
          A.e(i,j) = e(I.start()+i,I.start()+j);
 
      return A;
    }
 
  template <class AT>
    inline const RowVector<Ref,AT>  Matrix<General,Ref,Ref,AT>::row(int i) const {
 
      FMATVEC_ASSERT(i>=0, AT);
      FMATVEC_ASSERT(i<m, AT);
 
      RowVector<Ref,AT> x(n,NONINIT);
 
      for(int j=0; j<n; j++)
        x.e(j) = e(i,j);
 
      return x;
    }
 
  template <class AT>
    inline const Vector<Ref,AT>  Matrix<General,Ref,Ref,AT>::col(int j) const {
 
      FMATVEC_ASSERT(j>=0, AT);
      FMATVEC_ASSERT(j<n, AT);
 
      Vector<Ref,AT> x(m,NONINIT);
 
      for(int i=0; i<m; i++)
        x.e(i) = e(i,j);
 
      return x;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>::operator std::vector<std::vector<AT>>() const {
      std::vector<std::vector<AT>> ret(rows(),std::vector<AT>(cols()));
      for(int r=0; r<rows(); r++) {
        for(int c=0; c<cols(); c++)
          ret[r][c]= e(r,c);
      }
      return ret;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>::Matrix(const std::vector<std::vector<AT>> &m) : memory(m.size()>0?m.size()*(m.empty()?0:m[0].size()):0), ele((AT*)memory.get()), m(static_cast<int>(m.size())), n(static_cast<int>(m.size()>0?(m.empty()?0:m[0].size()):0)), lda(static_cast<int>(m.size())) {
      for(int r=0; r<rows(); r++) {
        if(static_cast<int>(m[r].size())!=cols())
          throw std::runtime_error("The rows of the input have different length.");
        for(int c=0; c<cols(); c++)
          e(r,c)=m[r][c];
      }
    }
 
 
  template <class AT> template <class Type, class Row, class Col>
    inline Matrix<General,Ref,Ref,AT>& Matrix<General,Ref,Ref,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
      for(int i=0; i<m; i++)
        for(int j=0; j<n; j++)
          e(i,j) = A.e(i,j);
      return *this;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>& Matrix<General,Ref,Ref,AT>::copy(const Matrix<Sparse,Ref,Ref,AT> &A) {
      for(int i=0; i<m; i++) {
    for(int j=0; j<n; j++)
      e(i,j) = 0;
        for(int j=A.Ip()[i]; j<A.Ip()[i+1]; j++)
          e(i,A.Jp()[j]) = A()[j];
      }
      return *this;
    }
 
  template <class AT>
    inline Matrix<General,Ref,Ref,AT>& Matrix<General,Ref,Ref,AT>::copy(const Matrix<SymmetricSparse,Ref,Ref,AT> &A) {
    for(int i=0; i<m; i++) {
      for(int j=0; j<n; j++)
    e(i,j) = 0;
      e(i,A.Jp()[A.Ip()[i]]) = A()[A.Ip()[i]];
    }
    for(int i=0; i<m; i++) {
      for(int j=A.Ip()[i]+1; j<A.Ip()[i+1]; j++) {
    e(i,A.Jp()[j]) = A()[j];
    e(A.Jp()[j],i) = A()[j];
      }
    }
    return *this;
  }
 
  template <class AT>
    inline const Matrix<General,Ref,Ref,AT> Matrix<General,Ref,Ref,AT>::T() const {
      Matrix<General,Ref,Ref,AT> A(n,m,NONINIT);
      for(int i=0; i<n; i++)
        for(int j=0; j<m; j++)
          A.e(i,j) = e(j,i);
      return A;
    }
 
  template <class AT> template <class Row>
    inline void Matrix<General,Ref,Ref,AT>::set(int j, const Vector<Row,AT> &x) {
      FMATVEC_ASSERT(j<cols(), AT);
      FMATVEC_ASSERT(rows()==x.size(), AT);
      for(int i=0; i<rows(); i++)
        e(i,j) = x.e(i);
    }
 
  template <class AT> template <class Col>
    inline void Matrix<General,Ref,Ref,AT>::set(int i, const RowVector<Col,AT> &x) {
      FMATVEC_ASSERT(i<rows(), AT);
      FMATVEC_ASSERT(cols()==x.size(), AT);
      for(int j=0; j<cols(); j++)
        e(i,j) = x.e(j);
    }
 
  template <class AT> template<class Type, class Row, class Col>
    inline void Matrix<General,Ref,Ref,AT>::set(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A) {
 
      FMATVEC_ASSERT(I.end()<rows(), AT);
      FMATVEC_ASSERT(J.end()<cols(), AT);
      FMATVEC_ASSERT(I.size()==A.rows(), AT);
      FMATVEC_ASSERT(J.size()==A.cols(), AT);
 
      for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
        for(int j=J.start(), jj=0; j<=J.end(); j++, jj++)
          e(i,j) = A.e(ii,jj);
    }
 
  template <class AT> template <class Row>
    inline void Matrix<General,Ref,Ref,AT>::add(int j, const Vector<Row,AT> &x) {
      FMATVEC_ASSERT(j<cols(), AT);
      FMATVEC_ASSERT(rows()==x.size(), AT);
      for(int i=0; i<rows(); i++)
        e(i,j) += x.e(i);
    }
 
  template <class AT> template <class Col>
    inline void Matrix<General,Ref,Ref,AT>::add(int i, const RowVector<Col,AT> &x) {
      FMATVEC_ASSERT(i<rows(), AT);
      FMATVEC_ASSERT(cols()==x.size(), AT);
      for(int j=0; j<cols(); j++)
        e(i,j) += x.e(j);
    }
 
  template <class AT> template<class Type, class Row, class Col>
    inline void Matrix<General,Ref,Ref,AT>::add(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A) {
 
      FMATVEC_ASSERT(I.end()<rows(), AT);
      FMATVEC_ASSERT(J.end()<cols(), AT);
      FMATVEC_ASSERT(I.size()==A.rows(), AT);
      FMATVEC_ASSERT(J.size()==A.cols(), AT);
 
      for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
        for(int j=J.start(), jj=0; j<=J.end(); j++, jj++)
          e(i,j) += A.e(ii,jj);
    }
 
  template <class AT> template<class Type, class Row, class Col>
    inline void Matrix<General,Ref,Ref,AT>::sub(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A) {
 
      FMATVEC_ASSERT(I.end()<rows(), AT);
      FMATVEC_ASSERT(J.end()<cols(), AT);
      FMATVEC_ASSERT(I.size()==A.rows(), AT);
      FMATVEC_ASSERT(J.size()==A.cols(), AT);
 
      for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
        for(int j=J.start(), jj=0; j<=J.end(); j++, jj++)
          e(i,j) -= A.e(ii,jj);
    }
 
  template <class AT>
    inline const Matrix<General,Ref,Ref,AT> Matrix<General,Ref,Ref,AT>::operator()(const Indices &I, const Indices &J) const {
      FMATVEC_ASSERT(I.max()<rows(), AT);
      FMATVEC_ASSERT(J.max()<cols(), AT);
 
      Matrix<General,Ref,Ref,AT> A(I.size(),J.size(),NONINIT);
 
      for(int i=0; i<A.rows(); i++)
        for(int j=0; j<A.cols(); j++)
      A.e(i,j) = e(I[i],J[j]);
 
      return A;
    }
 
  template <class AT> template <class Type, class Row, class Col>
    inline void Matrix<General,Ref,Ref,AT>::set(const Indices &I, const Indices &J, const Matrix<Type,Row,Col,AT> &A) {
      FMATVEC_ASSERT(I.max()<rows(), AT);
      FMATVEC_ASSERT(J.max()<cols(), AT);
      FMATVEC_ASSERT(I.size()==A.rows(), AT);
      FMATVEC_ASSERT(J.size()==A.cols(), AT);
      for(int i=0; i<I.size(); i++)
    for(int j=0; j<J.size(); j++)
      e(I[i],J[j]) = A.e(i,j);
    }
 
  template <class AT> template<class Row>
    inline void Matrix<General,Ref,Ref,AT>::set(const Indices &I, int j, const Vector<Row,AT> &x) {
      FMATVEC_ASSERT(I.max()<rows(), AT);
      FMATVEC_ASSERT(j<cols(), AT);
      FMATVEC_ASSERT(I.size()==x.size(), AT);
      for(int i=0; i<I.size(); i++)
    e(I[i],j) = x.e(i);
    }
 
  template <class AT>
    inline void Matrix<General,Ref,Ref,AT>::ref(Matrix<General,Ref,Ref,AT> &A, const Range<Var,Var> &I, const Range<Var,Var> &J) {
      FMATVEC_ASSERT(I.end()<A.rows(), AT);
      FMATVEC_ASSERT(J.end()<A.cols(), AT);
      m=I.size();
      n=J.size();
      memory = A.memory;
      ele = A.elePtr(I.start(),J.start());
      lda = A.lda;
    }
 
  template <class AT>
    inline void Matrix<General,Ref,Ref,AT>::ref(Matrix<Symmetric,Ref,Ref,AT> &A, const Range<Var,Var> &I, const Range<Var,Var> &J) {
      FMATVEC_ASSERT(I.end()<A.size(), AT);
      FMATVEC_ASSERT(J.end()<A.size(), AT);
      FMATVEC_ASSERT(I.start()>=J.start(), AT);
      FMATVEC_ASSERT(J.end()<=I.start(), AT);
      m=I.size();
      n=J.size();
      memory = A.memory;
      ele = A.elePtr(I.start(),J.start());
      lda = A.lda;
    }
 
  template <class AT>
    int Matrix<General,Ref,Ref,AT>::nonZeroElements() const {
    int k = 0;
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (fabs(e(i, j)) > 1e-16 || i==j) {
          k++;
        }
      }
    }
    return k;
  }
 
 
}
 
#endif