math/matrix.h

#ifndef RTS_MATRIX_H
#define RTS_MATRIX_H

//#include "rts/vector.h"
#include <string.h>
#include <iostream>
#include "vector.h"
#include "../cuda/callable.h"

namespace stim{

template <class T, int N>
struct matrix
{
	//the matrix will be stored in column-major order (compatible with OpenGL)
	T M[N*N];

	CUDA_CALLABLE matrix()
	{
		for(int r=0; r<N; r++)
			for(int c=0; c<N; c++)
				if(r == c)
					(*this)(r, c) = 1;
				else
					(*this)(r, c) = 0;
	}

	CUDA_CALLABLE T& operator()(int row, int col)
	{
		return M[col * N + row];
	}

	CUDA_CALLABLE matrix<T, N> operator=(T rhs)
	{
		int Nsq = N*N;
		for(int i=0; i<Nsq; i++)
			M[i] = rhs;

		return *this;
	}

	template<typename Y>
	CUDA_CALLABLE vec<Y, N> operator*(vec<Y, N> rhs)
	{
		vec<Y, N> result;

		for(int r=0; r<N; r++)
			for(int c=0; c<N; c++)
				result[r] += (*this)(r, c) * rhs[c];

		return result;
	}

	CUDA_CALLABLE std::string toStr()
	{
		std::stringstream ss;

		for(int r = 0; r < N; r++)
		{
			ss<<"| ";
			for(int c=0; c<N; c++)
			{
				ss<<(*this)(r, c)<<" ";
			}
			ss<<"|"<<std::endl;
		}

		return ss.str();
	}


};

}	//end namespace rts

template <typename T, int N>
std::ostream& operator<<(std::ostream& os, stim::matrix<T, N> M)
{
    os<<M.toStr();
    return os;
}

//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
//template<class T, int N> using rtsMatrix = rts::matrix<T, N>;
//#endif

#endif