sparse_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 sparse_matrix_h
#define sparse_matrix_h
 
#include "matrix.h"
#include "square_matrix.h"
#include "types.h"
#include "_memory.h"
 
namespace fmatvec {
 
  template <class AT>
    class Matrix<Sparse,Ref,Ref,AT> {
 
      protected:
 
    
    Memory<AT> memEle;
    Memory<int> memI, memJ;
    AT *ele;
    int *I{nullptr}, *J{nullptr};
    int m{0}, n{0}, k{0};
 
        template <class Type, class Row, class Col> inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
    template <class Row> inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Symmetric,Row,Row,AT> &A);
    inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Sparse,Ref,Ref,AT> &A);
 
 
 
      public:
        static constexpr bool isVector {false};
 
        using value_type = AT;
        using shape_type = Sparse;
 
    explicit Matrix() : memEle(),  ele(0) { }
 
        explicit Matrix(int m_, int n_, int k_, Noinit) : memEle(k_), memI(m_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(m_), n(n_), k(k_) { }
        explicit Matrix(int m_, int n_, int k_, Init ini=INIT, const AT &a=AT()) : memEle(k_), memI(m_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(m_), n(n_), k(k_) { init(a); }
 
    Matrix(const Matrix<Sparse,Ref,Ref,AT> &A) : memEle(A.k), memI(A.n+1), memJ(A.k), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(A.m), n(A.n), k(A.k) {
          copy(A);
    }
 
        explicit Matrix(int m_, int n_, int k_, int *I_, int *J_, Init ini=INIT, const AT &a=AT()) : memEle(k_), memI(), memJ(), ele((AT*)memEle.get()), I(I_), J(J_), m(m_), n(n_), k(k_) {
      init(a);
        }
 
    ~Matrix() = default;
 
        Matrix<Sparse,Ref,Ref,AT>& resize(int m_, int n_, int k_, Noinit) {
            m = m_; n = n_; k = k_;
            memEle.resize(k);
            memI.resize(m+1);
            memJ.resize(k);
            ele = (AT*)memEle.get();
            I = (int*)memI.get();
            J = (int*)memJ.get();
            return *this;
        }
 
        Matrix<Sparse,Ref,Ref,AT>& resize(int m, int n, int k, Init ini=INIT, const AT &a=AT()) { return resize(m,n,k,Noinit()).init(a); }
 
        void resize(int m, int n) {
          resize(m,n,m*n);
        }
 
    inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Sparse,Ref,Ref,AT> &A) {
          FMATVEC_ASSERT(m == A.m, AT);
          FMATVEC_ASSERT(n == A.n, AT);
          FMATVEC_ASSERT(k == A.k, AT);
          return copy(A);
        }
 
        template<class Type, class Row, class Col>
    inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
          FMATVEC_ASSERT(m == A.rows(), AT);
          FMATVEC_ASSERT(n == A.cols(), AT);
          FMATVEC_ASSERT(k == A.nonZeroElements(), AT);
          return copy(A);
        }
 
        inline Matrix<Sparse,Ref,Ref,AT>& operator&=(Matrix<Sparse,Ref,Ref,AT> &A) {
          m=A.m;
          n=A.n;
          k=A.k;
          memEle = A.memEle;
          memI = A.memI;
          memJ = A.memJ;
          ele = A.ele;
          I = A.I;
          J = A.J;
          k = A.k;
          return *this;
        }
 
        template<class Type, class Row, class Col>
    inline Matrix<Sparse,Ref,Ref,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
          int k_ = A.nonZeroElements();
          if(m!=A.rows() || k!=k_) resize(A.rows(),A.cols(),k_,NONINIT);
          return copy(A);
        }
 
    const AT* operator()() const {
      return ele;
    }
 
    AT* operator()() {
      return ele;
    }
 
        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[pos(i,j)];
    }
 
    const AT& e(int i, int j) const {
      return ele[pos(i,j)];
    }
 
    int pos(int i, int j) const {
      for(int k=I[i]; k<I[i+1]; k++)
        if(J[k]==j)
          return k;
      return 0;
    }
 
        int ldim() {
          throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::ldim() cannot be called.");
        }
 
    const int* Ip() const {
      return I;
    }
 
    int* Ip() {
      return I;
    }
 
    const int* Jp() const {
      return J;
    }
 
    int* Jp() {
      return J;
    }
 
    int rows() const {return m;}
 
    int cols() const {return n;}
 
        CBLAS_ORDER blasOrder() const {
          throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::blasOrder() cannot be called.");
        }
 
    Matrix<Sparse,Ref,Ref,AT>& init(const AT &val);
        inline Matrix<Sparse,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
        inline Matrix<Sparse,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
 
        int nonZeroElements() const { return k; }
    };
  // ------------------------- Constructors -------------------------------------
  // ----------------------------------------------------------------------------
 
  template <class AT> template <class Type, class Row, class Col>
    inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
      int k=0;
      int i;
      for(i=0; i<A.rows(); i++) {
    ele[k]=A.e(i,i);
    J[k]=i;
    I[i]=k++;
    for(int j=0; j<i; j++) {
      // TODO eps
      if(fabs(A.e(i,j))>1e-16) {
        ele[k]=A.e(i,j);
        J[k++]=j;
      }
    }
    for(int j=i+1; j<A.cols(); j++) {
      // TODO eps
      if(fabs(A.e(i,j))>1e-16) {
        ele[k]=A.e(i,j);
        J[k++]=j;
      }
    }
      }
      if(n) I[i]=k;
      return *this;
    }
 
  template <class AT> template <class Row> 
    inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Symmetric,Row,Row,AT> &A) {
      int k=0;
      int i;
      for(i=0; i<A.size(); i++) {
        ele[k]=A.ej(i,i);
        J[k]=i;
        I[i]=k++;
        for(int j=0; j<i; j++) {
          // TODO eps
          if(fabs(A.ei(i,j))>1e-16) {
            ele[k]=A.ei(i,j);
            J[k++]=j;
          }
        }
        for(int j=i+1; j<A.size(); j++) {
          // TODO eps
          if(fabs(A.ej(i,j))>1e-16) {
            ele[k]=A.ej(i,j);
            J[k++]=j;
          }
        }
      }
      if(n) I[i]=k;
      return *this;
  }
 
  template <class AT>
    inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Sparse,Ref,Ref,AT> &A) {
      for(int i=0; i<=m; i++) {
    I[i] = A.I[i];
      }
      for(int i=0; i<k; i++) {
    ele[i] = A.ele[i];
    J[i] = A.J[i];
      }
      return *this;
    }
 
  template <class AT>
    Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::init(const AT& val) {
      for(int i=0; i<k; i++) {
    ele[i] = val;
      }
      return *this;
    }
 
}
 
#endif