band_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_band_matrix_h
#define general_band_matrix_h
 
#include "vector.h"
#include "types.h"
#include "_memory.h"
#include <cstdlib>
 
namespace fmatvec {
 
  template <class AT> class Matrix<GeneralBand,Ref,Ref,AT> {
 
 
    protected:
 
      Memory<AT> memory;
      AT *ele;
      int n{0};
      int kl{0};
      int ku{0};
 
      inline Matrix<GeneralBand,Ref,Ref,AT>& copy(const Matrix<GeneralBand,Ref,Ref,AT> &A);
 
 
    public:
      static constexpr bool isVector {false};
 
      using value_type = AT;
      using shape_type = GeneralBand;
 
      explicit Matrix() : memory(), ele(0) {
      }
 
//      template<class Ini=All<AT>>
//        Matrix(int n_, int kl_, int ku_, Ini ini=All<AT>()) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) {  
//          init(ini);
//        }
//      template<class Ini=All<AT>>
//        Matrix(int m_, int n_, Ini ini=All<AT>()) : memory(n_*(n_+n_+1)), ele((AT*)memory.get()), n(n_), kl(n_), ku(n_) {  
//          init(ini);
//        }
 
      explicit Matrix(int n_, int kl_, int ku_, Noinit) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { }
      explicit Matrix(int n_, int kl_, int ku_, Init ini=INIT, const AT &a=AT()) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { init(a); }
      explicit Matrix(int n_, int kl_, int ku_, Eye ini, const AT &a=1) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { init(ini,a); }
      explicit Matrix(int m_, int n_, int kl_, int ku_, Noinit) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { }
      explicit Matrix(int m_, int n_, int kl_, int ku_, Init ini=INIT, const AT &a=AT()) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { init(a); }
      explicit Matrix(int m_, int n_, int kl_, int ku_, Eye ini, const AT &a=1) : memory(n_*(kl_+ku_+1)), ele((AT*)memory.get()), n(n_), kl(kl_), ku(ku_) { init(ini,a); }
 
      Matrix(const Matrix<GeneralBand,Ref,Ref,AT> &A) : memory(A.n*(A.kl+A.ku+1)), ele((AT*)memory.get()), n(A.n), kl(A.kl), ku(A.ku) {
        copy(A);
      }
 
      template<class Type, class Row, class Col>
      explicit Matrix(const Matrix<Type,Row,Col,AT> &A) : memory(A.rows()*(2*A.rows()-1)), ele((AT*)memory.get()), n(A.rows()), kl(n-1), ku(n-1) {
        FMATVEC_ASSERT(A.rows() == A.cols(), AT);
    copy(A);
      }
 
      explicit Matrix(int n_, int kl_, int ku_, AT* ele_) : memory(), ele(ele_), n(n_), kl(kl_), ku(ku_) { 
      }
 
      ~Matrix() = default;
 
      Matrix<GeneralBand,Ref,Ref,AT>& resize(int n_, int kl_, int ku_, Noinit) { 
          n = n_; kl = kl_; ku = ku_;
          memory.resize(n*(kl+ku+1));
          ele = (AT*)memory.get();
          return *this;
      }
 
      Matrix<GeneralBand,Ref,Ref,AT>& resize(int n, int kl, int ku=0, Init ini=INIT, const AT &a=AT()) { return resize(n,kl,ku,Noinit()).init(a); }
 
      Matrix<GeneralBand,Ref,Ref,AT>& resize(int n, int kl, int ku, Eye ini, const AT &a=1) { return resize(n,kl,ku,Noinit()).init(a); }
 
      inline Matrix<GeneralBand,Ref,Ref,AT>& operator=(const Matrix<GeneralBand,Ref,Ref,AT> &A) {
        FMATVEC_ASSERT(n == A.n, AT);
        return copy(A);
      }
 
      template<class Type, class Row, class Col>
      inline Matrix<GeneralBand,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<GeneralBand,Ref,Ref,AT>& operator&=(Matrix<GeneralBand,Ref,Ref,AT> &A) {
        n = A.n;
        memory = A.memory;
        ele = A.ele;
        return *this;
      }
 
      template<class Type, class Row, class Col>
      inline Matrix<GeneralBand,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);
    //     FMATVEC_ASSERT(i-j<=kl, AT);
    //      FMATVEC_ASSERT(i-j>=-ku, AT);
    static AT zero=0;
 
    return ((i-j>kl) || (i-j<-ku)) ? zero : ele[ku+i+j*(kl+ku)];
      }
 
      AT* operator()() {return ele;}
 
      const AT* operator()() const {return ele;}
 
      int size() const {return n;}
 
      int rows() const {return n;}
 
      int cols() const {return n;}
 
      CBLAS_ORDER blasOrder() const {
        return CblasColMajor;
      }
 
//      /*! \brief Diagonal operator
//       *
//       * See operator()(int)
//       * */
//      inline const Vector<Ref,AT> operator()(int i) const;
//
//      /*! \brief Diagonal operator
//       *
//       * If i<0, the i-th subdiagonal is returned. Otherwise the i-th
//       * superdiagonal is returned.
//       * \param i The i-th super- and subdiagonal,
//       * respectively.
//       * */
//      inline Vector<Ref,AT> operator()(int i);
 
      inline Matrix<GeneralBand,Ref,Ref,AT>& init(const AT &val=AT());
      inline Matrix<GeneralBand,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
      inline Matrix<GeneralBand,Ref,Ref,AT>& init(Eye, const AT &val=1);
      inline Matrix<GeneralBand,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
 
      explicit inline operator std::vector<std::vector<AT>>() const;
  };
 
  template <class AT>
    inline Matrix<GeneralBand,Ref,Ref,AT>&  Matrix<GeneralBand,Ref,Ref,AT>::init(const AT& val) {
      for(int i=0; i<kl+ku+1; i++)
    for(int j=0; j<n; j++)
      ele[i+j*(kl+ku+1)] = val;
      return *this;
    }
 
  template <class AT>
    inline Matrix<GeneralBand,Ref,Ref,AT>&  Matrix<GeneralBand,Ref,Ref,AT>::init(Eye, const AT &val) {
      for(int i=0; i<n; i++)
        ele[i] = val;
      return *this;
    }
 
//  template <class AT>
//    inline Vector<Ref,AT> Matrix<GeneralBand,Ref,Ref,AT>::operator()(int i) {
//      FMATVEC_ASSERT(i<=ku, AT);
//      FMATVEC_ASSERT(i>=-kl, AT);
//
//      //return Vector<Ref,AT>(n-i,ku+kl+1,memory,ele[ku-i + (i>0?i:0)*(kl+ku+1)]);
//      return Vector<Ref,AT>(n-abs(i),ku+kl+1,true,memory,ele+ku-i + (i>0?i:0)*(kl+ku+1));
//      //(ku-i,i>0?i:0));
//    }
//
//  template <class AT>
//    inline const Vector<Ref,AT> Matrix<GeneralBand,Ref,Ref,AT>::operator()(int i) const {
//      FMATVEC_ASSERT(i<=ku, AT);
//      FMATVEC_ASSERT(i>=-kl, AT);
//
//      return Vector<Ref,AT>(n-abs(i),ku+kl+1,true,memory,ele+ku-i + (i>0?i:0)*(kl+ku+1));
//      //(ku-i,i>0?i:0));
//    }
 
  template <class AT>
    inline Matrix<GeneralBand,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]=operator()(r,c);
      }
      return ret;
    }
 
 
  template <class AT>
    inline Matrix<GeneralBand,Ref,Ref,AT>& Matrix<GeneralBand,Ref,Ref,AT>::copy(const Matrix<GeneralBand,Ref,Ref,AT> &A) {
      for(int i=0; i<kl+ku+1; i++)
    for(int j=0; j<n; j++)
      ele[i+j*(kl+ku+1)] = A.ele[i+j*(kl+ku+1)];
      return *this;
    }
 
 
}
 
#endif