PhxMatrix.h

Go to the documentation of this file.

#ifndef PHX_MATRIX_H
#define PHX_MATRIX_H

#include <Phx/PhxConfig.h>
#include <Phx/PhxTypes.h>
#include <math.h>
#include <iostream>

namespace Phx {
  
  template <uint16_t Rows, uint16_t Cols, typename value_type = double>
  class Matrix {
    private:
    value_type mData[Rows*Cols]; // data in row major order
    
    public:
    Matrix() { // constructs zero matrix
      for (uint32_t i = 0; i < Rows*Cols; ++i) {
        mData[i] = 0;
      }
    }

    Matrix(const value_type* data) { // constructs from array
      if (!data) {
        std::cerr << "Constructing matrix from 0 array." << std::endl;
        throw RangeException();
      }

      for (uint32_t i = 0; i < Rows*Cols; ++i) {
        mData[i] = data[i];
      }
    }

    uint32_t rows() const { return Rows; }
    uint32_t columns() const { return Cols; }

    value_type element(uint16_t row, uint16_t col) const { 
      if (row >= rows() || col >= columns()) {
        std::cerr << "Illegal read access: " << row << ", " << col
                  << " in " << Rows << "x" << Cols << " matrix." << std::endl;
        return mData[0];
      }
      return mData[row*Cols+col];
    }

    value_type& element(uint16_t row, uint16_t col) {
      if (row >= rows() || col >= columns()) {
        std::cerr << "Illegal read access: " << row << ", " << col
                  << " in " << Rows << "x" << Cols << " matrix." << std::endl;
        return mData[0];
      }
      return mData[row*Cols+col];
    }

    void element(uint16_t row, uint16_t col, value_type value) {
      if (row >= rows() || col >= columns()) {
        std::cerr << "Illegal write access: " << row << ", " << col
                  << " in " << Rows << "x" << Cols << " matrix." << std::endl;
        return ;
      }
      mData[row*Cols+col] = value;
    }

    template <uint16_t RowRangeSize, uint16_t ColRangeSize>
    Matrix<RowRangeSize, ColRangeSize, value_type> 
    subMatrix(const Matrix<1, RowRangeSize, uint16_t>& rowRange, 
              const Matrix<1, ColRangeSize, uint16_t>& colRange) {
      Matrix<RowRangeSize, ColRangeSize, value_type> result;
      for (uint32_t i = 0; i < RowRangeSize; i++) {
        for (uint32_t j = 0; j < ColRangeSize; j++) {
          result.element(i,j, element(rowRange.element(0, i), colRange.element(0,j)));
        }
      }
      return result;
    }

    // Double/single indexing operators.
    /*value_type& operator()(uint16_t i, uint16_t j) {
      return element(i,j);
    }
    value_type& operator()(uint16_t i) {
      return element(i / Rows, i % Rows);
    }
    value_type operator()(uint16_t i, uint16_t j) const {
      return element(i,j);
    }
    value_type operator()(uint16_t i) const {
      return element(i / Rows, i % Rows);
      }*/

    template <uint32_t RowRangeSize, uint32_t ColRangeSize>
    Matrix<RowRangeSize, ColRangeSize, value_type> 
    operator()(const Matrix<1, RowRangeSize, uint16_t>& rowRange, 
               const Matrix<1, ColRangeSize, uint16_t>& colRange) {
      return subMatrix(rowRange, colRange);
    }

    template <uint16_t RhsCols>
    Matrix<Rows, RhsCols, value_type> operator*(const Matrix<Cols, RhsCols, value_type>& rhs) const {
      Matrix<Rows, RhsCols, value_type> result;
      for (uint32_t i = 0; i < Rows; i++) {
        for (uint32_t j = 0; j < RhsCols; j++) {
          for (uint32_t k = 0; k < Cols; k++) {
            result.element(i,j) += element(i,k)*rhs.element(k,j);
          }
        }
      }
      return result;
    }

    Matrix<Rows, Cols, value_type> operator+(const Matrix<Rows, Cols, value_type>& rhs) const {
      Matrix<Rows, Cols, value_type> result;
      for (uint32_t i = 0;  i < Rows*Cols;  i++) {
        result.mData[i] = mData[i] + rhs.data()[i];
      }
      return result;
    }

    Matrix<Rows, Cols, value_type> operator-(const Matrix<Rows, Cols, value_type>& rhs) const {
      Matrix<Rows, Cols, value_type> result;
      for (uint32_t i = 0;  i < Rows*Cols;  i++) {
        result.mData[i] = mData[i] - rhs.data()[i];
      }
      return result;
    }

    Matrix<Rows, Cols, value_type> operator*(value_type s) const {
      Matrix<Rows, Cols, value_type> result;
      for (uint32_t i = 0;  i < Rows*Cols;  i++) {
        result.mData[i] = mData[i] * s;
      }
      return result;      
    }

    Matrix<Rows, Cols, value_type> operator/(value_type s) const {
      Matrix<Rows, Cols, value_type> result;
      for (uint32_t i = 0;  i < Rows*Cols;  i++) {
        result.mData[i] = mData[i] / s;
      }
      return result;      
    }

    Matrix<Rows, Cols, value_type> operator-(void) const {
      Matrix<Rows, Cols, value_type> result;
      for (uint32_t i = 0;  i < Rows*Cols;  i++) {
        result.mData[i] = -mData[i];
      }
      return result;
    }

    Matrix<Cols, Rows, value_type> transpose(void) const {
      Matrix<Cols, Rows, value_type> result;
      for (uint16_t i = 0;  i < Rows;  i++) {
        for (uint16_t j = 0;  j < Cols;  j++) {
          result.element(j,i, element(i,j));
        }
      }
      return result;
    }

    static Matrix<Rows, Cols, value_type> identity() {
      Matrix<Rows, Cols, value_type> id;
      for (uint16_t i = 0; i < Rows && i < Cols; i++) {
        id.element(i,i) = 1;
      }
      return id;
    }

    static Matrix<Rows, Cols, value_type> zero() {
      return Matrix<Rows, Cols, value_type>();
    }

    value_type* data() { return mData; }
    const value_type* data() const { return mData; }

  private:

    void print(std::ostream& os) const {
      for (uint16_t i = 0; i < Rows; i++) {
        for (uint16_t j = 0; j < Cols; j++) {
          os.width(11);
          os.precision(11);
          os << element(i,j) << " ";
        }
        os << std::endl;
      }
    }

    void input(std::istream& is) {
      for (uint16_t i = 0; i < Rows; i++) {
        for (uint16_t j = 0; j < Cols; j++) {
          if (is.good()) {
            is >> element(i,j);
          } else {
            std::cerr << "Unable to read matrix from input stream." << std::endl;
            throw IOException();
          }
        }
      }
    }

    template <uint16_t Rows2, uint16_t Cols2, typename value_type2>
    friend std::ostream& operator<<(std::ostream& os, const Matrix<Rows2, Cols2, value_type2>& rhs);

    template <uint16_t Rows2, uint16_t Cols2, typename value_type2>
    friend std::istream& operator>>(std::istream& is, Matrix<Rows2, Cols2, value_type2>& rhs);
  };


  typedef Matrix<2,2>  Matrix2;
  typedef Matrix<3,3>  Matrix3;
  typedef Matrix<4,4>  Matrix4;
  typedef Matrix<2,1>  Vector2;
  typedef Matrix<3,1>  Vector3;
  typedef Matrix<4,1>  Vector4;

  // input / output operators
  template <uint16_t Rows, uint16_t Cols, typename value_type>
    std::ostream& operator<<(std::ostream& os, const Matrix<Rows, Cols, value_type>& rhs) {
    rhs.print(os);
    return os;
  }

  template <uint16_t Rows, uint16_t Cols, typename value_type>
    std::istream& operator>>(std::istream& is, Matrix<Rows, Cols, value_type>& rhs) {
    rhs.input(is);
    return is;
  }

}; // namespace Phx


#endif /* PHX_MATRIX_H */

---