mthMatrix.h

Go to the documentation of this file.

/******************************************************************************

  Copyright 2015 Scientific Computation Research Center,
      Rensselaer Polytechnic Institute. All rights reserved.

  This work is open source software, licensed under the terms of the
  BSD license as described in the LICENSE file in the top-level directory.

*******************************************************************************/

#ifndef MTH_MATRIX_H
#define MTH_MATRIX_H

#include "mthVector.h"

namespace mth {

template <class T, unsigned M=0, unsigned N=0>
class Matrix : public can::Array<Vector<T,N>,M>
{
  public:
    Matrix() {}
    Matrix(unsigned m, unsigned n) { (void) m; (void) n; }
    unsigned rows() const {return M;}
    unsigned cols() const {return N;}
    T& operator()(unsigned i, unsigned j) {return (*this)[i][j];}
    T const& operator()(unsigned i, unsigned j) const {return (*this)[i][j];}
    Matrix<T,M,N> operator+(Matrix<T,M,N> const& b) const
    {
      Matrix<T,M,N> r;
      for (unsigned i=0; i < M; ++i)
        r.elems[i] = this->elems[i] + b.elems[i];
      return r;
    }
    Matrix<T,M,N> operator-(Matrix<T,M,N> const& b) const
    {
      Matrix<T,M,N> r;
      for (unsigned i=0; i < M; ++i)
        r.elems[i] = this->elems[i] - b.elems[i];
      return r;
    }
    Matrix<T,M,N> operator*(T const& s) const
    {
      Matrix<T,M,N> r;
      for (unsigned i=0; i < M; ++i)
        r.elems[i] = this->elems[i] * s;
      return r;
    }
    Matrix<T,M,N> operator/(T const& s) const
    {
      Matrix<T,M,N> r;
      for (unsigned i=0; i < M; ++i)
        r.elems[i] = this->elems[i] / s;
      return r;
    }
    Vector<T,M> operator*(Vector<T,N> const& b) const
    {
      Vector<T,M> r;
      for (unsigned i=0; i < M; ++i)
        r[i] = this->elems[i] * b;
      return r;
    }
    template <unsigned O>
    Matrix<T,M,O> operator*(Matrix<T,N,O> const& b) const
    {
      Matrix<T,M,O> r;
      for (unsigned i=0; i < M; ++i)
      for (unsigned j=0; j < O; ++j)
      {
        r[i][j] = this->elems[i][0]*b[0][j];
        for (unsigned k=1; k < N; ++k)
          r[i][j] += this->elems[i][k]*b[k][j];
      }
      return r;
    }
    void zero()
    {
      for (unsigned i=0; i < M; ++i)
        this->elems[i].zero();
    }
    void resize(unsigned m, unsigned n)
    {
      (void) m;
      (void) n;
    }
};

template <class T>
class Matrix<T,0,0>
{
  public:
    Matrix() : columns(0) {}
    Matrix(unsigned m, unsigned n) : columns(n), elems(m*n) {}
    unsigned rows() const {return elems.size()/columns;}
    unsigned cols() const {return columns;}
    void resize(unsigned m, unsigned n)
    {
      columns = n;
      elems.resize(m*n);
    }
    T& operator()(unsigned i, unsigned j)
    {
      return elems[i*columns + j];
    }
    T const& operator()(unsigned i, unsigned j) const
    {
      return elems[i*columns + j];
    }
    Matrix<T>& operator+=(Matrix<T> const& b)
    {
      for (unsigned i=0; i < this->elems.size(); ++i)
        this->elems[i] += b.elems[i];
      return *this;
    }
    Matrix<T>& operator-=(Matrix<T> const& b)
    {
      for (unsigned i=0; i < this->elems.size(); ++i)
        this->elems[i] -= b.elems[i];
      return *this;
    }
    Matrix<T>& operator*=(T const& s)
    {
      for (unsigned i=0; i < this->elems.size(); ++i)
        this->elems[i] *= s;
      return *this;
    }
    Matrix<T>& operator/=(T const& s)
    {
      for (unsigned i=0; i < this->elems.size(); ++i)
        this->elems[i] /= s;
      return *this;
    }
    void zero()
    {
      for (unsigned i=0; i < cols() * rows(); ++i)
        this->elems[i] = (T)0.0;
    }
  protected:
    unsigned columns;
    can::Array<T> elems;
};

template <class T>
class Matrix3x3 : public Matrix<T,3,3>
{
  public:
    Matrix3x3() {}
    Matrix3x3(
        T const& a11,T const& a12,T const& a13,
        T const& a21,T const& a22, T const& a23,
        T const& a31,T const& a32, T const& a33)
    {
      (*this)[0] = Vector3<T>(a11,a12,a13);
      (*this)[1] = Vector3<T>(a21,a22,a23);
      (*this)[2] = Vector3<T>(a31,a32,a33);
    }
    Matrix3x3(Matrix<T,3,3> const& other) : Matrix<T,3,3>(other) {}
    void toArray(T (*array)[3]) const
    {
      for (unsigned i=0; i < 3; ++i)
      for (unsigned j=0; j < 3; ++j)
        array[i][j] = (*this)[i][j];
    }
};

}

template <class T, unsigned M, unsigned N>
std::ostream& operator<<(std::ostream& s, mth::Matrix<T,M,N> const& a)
{
  for (unsigned i=0; i < a.rows(); ++i) {
    for (unsigned j=0; j < a.cols(); ++j)
      s << a(i,j) << ' ';
    s << '\n';
  }
  return s;
}

#endif