■

勉強がてらに、C++ で行列演算を表現するクラスと演算子と関数を約160行で作ってみた。遅延評価とかメモリ節約とか一切考えていない、最低限の機能のみの超シンプルなクラス。

  /*! \class Matrix
	  \brief a matrix with fixed size, which is (ROWS x COLS).
	  	  
	  \attention You may not be able to use Matrix class with a large size,
          because Matrix class uses stack memory.
  */
  template <size_t ROWS, size_t COLS, typename T>
  class Matrix {
  private:
	T elem[ROWS][COLS];
	
  public:
	Matrix() throw(){}
	template <typename S>
	explicit Matrix(const Matrix<ROWS,COLS,S> &a) throw(){  (*this) = a;}
	explicit Matrix(const T a[ROWS][COLS]) throw()
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		elem[i][j] = a[i][j];
	}

	//-----------------------------------------------------------------------
	//@{
	/** access operators to elem */
	const T& operator()(size_t i, size_t j) const{  return elem[i][j];}
	T& operator()(size_t i, size_t j){  return elem[i][j];}
	//@}

	//-----------------------------------------------------------------------
	bool operator==(const Matrix<ROWS,COLS,T> &a) const
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		if (elem[i][j] != a.elem[i][j])
		  return false;
	  return true;
	}
	bool operator!=(const Matrix<ROWS,COLS,T> &a) const{
	  return !(*this==a);}
	
	//-----------------------------------------------------------------------
	//@{
	/** assignment operators */
	template <typename S>
	Matrix<ROWS,COLS,T>& operator=(const Matrix<ROWS,COLS,S> &a)
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		elem[i][j] = a(i,j);
	  return *this;
	}

	Matrix<ROWS,COLS,T>& operator+=(const Matrix<ROWS,COLS,T> &a)
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		elem[i][j] += a(i,j);
	  return *this;
	}

	Matrix<ROWS,COLS,T>& operator-=(const Matrix<ROWS,COLS,T> &a)
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		elem[i][j] -= a(i,j);
	  return *this;
	}

	Matrix<ROWS,COLS,T>& operator*=(T a)
	{
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS; ++j)
		elem[i][j] *= a;
	  return *this;
	}

	Matrix<ROWS,COLS,T>& operator*=(const Matrix<ROWS,COLS,T> &a){
	  STATIC_ASSERT(ROWS==COLS);
	  *this = (*this) * a;
	  return *this;
	}
	//@}

	//-----------------------------------------------------------------------
	//@{
	/** binary operators */
	const Matrix<ROWS,COLS,T> operator+(const Matrix<ROWS,COLS,T> &a) const
	{
	  Matrix<ROWS,COLS,T> result(*this);
	  result += a;
	  return result;
	}

	const Matrix<ROWS,COLS,T> operator-(const Matrix<ROWS,COLS,T> &a) const
	{
	  Matrix<ROWS,COLS,T> result(*this);
	  result -= a;
	  return result;
	}

	const Matrix<ROWS,COLS,T> operator*(T a) const
	{
	  Matrix<ROWS,COLS,T> result(*this);
	  result *= a;
	  return result;
	}
	
	template <size_t COLS2>
	const Matrix<ROWS,COLS2,T> operator*(const Matrix<COLS,COLS2,T> &a) const
	{
	  Matrix<ROWS,COLS2,T> result;
	  for (size_t i=0; i<ROWS; ++i)
	  for (size_t j=0; j<COLS2; ++j){
		result(i,j) = 0;
		for (size_t k=0; k<COLS; ++k)
		  result(i,j) += elem[i][k] * a(k,j);
	  }
	  return result;
	}
	//@}
  };

  //-------------------------------------------------------------------
  //@{
  /** associate operators for Matrix */
  //-----------------------------------------------------------------------
  template <size_t ROWS, size_t COLS, typename T> inline
  const Matrix<ROWS,COLS,T> operator*(T a, const Matrix<ROWS,COLS,T> &m){
 	return m * a;}

  template <size_t ROWS, size_t COLS, typename T>
  std::ostream& operator<<(std::ostream& os, const Matrix<ROWS,COLS,T> &a){
	for (size_t i=0; i<ROWS; ++i){
	  for (size_t j=0; j<COLS; ++j)
		os << a(i,j) << " ";
	  os << endl;
	}
	return os;
  }
  //@}


  //-----------------------------------------------------------------------
  //@{
  /** associate functions for Matrix */
  //-----------------------------------------------------------------------
  template <size_t SIZE, typename T>
  const T trace(const Matrix<SIZE,SIZE,T> &a)
  {
	T result=0;
	for (size_t i=0; i<SIZE; ++i)
	  result += a(i,i);
	return result;
  }

  //-----------------------------------------------------------------------
  template <size_t ROWS, size_t COLS, typename T>
  const Matrix<COLS,ROWS,T> transpose(const Matrix<ROWS,COLS,T> &a)
  {
	Matrix<COLS,ROWS,T> result;
	for (size_t i=0; i<COLS; ++i)
	for (size_t j=0; j<ROWS; ++j)
	  result(i,j) = a(j,i);
	return result;
  }