diagonal_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 diagonal_matrix_h
#define diagonal_matrix_h
 
#include "matrix.h"
#include "types.h"
#include "_memory.h"
 
namespace fmatvec {
 
  template <class AT> class Matrix<Diagonal,Ref,Ref,AT> {
 
    protected:
 
 
      Memory<AT> memory;
      AT *ele;
      int n{0};
 
      inline Matrix<Diagonal,Ref,Ref,AT>& copy(const Matrix<Diagonal,Ref,Ref,AT> &A);
 
      const AT* elePtr(int i) const {
        return ele+i;
      }
 
      AT* elePtr(int i) {
        return ele+i;
      }
 
 
    public:
      static constexpr bool isVector {false};
 
      using value_type = AT;
      using shape_type = Diagonal;
 
      explicit Matrix() : memory(), ele(0) { }
 
      //        template<class Ini=All<AT>>
      //          Matrix(int n_, Ini ini=All<AT>()) : memory(n_), ele((AT*)memory.get()), n(n_) { init(ini); }
      //        template<class Ini=All<AT>>
      //          Matrix(int m_, int n_, Ini ini=All<AT>()) : memory(n_), ele((AT*)memory.get()), n(n_) { init(ini); }
 
      explicit Matrix(int n_, Noinit) : memory(n_), ele((AT*)memory.get()), n(n_) { }
      explicit Matrix(int n_, Init ini=INIT, const AT &a=AT()) : memory(n_), ele((AT*)memory.get()), n(n_) { init(a); }
      explicit Matrix(int n_, Eye ini, const AT &a=1) : memory(n_), ele((AT*)memory.get()), n(n_) { init(ini,a); }
      explicit Matrix(int m_, int n_, Noinit) : memory(n_), ele((AT*)memory.get()), n(n_) { }
 
      Matrix(const Matrix<Diagonal,Ref,Ref,AT> &A) : memory(A.n), ele((AT*)memory.get()), n(A.n) {
        copy(A);
      }
 
      template<class Type, class Row, class Col>
        Matrix(const Matrix<Type,Row,Col,AT> &A) : memory(A.rows()), ele((AT*)memory.get()), n(A.rows()) {
          FMATVEC_ASSERT(A.rows() == A.cols(), AT);
          copy(A);
        }
 
      ~Matrix() = default;
 
      Matrix<Diagonal,Ref,Ref,AT>& resize(int n_, Noinit) {
        n = n_;
        memory.resize(n);
        ele = (AT*)memory.get();
        return *this;
      }
 
      Matrix<Diagonal,Ref,Ref,AT>& resize(int n, Init ini=INIT, const AT &a=AT()) { resize(n,Noinit()).init(a); }
      Matrix<Diagonal,Ref,Ref,AT>& resize(int n, Eye ini, const AT &a=1) { resize(n,Noinit()).init(ini,a); }
 
      inline Matrix<Diagonal,Ref,Ref,AT>& operator=(const Matrix<Diagonal,Ref,Ref,AT> &A) {
        FMATVEC_ASSERT(n == A.n, AT);
        return copy(A);
      }
 
      template<class Type, class Row, class Col>
        inline Matrix<Diagonal,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
          FMATVEC_ASSERT(n == A.rows(), AT);
          FMATVEC_ASSERT(n == A.cols(), AT);
          return copy(A);
        }
 
      inline Matrix<Diagonal,Ref,Ref,AT>& operator&=(Matrix<Diagonal,Ref,Ref,AT> &A) {
        n=A.n;
        memory = A.memory;
        ele = A.ele;
        return *this;
      }
 
      template<class Type, class Row, class Col>
        inline Matrix<Diagonal,Ref,Ref,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
          if(n!=A.n) resize(A.rows(),NONINIT);
          return copy(A);
        }
 
      const AT& operator()(int i, int j) const {
 
        FMATVEC_ASSERT(i>=0, AT);
        FMATVEC_ASSERT(j>=0, AT);
        FMATVEC_ASSERT(i<n, AT);
        FMATVEC_ASSERT(j<n, AT);
 
        return e(i,j);
      }
 
      const AT& operator()(int i) const {
 
        FMATVEC_ASSERT(i>=0, AT);
        FMATVEC_ASSERT(i<n, AT);
 
        return e(i);
      }
 
      AT& operator()(int i) {
 
        FMATVEC_ASSERT(i>=0, AT);
        FMATVEC_ASSERT(i<n, AT);
 
        return e(i);
      }
 
      const AT& e(int i, int j) const {
        static AT zero=0;
        return i==j ? ele[i] : zero;
      }
 
      const AT& e(int i) const {
        return ele[i];
      }
 
      AT& e(int i) {
        return ele[i];
      }
 
      const AT* operator()() const {
        return ele;
      }
 
      AT* operator()() {
        return ele;
      }
 
      int rows() const {return n;}
 
      int cols() const {return n;}
 
      int size() const {return n;}
 
      CBLAS_ORDER blasOrder() const {
        return CblasColMajor;
      }
 
      Matrix<Diagonal,Ref,Ref,AT>& init(const AT &val=AT());
      inline Matrix<Diagonal,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
      inline Matrix<Diagonal,Ref,Ref,AT>& init(Eye, const AT &a=1) { return init(a); }
      inline Matrix<Diagonal,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
 
      explicit operator std::vector<std::vector<AT>>() const;
 
      explicit 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(n==1, AT);
//        return ele[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()), n(1) { ele[0] = x; }
  };
 
  template <class AT>
    Matrix<Diagonal,Ref,Ref,AT>::operator std::vector<std::vector<AT>>() const {
      std::vector<std::vector<AT>> ret(rows(),std::vector<AT>(cols(),AT()));
      for(int r=0; r<rows(); r++) {
    ret[r][r]=e(r);
      }
      return ret;
    }
 
  template <class AT>
    Matrix<Diagonal,Ref,Ref,AT>::Matrix(const std::vector<std::vector<AT>> &m) :  memory(m.size()), ele((AT*)memory.get()), n(static_cast<int>(m.size())) {
      for(size_t r=0; r<m.size(); ++r) {
        if(n!=static_cast<int>(m[r].size()))
          throw std::runtime_error("The rows of the input have different lenght.");
        for(size_t c=0; c<m[r].size(); ++c) {
          if(r==c)
            e(r)=m[r][r];
          if(r!=c && (fabs(m[r][c])>1e-13 || fabs(m[r][c])>1e-13))
            throw std::runtime_error("The input is not diagonal.");
        }
      }
    }
 
  template <class AT>
    Matrix<Diagonal,Ref,Ref,AT>&  Matrix<Diagonal,Ref,Ref,AT>::init(const AT &val) {
      for(int i=0; i<rows(); i++) 
        e(i) = val;
      return *this;
    }
 
  template <class AT>
    inline Matrix<Diagonal,Ref,Ref,AT>& Matrix<Diagonal,Ref,Ref,AT>::copy(const Matrix<Diagonal,Ref,Ref,AT> &A) {
      for(int i=0; i<n; i++)
        e(i) = A.e(i);
      return *this;
    }
 
}
 
#endif