---

Matrix.h

Go to the documentation of this file.

#ifndef FASTSHOWER_MATRIX_H
#define FASTSHOWER_MATRIX_H
#include <valarray>
#include "AtlfastUtils/Slice_iter.h"
#include "AtlfastUtils/Cslice_iter.h"
#include <iostream>
#include <cstddef>
#include <assert.h>
//=================================================
//
// Matrix
//
// Matrix class from Sjoustrup's "the C++ programming language" ch 22.4.6
//
//


namespace FastShower {
  using std::valarray;
  using std::slice;

  template<class T> 
    class Matrix{
  private:
    valarray<T>* m_v;
    size_t m_d1, m_d2;
  public:
    Matrix(size_t x, size_t y);
    Matrix(size_t x, size_t y, T init);
    Matrix(const Matrix&);
    Matrix& operator=(const Matrix&);
    ~Matrix(){delete m_v;}
    
    size_t size() const {return m_d1*m_d2;}
    size_t dim1() const {return m_d1;}
    size_t dim2() const {return m_d2;}
    
    Slice_iter<T> row(size_t i);
    Cslice_iter<T> row(size_t i)const;

    Slice_iter<T> column(size_t i);
    Cslice_iter<T> column(size_t i)const;

    T& operator()(size_t x, size_t y);
    T operator()(size_t x, size_t y)const;

    Slice_iter<T> operator()(size_t i) {return column(i);}
    Cslice_iter<T> operator()(size_t i) const {return column(i);}

    Slice_iter<T> operator[](size_t i) {return column(i);}
    Cslice_iter<T> operator[](size_t i) const {return column(i);}

    Matrix& operator*=(T);

    valarray<T>& array() {return *m_v;}

    Matrix sub(int minRow,int maxRow,int minCol,int maxCol) const;
    void sub(int minRow, int minCol, Matrix& dest) const;
  };




  template<class T> Matrix<T>::Matrix(size_t x, size_t y){
    m_d1=x;
    m_d2=y;
    m_v=new valarray<T>(x*y);
  }

  template<class T> Matrix<T>::Matrix(size_t x, size_t y, T init){
    m_d1=x;
    m_d2=y;
    m_v=new valarray<T>(x*y);
    for(size_t i=0; i<m_d1*m_d2; ++i)  m_v[i]=init;
  }

  template<class T> T& Matrix<T>::operator()(size_t x, size_t y){
    return column(x)[y];
  }

  template<class T> T Matrix<T>::operator() (size_t x, size_t y) const {
    return column(x)[y];
  }

  template<class T> 
  T mul(const Cslice_iter<T>& v1, 
             const valarray<T>& v2){
    T res=0;
    for( size_t i=0; i<v2.size(); i++) res+=v1[i]*v2[i];
    return res;
  }

  template<class T> 
  valarray<T> operator*(const Matrix<T>&m, const valarray<T>&v){
    valarray<T> res(m.dim2());
    for(size_t i=0; i<m.dim2(); i++) res[i]=mul(m.row(i), v);
    return res;
  }

  template<class T> 
  Matrix<T>& Matrix<T>::operator*=(T d){
    (*m_v) *= d;
    return *this;
  }

  template<class T> 
  Matrix<T> Matrix<T>::sub(int minRow, int maxRow, 
                     int minCol, int maxCol) const{
    int endRow=maxRow+1;
    int endCol=maxCol+1;

    assert(minRow>=0);
    assert(minCol>=0);
    assert(endRow<=m_d1);
    assert(endCol<=m_d2);

    int nRow=endRow-minRow;
    assert(nRow>0);
    
    int nCol=endCol-minCol;
    assert(nCol>0);

    Matrix<T> m(nRow, nCol);
    int ii=0;
    int jj;
    for(int i = minRow; i<endRow; ++i){
      jj=0;
      for(int j = minCol; j<endCol; ++j) m[ii][jj++]=(*this)[i][j];
      ++ii;
    }
    assert(ii == endRow);
    assert(jj == endCol);
    return m;
  }
  template<class T>
  void Matrix<T>::sub(int minRow, int minCol, Matrix& m) const{
    int endRow=minRow+m.m_d1;
    int endCol=minCol+m.m_d2;

    assert(minRow>=0);
    assert(minCol>=0);
    assert(endRow<=m_d1);
    assert(endCol<=m_d2);

    int nRow=endRow-minRow;
    assert(nRow>0);
    
    int nCol=endCol-minCol;
    assert(nCol>0);

    int ii=0;
    int jj;
    for(int i = minRow; i<endRow; ++i){
      jj=0;
      for(int j = minCol; j<endCol; ++j) m[ii][jj++]=(*this)[i][j];
      ++ii;
    }
    assert(ii == endRow);
    assert(jj == endCol);
    return;
  }

  template<class T> inline Slice_iter<T> Matrix<T>::row(size_t i){
    return Slice_iter<T>(m_v, slice(i,m_d1,m_d2));
  }

  template<class T> Cslice_iter<T> Matrix<T>::row(size_t i) const {
    return Cslice_iter<T>(m_v, slice(i,m_d1,m_d2));
  }

  template<class T> inline Slice_iter<T> Matrix<T>::column(size_t i){
    return Slice_iter<T>(m_v, slice(i*m_d2,m_d2,1));
  }

  template<class T> inline Cslice_iter<T> Matrix<T>::column(size_t i) const {
    return Cslice_iter<T>(m_v, slice(i*m_d2,m_d2,1));
  }
  template<class T>
  ostream& operator <<(ostream& o,const Atlfast::Matrix<T>& m){
    size_t end1 = m.dim1();
    size_t end2 = m.dim2();
    size_t i = 0;
    size_t j = 0;
    for(; i!= end1; ++i){
      for(j=0; j!= end2; ++j){
        cout<<m(i,j)<<" ";
      }
      cout<<endl;
    }
    return o;
  }
      
}//namespace

#endif

---