#ifndef RTS_MATRIX_H
#define RTS_MATRIX_H

//#include "rts/vector.h"
#include <string.h>
#include <iostream>
#include <stim/math/vector.h>
#include <stim/math/vec3.h>
#include <stim/cuda/cudatools/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 matrix(T rhs[N*N])
	{
		memcpy(M,rhs, sizeof(T)*N*N);
	}

	CUDA_CALLABLE matrix<T,N> set(T rhs[N*N])
	{
		memcpy(M, rhs, sizeof(T)*N*N);
		return *this;
	}

	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>
	vec<Y> operator*(vec<Y> rhs){
		unsigned int N = rhs.size();

		vec<Y> result;
		result.resize(N);

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

		return result;
	}

	template<typename Y>
	CUDA_CALLABLE vec3<Y> operator*(vec3<Y> rhs){
		vec3<Y> result = 0;
		for(int r=0; r<3; r++)
			for(int c=0; c<3; c++)
				result[r] += (*this)(r, c) * rhs[c];

		return result;
	}

	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