#ifndef __CUVECT3_HPP__
#define __CUVECT3_HPP__

namespace amscuda
{

    class cuvect3
    {
        public:
        double x;
        double y;
        double z;

        __host__ __device__ cuvect3();
        __host__ __device__ ~cuvect3();
        __host__ __device__ cuvect3(double _x, double _y, double _z);
        
        
        __host__ __device__ double& operator[](const int I);
        __host__ __device__ const double& operator[](const int I) const;

        __host__ __device__ cuvect3 operator+(cuvect3 lhs);
        __host__ __device__ cuvect3 operator-(cuvect3 lhs);
        __host__ __device__ cuvect3 operator*(double lhs);
        __host__ __device__ cuvect3 operator/(double lhs);
        __host__ __device__ friend cuvect3 operator-(cuvect3 rhs);
    };

    class cumat3
    {
        public:
        double dat[9];

        __host__ __device__ cumat3();
        __host__ __device__ ~cumat3();
        __host__ __device__ double& operator[](const int I);
        __host__ __device__ double& operator()(const int I, const int J);
        __host__ __device__ double& at(const int I, const int J);

        __host__ __device__ cumat3 operator+(cumat3 lhs);
        __host__ __device__ cumat3 operator-(cumat3 lhs);
        __host__ __device__ cumat3 operator*(double lhs);
        __host__ __device__ cumat3 operator/(double lhs);
        __host__ __device__ cuvect3 operator*(cuvect3 lhs);
        __host__ __device__ cumat3 operator*(cumat3 lhs);
        __host__ __device__ friend cumat3 operator-(cumat3 rhs);

        __host__ __device__ double det();
        __host__ __device__ cumat3 transpose();
        __host__ __device__ cumat3 inverse();
    };

    __host__ __device__ double cuvect3_dot(cuvect3 a, cuvect3 b);
    __host__ __device__ cuvect3 cuvect3_cross(cuvect3 a, cuvect3 b);
    __host__ __device__ double cuvect3_norm(cuvect3 a);
    __host__ __device__ cuvect3 cuvect3_normalize(cuvect3 a);
    __host__ __device__ cuvect3 cuvect3_proj(cuvect3 a, cuvect3 b);

    //3x3 matrix operations
    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
    
    //transposes a 3x3 (9 element) matrix
    __host__ __device__ void mat3_transpose(double *mat3inout);

    //copies src to dest
    __host__ __device__ void mat3_copy(double *mat3_dest, const double *mat3_src);

    //returns determinant of 3x3 matrix
    __host__ __device__ double mat3_det(double *mat3in);

    //inverts a 3x3 (9 element) matrix
    __host__ __device__ void mat3_inverse(double *mat3inout);

    __host__ __device__ cuvect3 mat3_mult(double *mat3in, cuvect3 cvin);
    __host__ __device__ void mat3_mult(double *matina, double *matinb, double *matout);

    __host__ __device__ void mat3_hodgedual(cuvect3 vecin, double *matout);
    __host__ __device__ void mat3_hodgedual(double *matin, cuvect3 vecout);

    //returns direction cosine rotation matrix from axis and angle
    __host__ __device__ void mat3_rot_from_axisangle(cuvect3 axis, double angle, double *matout);

    __host__ void test_cudavect_logic1();

}; //end namespace amscuda

#endif