__host__ __device__ cuvec4::cuvec4() 
	{
		x = 0; y = 0; z = 0; w = 0; 
		return;
	}

	__host__ __device__ cuvec4::~cuvec4() 
	{
		x = 0; y = 0; z = 0; w = 0; 
		return;
	}

	__host__ __device__ cuvec4::cuvec4(const double &_x, const double &_y, const double &_z, const double &_w) 
	{
		x = _x; y = _y; z = _z; w = _w; 
		return;
	}

	__host__ __device__ double& cuvec4::operator[](const int &I)
	{
		switch(I)
		{
			case 0:
				return x;
			case 1:
				return y;
			case 2:
				return z;
			case 3:
				return w;
		}

		return x;
	}

	__host__ __device__ const double& cuvec4::operator[](const int &I) const
	{
		switch(I)
		{
			case 0:
				return x;
			case 1:
				return y;
			case 2:
				return z;
			case 3:
				return w;
		}

		return x;
	}

	__host__ __device__ cuvec4 cuvec4::operator+(const cuvec4& rhs) const
	{
		cuvec4 ret;
		ret.x = x + rhs.x;
		ret.y = y + rhs.y;
		ret.z = z + rhs.z;
		ret.w = w + rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 cuvec4::operator-(const cuvec4& rhs) const
	{
		cuvec4 ret;
		ret.x = x - rhs.x;
		ret.y = y - rhs.y;
		ret.z = z - rhs.z;
		ret.w = w - rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 cuvec4::operator*(const cuvec4& rhs) const
	{
		//Elementwise product
		cuvec4 ret;
		ret.x = x * rhs.x;
		ret.y = y * rhs.y;
		ret.z = z * rhs.z;
		ret.w = w * rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 cuvec4::operator/(const cuvec4& rhs) const
	{
		//Elementwise division
		cuvec4 ret;
		ret.x = x / rhs.x;
		ret.y = y / rhs.y;
		ret.z = z / rhs.z;
		ret.w = w / rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 operator*(const cuvec4& lhs, const double& rhs)
	{
		cuvec4 ret;
		ret.x = lhs.x*rhs;
		ret.y = lhs.y*rhs;
		ret.z = lhs.z*rhs;
		ret.w = lhs.w*rhs;
		return ret;
	}

	__host__ __device__ cuvec4 operator*(const double& lhs, const cuvec4& rhs)
	{
		cuvec4 ret;
		ret.x = lhs*rhs.x;
		ret.y = lhs*rhs.y;
		ret.z = lhs*rhs.z;
		ret.w = lhs*rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 operator/(const cuvec4& lhs, const double& rhs)
	{
		cuvec4 ret;
		ret.x = lhs.x/rhs;
		ret.y = lhs.y/rhs;
		ret.z = lhs.z/rhs;
		ret.w = lhs.w/rhs;
		return ret;
	}

	__host__ __device__ cuvec4 operator/(const double& lhs, const cuvec4& rhs)
	{
		cuvec4 ret;
		ret.x = lhs/rhs.x;
		ret.y = lhs/rhs.y;
		ret.z = lhs/rhs.z;
		ret.w = lhs/rhs.w;
		return ret;
	}

	__host__ __device__ cuvec4 operator-(const cuvec4& other)
	{
		cuvec4 ret;
		ret.x = -other.x;
		ret.y = -other.y;
		ret.z = -other.z;
		ret.w = -other.w;
		return ret;
	}

	__host__ __device__ cuvec4& cuvec4::operator+=(const cuvec4& rhs)
	{
		x += rhs.x;
		y += rhs.y;
		z += rhs.z;
		w += rhs.w;
		return *this;
	}

	__host__ __device__ cuvec4& cuvec4::operator-=(const cuvec4& rhs)
	{
		x -= rhs.x;
		y -= rhs.y;
		z -= rhs.z;
		w -= rhs.w;
		return *this;
	}

	__host__ __device__ cuvec4& cuvec4::operator*=(const double& rhs)
	{
		x *= rhs;
		y *= rhs;
		z *= rhs;
		w *= rhs;
		return *this;
	}

	__host__ __device__ cuvec4& cuvec4::operator/=(const double& rhs)
	{
		x /= rhs;
		y /= rhs;
		z /= rhs;
		w /= rhs;
		return *this;
	}


//
//Matrix Header Stuff
//

	__host__ __device__ cumat4::cumat4() 
	{
		m00 = 0;
		m01 = 0;
		m02 = 0;
		m03 = 0;

		m10 = 0;
		m11 = 0;
		m12 = 0;
		m13 = 0;

		m20 = 0;
		m21 = 0;
		m22 = 0;
		m23 = 0;

		m30 = 0;
		m31 = 0;
		m32 = 0;
		m33 = 0;

		return;
	}

	__host__ __device__ cumat4::~cumat4() 
	{
		//m00 = 0;
		//m01 = 0;
		//m02 = 0;
		//m03 = 0;

		//m10 = 0;
		//m11 = 0;
		//m12 = 0;
		//m13 = 0;

		//m20 = 0;
		//m21 = 0;
		//m22 = 0;
		//m23 = 0;

		//m30 = 0;
		//m31 = 0;
		//m32 = 0;
		//m33 = 0;

		return;
	}

	__host__ __device__ cumat4::cumat4(
		const double& _m00, const double& _m01, const double& _m02, const double& _m03, 
		const double& _m10, const double& _m11, const double& _m12, const double& _m13, 
		const double& _m20, const double& _m21, const double& _m22, const double& _m23, 
		const double& _m30, const double& _m31, const double& _m32, const double& _m33
	)
	{
		m00 = _m00; 
		m10 = _m10; 
		m20 = _m20; 
		m30 = _m30; 

		m01 = _m01; 
		m11 = _m11; 
		m21 = _m21; 
		m31 = _m31; 

		m02 = _m02; 
		m12 = _m12; 
		m22 = _m22; 
		m32 = _m32; 

		m03 = _m03; 
		m13 = _m13; 
		m23 = _m23; 
		m33 = _m33; 


		return;
	}

	__host__ __device__ cumat4::cumat4(const double* data16)
	{
		m00 = data16[0]; 
		m10 = data16[1]; 
		m20 = data16[2]; 
		m30 = data16[3]; 

		m01 = data16[4]; 
		m11 = data16[5]; 
		m21 = data16[6]; 
		m31 = data16[7]; 

		m02 = data16[8]; 
		m12 = data16[9]; 
		m22 = data16[10]; 
		m32 = data16[11]; 

		m03 = data16[12]; 
		m13 = data16[13]; 
		m23 = data16[14]; 
		m33 = data16[15]; 


		return;
	}

	__host__ __device__ double& cumat4::operator[](const int &I)
	{
		switch(I)
		{
			case 0:
				return m00;
			case 1:
				return m10;
			case 2:
				return m20;
			case 3:
				return m30;
			case 4:
				return m01;
			case 5:
				return m11;
			case 6:
				return m21;
			case 7:
				return m31;
			case 8:
				return m02;
			case 9:
				return m12;
			case 10:
				return m22;
			case 11:
				return m32;
			case 12:
				return m03;
			case 13:
				return m13;
			case 14:
				return m23;
			case 15:
				return m33;
		}

		return m00;
	}

	__host__ __device__ const double& cumat4::operator[](const int &I) const
	{
		switch(I)
		{
			case 0:
				return m00;
			case 1:
				return m10;
			case 2:
				return m20;
			case 3:
				return m30;
			case 4:
				return m01;
			case 5:
				return m11;
			case 6:
				return m21;
			case 7:
				return m31;
			case 8:
				return m02;
			case 9:
				return m12;
			case 10:
				return m22;
			case 11:
				return m32;
			case 12:
				return m03;
			case 13:
				return m13;
			case 14:
				return m23;
			case 15:
				return m33;
		}

		return m00;
	}

	__host__ __device__ double& cumat4::operator()(const int &I, const int &J)
	{
		return (*this)[I+4*J];
	}

	__host__ __device__ const double& cumat4::operator()(const int &I, const int &J) const
	{
		return (*this)[I+4*J];
	}

	__host__ __device__ double& cumat4::at(const int &I, const int &J)
	{
		return (*this)[I+4*J];
	}

	__host__ __device__ const double& cumat4::at(const int &I, const int &J) const
	{
		return (*this)[I+4*J];
	}

__host__ __device__ double* cumat4::data()
{
	return (double*)this;
}

__host__ __device__ const double* cumat4::data() const
{
	return (double*)this;
}

__host__ __device__ cumat4 cumat4::operator+(const cumat4& rhs) const
{
	cumat4 ret;
	ret.m00 = m00 + rhs.m00;
	ret.m10 = m10 + rhs.m10;
	ret.m20 = m20 + rhs.m20;
	ret.m30 = m30 + rhs.m30;

	ret.m01 = m01 + rhs.m01;
	ret.m11 = m11 + rhs.m11;
	ret.m21 = m21 + rhs.m21;
	ret.m31 = m31 + rhs.m31;

	ret.m02 = m02 + rhs.m02;
	ret.m12 = m12 + rhs.m12;
	ret.m22 = m22 + rhs.m22;
	ret.m32 = m32 + rhs.m32;

	ret.m03 = m03 + rhs.m03;
	ret.m13 = m13 + rhs.m13;
	ret.m23 = m23 + rhs.m23;
	ret.m33 = m33 + rhs.m33;

	return ret;
}

__host__ __device__ cumat4 cumat4::operator-(const cumat4& rhs) const
{
	cumat4 ret;
	ret.m00 = m00 - rhs.m00;
	ret.m10 = m10 - rhs.m10;
	ret.m20 = m20 - rhs.m20;
	ret.m30 = m30 - rhs.m30;

	ret.m01 = m01 - rhs.m01;
	ret.m11 = m11 - rhs.m11;
	ret.m21 = m21 - rhs.m21;
	ret.m31 = m31 - rhs.m31;

	ret.m02 = m02 - rhs.m02;
	ret.m12 = m12 - rhs.m12;
	ret.m22 = m22 - rhs.m22;
	ret.m32 = m32 - rhs.m32;

	ret.m03 = m03 - rhs.m03;
	ret.m13 = m13 - rhs.m13;
	ret.m23 = m23 - rhs.m23;
	ret.m33 = m33 - rhs.m33;

	return ret;
}

__host__ __device__ cumat4 cumat4::operator*(const cumat4& rhs) const
{
	cumat4 ret; //should be zeroed in constructor
	ret.m00 = m00*rhs.m00 + m01*rhs.m10 + m02*rhs.m20 + m03*rhs.m30;
	ret.m01 = m00*rhs.m01 + m01*rhs.m11 + m02*rhs.m21 + m03*rhs.m31;
	ret.m02 = m00*rhs.m02 + m01*rhs.m12 + m02*rhs.m22 + m03*rhs.m32;
	ret.m03 = m00*rhs.m03 + m01*rhs.m13 + m02*rhs.m23 + m03*rhs.m33;

	ret.m10 = m10*rhs.m00 + m11*rhs.m10 + m12*rhs.m20 + m13*rhs.m30;
	ret.m11 = m10*rhs.m01 + m11*rhs.m11 + m12*rhs.m21 + m13*rhs.m31;
	ret.m12 = m10*rhs.m02 + m11*rhs.m12 + m12*rhs.m22 + m13*rhs.m32;
	ret.m13 = m10*rhs.m03 + m11*rhs.m13 + m12*rhs.m23 + m13*rhs.m33;

	ret.m20 = m20*rhs.m00 + m21*rhs.m10 + m22*rhs.m20 + m23*rhs.m30;
	ret.m21 = m20*rhs.m01 + m21*rhs.m11 + m22*rhs.m21 + m23*rhs.m31;
	ret.m22 = m20*rhs.m02 + m21*rhs.m12 + m22*rhs.m22 + m23*rhs.m32;
	ret.m23 = m20*rhs.m03 + m21*rhs.m13 + m22*rhs.m23 + m23*rhs.m33;

	ret.m30 = m30*rhs.m00 + m31*rhs.m10 + m32*rhs.m20 + m33*rhs.m30;
	ret.m31 = m30*rhs.m01 + m31*rhs.m11 + m32*rhs.m21 + m33*rhs.m31;
	ret.m32 = m30*rhs.m02 + m31*rhs.m12 + m32*rhs.m22 + m33*rhs.m32;
	ret.m33 = m30*rhs.m03 + m31*rhs.m13 + m32*rhs.m23 + m33*rhs.m33;

	return ret;
}

__host__ __device__ cumat4 operator*(const cumat4& lhs, const double& rhs)
{
	cumat4 ret;
	ret.m00=lhs.m00*rhs;
	ret.m10=lhs.m10*rhs;
	ret.m20=lhs.m20*rhs;
	ret.m30=lhs.m30*rhs;
	ret.m01=lhs.m01*rhs;
	ret.m11=lhs.m11*rhs;
	ret.m21=lhs.m21*rhs;
	ret.m31=lhs.m31*rhs;
	ret.m02=lhs.m02*rhs;
	ret.m12=lhs.m12*rhs;
	ret.m22=lhs.m22*rhs;
	ret.m32=lhs.m32*rhs;
	ret.m03=lhs.m03*rhs;
	ret.m13=lhs.m13*rhs;
	ret.m23=lhs.m23*rhs;
	ret.m33=lhs.m33*rhs;
	return ret;
}

__host__ __device__ cumat4 operator/(const cumat4& lhs, const double& rhs)
{
	cumat4 ret;
	ret.m00=lhs.m00/rhs;
	ret.m10=lhs.m10/rhs;
	ret.m20=lhs.m20/rhs;
	ret.m30=lhs.m30/rhs;
	ret.m01=lhs.m01/rhs;
	ret.m11=lhs.m11/rhs;
	ret.m21=lhs.m21/rhs;
	ret.m31=lhs.m31/rhs;
	ret.m02=lhs.m02/rhs;
	ret.m12=lhs.m12/rhs;
	ret.m22=lhs.m22/rhs;
	ret.m32=lhs.m32/rhs;
	ret.m03=lhs.m03/rhs;
	ret.m13=lhs.m13/rhs;
	ret.m23=lhs.m23/rhs;
	ret.m33=lhs.m33/rhs;
	return ret;
}

__host__ __device__ cumat4 operator*(const double& lhs, const cumat4& rhs)
{
	cumat4 ret;
	ret.m00=lhs*rhs.m00;
	ret.m10=lhs*rhs.m10;
	ret.m20=lhs*rhs.m20;
	ret.m30=lhs*rhs.m30;
	ret.m01=lhs*rhs.m01;
	ret.m11=lhs*rhs.m11;
	ret.m21=lhs*rhs.m21;
	ret.m31=lhs*rhs.m31;
	ret.m02=lhs*rhs.m02;
	ret.m12=lhs*rhs.m12;
	ret.m22=lhs*rhs.m22;
	ret.m32=lhs*rhs.m32;
	ret.m03=lhs*rhs.m03;
	ret.m13=lhs*rhs.m13;
	ret.m23=lhs*rhs.m23;
	ret.m33=lhs*rhs.m33;
	return ret;
}

__host__ __device__ cuvec4 operator*(const cumat4& lhs, const cuvec4& rhs)
{
	cuvec4 ret;
	ret.x = lhs.m00*rhs.x + lhs.m01*rhs.y + lhs.m02*rhs.z + lhs.m03*rhs.w;
	ret.y = lhs.m10*rhs.x + lhs.m11*rhs.y + lhs.m12*rhs.z + lhs.m13*rhs.w;
	ret.z = lhs.m20*rhs.x + lhs.m21*rhs.y + lhs.m22*rhs.z + lhs.m23*rhs.w;
	ret.w = lhs.m30*rhs.x + lhs.m31*rhs.y + lhs.m32*rhs.z + lhs.m33*rhs.w;
	return ret;
}

__host__ __device__ cuvec4 operator*(const cuvec4& lhs, const cumat4& rhs)
{
	cuvec4 ret;
	ret.x = lhs.x*rhs.m00 + lhs.y*rhs.m10 + lhs.z*rhs.m20 + lhs.w*rhs.m30;
	ret.y = lhs.x*rhs.m01 + lhs.y*rhs.m11 + lhs.z*rhs.m21 + lhs.w*rhs.m31;
	ret.z = lhs.x*rhs.m02 + lhs.y*rhs.m12 + lhs.z*rhs.m22 + lhs.w*rhs.m32;
	ret.w = lhs.x*rhs.m03 + lhs.y*rhs.m13 + lhs.z*rhs.m23 + lhs.w*rhs.m33;
	return ret;
}

__host__ __device__ cumat4 operator-(const cumat4& rhs)
{
	cumat4 ret;
	ret.m00 = -rhs.m00;
	ret.m10 = -rhs.m10;
	ret.m20 = -rhs.m20;
	ret.m30 = -rhs.m30;
	ret.m01 = -rhs.m01;
	ret.m11 = -rhs.m11;
	ret.m21 = -rhs.m21;
	ret.m31 = -rhs.m31;
	ret.m02 = -rhs.m02;
	ret.m12 = -rhs.m12;
	ret.m22 = -rhs.m22;
	ret.m32 = -rhs.m32;
	ret.m03 = -rhs.m03;
	ret.m13 = -rhs.m13;
	ret.m23 = -rhs.m23;
	ret.m33 = -rhs.m33;
	return ret;
}

__host__ __device__ cumat4& cumat4::operator+=(const cumat4& rhs)
{
	m00 += rhs.m00;
	m10 += rhs.m10;
	m20 += rhs.m20;
	m30 += rhs.m30;
	m01 += rhs.m01;
	m11 += rhs.m11;
	m21 += rhs.m21;
	m31 += rhs.m31;
	m02 += rhs.m02;
	m12 += rhs.m12;
	m22 += rhs.m22;
	m32 += rhs.m32;
	m03 += rhs.m03;
	m13 += rhs.m13;
	m23 += rhs.m23;
	m33 += rhs.m33;
	return *this;
}

__host__ __device__ cumat4& cumat4::operator-=(const cumat4& rhs)
{
	m00 -= rhs.m00;
	m10 -= rhs.m10;
	m20 -= rhs.m20;
	m30 -= rhs.m30;
	m01 -= rhs.m01;
	m11 -= rhs.m11;
	m21 -= rhs.m21;
	m31 -= rhs.m31;
	m02 -= rhs.m02;
	m12 -= rhs.m12;
	m22 -= rhs.m22;
	m32 -= rhs.m32;
	m03 -= rhs.m03;
	m13 -= rhs.m13;
	m23 -= rhs.m23;
	m33 -= rhs.m33;
	return *this;
}

__host__ __device__ cumat4& cumat4::operator*=(const double& rhs)
{
	m00 *= rhs;
	m10 *= rhs;
	m20 *= rhs;
	m30 *= rhs;
	m01 *= rhs;
	m11 *= rhs;
	m21 *= rhs;
	m31 *= rhs;
	m02 *= rhs;
	m12 *= rhs;
	m22 *= rhs;
	m32 *= rhs;
	m03 *= rhs;
	m13 *= rhs;
	m23 *= rhs;
	m33 *= rhs;
	return *this;
}

__host__ __device__ cumat4& cumat4::operator/=(const double& rhs)
{
	m00 /= rhs;
	m10 /= rhs;
	m20 /= rhs;
	m30 /= rhs;
	m01 /= rhs;
	m11 /= rhs;
	m21 /= rhs;
	m31 /= rhs;
	m02 /= rhs;
	m12 /= rhs;
	m22 /= rhs;
	m32 /= rhs;
	m03 /= rhs;
	m13 /= rhs;
	m23 /= rhs;
	m33 /= rhs;
	return *this;
}

__host__ __device__ cumat4& cumat4::operator*=(const cumat4& rhs)
{
	cumat4 tmp = *this;
	m00 = tmp.m00*rhs.m00 + tmp.m01*rhs.m10 + tmp.m02*rhs.m20 + tmp.m03*rhs.m30;
	m01 = tmp.m00*rhs.m01 + tmp.m01*rhs.m11 + tmp.m02*rhs.m21 + tmp.m03*rhs.m31;
	m02 = tmp.m00*rhs.m02 + tmp.m01*rhs.m12 + tmp.m02*rhs.m22 + tmp.m03*rhs.m32;
	m03 = tmp.m00*rhs.m03 + tmp.m01*rhs.m13 + tmp.m02*rhs.m23 + tmp.m03*rhs.m33;
	m10 = tmp.m10*rhs.m00 + tmp.m11*rhs.m10 + tmp.m12*rhs.m20 + tmp.m13*rhs.m30;
	m11 = tmp.m10*rhs.m01 + tmp.m11*rhs.m11 + tmp.m12*rhs.m21 + tmp.m13*rhs.m31;
	m12 = tmp.m10*rhs.m02 + tmp.m11*rhs.m12 + tmp.m12*rhs.m22 + tmp.m13*rhs.m32;
	m13 = tmp.m10*rhs.m03 + tmp.m11*rhs.m13 + tmp.m12*rhs.m23 + tmp.m13*rhs.m33;
	m20 = tmp.m20*rhs.m00 + tmp.m21*rhs.m10 + tmp.m22*rhs.m20 + tmp.m23*rhs.m30;
	m21 = tmp.m20*rhs.m01 + tmp.m21*rhs.m11 + tmp.m22*rhs.m21 + tmp.m23*rhs.m31;
	m22 = tmp.m20*rhs.m02 + tmp.m21*rhs.m12 + tmp.m22*rhs.m22 + tmp.m23*rhs.m32;
	m23 = tmp.m20*rhs.m03 + tmp.m21*rhs.m13 + tmp.m22*rhs.m23 + tmp.m23*rhs.m33;
	m30 = tmp.m30*rhs.m00 + tmp.m31*rhs.m10 + tmp.m32*rhs.m20 + tmp.m33*rhs.m30;
	m31 = tmp.m30*rhs.m01 + tmp.m31*rhs.m11 + tmp.m32*rhs.m21 + tmp.m33*rhs.m31;
	m32 = tmp.m30*rhs.m02 + tmp.m31*rhs.m12 + tmp.m32*rhs.m22 + tmp.m33*rhs.m32;
	m33 = tmp.m30*rhs.m03 + tmp.m31*rhs.m13 + tmp.m32*rhs.m23 + tmp.m33*rhs.m33;
	return *this;
}

__host__ __device__ cumat4 cumat4::transpose() const
{
	cumat4 ret;
	ret.m00 = m00;
	ret.m10 = m01;
	ret.m20 = m02;
	ret.m30 = m03;
	ret.m01 = m10;
	ret.m11 = m11;
	ret.m21 = m12;
	ret.m31 = m13;
	ret.m02 = m20;
	ret.m12 = m21;
	ret.m22 = m22;
	ret.m32 = m23;
	ret.m03 = m30;
	ret.m13 = m31;
	ret.m23 = m32;
	ret.m33 = m33;
	return ret;
}