/*! \file inline-matrix.h
 *
 * \author John Reppy & Lamont Samuels
 */

/*
 * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
 * All rights reserved.
 */

#ifndef _DIDEROT_CL_INLINE_MATRIX_H_
#define _DIDEROT_CL_INLINE_MATRIX_H_

#ifndef _DIDEROT_CL_TYPES_H_
#include "cl-types.h"
#endif

/********** 2x2 matrix functions **********/


inline void zero2x2 (Diderot_Mat2x2_t dst)
{
    dst[0] = (Diderot_vec2_t)(0.0, 0.0);
    dst[1] = (Diderot_vec2_t)(0.0, 0.0);
}

inline void identity2x2 (Diderot_Mat2x2_t dst)
{
    dst[0] = (Diderot_vec2_t)(1.0, 0.0);
    dst[1] = (Diderot_vec2_t)(0.0, 1.0);
}

inline void copy2x2 (Diderot_Mat2x2_t dst, Diderot_Mat2x2_t src)
{
    dst[0] = src[0];
    dst[1] = src[1];
}

inline void scale2x2 (Diderot_Mat2x2_t dst, float s, Diderot_Mat2x2_t src)
{
    Diderot_vec2_t scale = (Diderot_vec2_t)(s, s);
    dst[0] = scale * src[0];
    dst[1] = scale * src[1];
}

inline void add2x2 (Diderot_Mat2x2_t dst, Diderot_Mat2x2_t a, Diderot_Mat2x2_t b)
{
    dst[0] = a[0] + b[0];
    dst[1] = a[1] + b[1];
}

inline void sub2x2 (Diderot_Mat2x2_t dst, Diderot_Mat2x2_t a, Diderot_Mat2x2_t b)
{
    dst[0] = a[0] - b[0];
    dst[1] = a[1] - b[1];
}

inline Diderot_vec2_t mulVec3Mat2x2 (Diderot_vec2_t v, Diderot_Mat2x2_t m)
{
    return (Diderot_vec2_t)(
	dot(v, (Diderot_vec2_t)(m[0].s0, m[1].s0)),
	dot(v, (Diderot_vec2_t)(m[0].s1, m[1].s1)));
}

inline Diderot_vec2_t mulMat2x2Vec2 (Diderot_Mat2x2_t m, Diderot_vec2_t v)
{
    return (Diderot_vec2_t)(dot(m[0], v), dot(m[1], v));
}

inline void mulMat2x2Mat2x2 (Diderot_Mat2x2_t dst,  Diderot_Mat2x2_t m1,   Diderot_Mat2x2_t m2)
{
    dst[0] = (Diderot_vec2_t)(
	dot(m1[0], (Diderot_vec2_t)(m2[0].s0, m2[1].s0)),
	dot(m1[0], (Diderot_vec2_t)(m2[0].s1, m2[1].s1)));
    dst[1] = (Diderot_vec2_t)(
	dot(m1[1], (Diderot_vec2_t)(m2[0].s0, m2[1].s0)),
	dot(m1[1], (Diderot_vec2_t)(m2[0].s1, m2[1].s1)));
}

inline void transpose2x2 (Diderot_Mat2x2_t dst, Diderot_Mat2x2_t src)
{
    dst[0] = (Diderot_vec2_t)(src[0].s0, src[1].s0);
    dst[1] = (Diderot_vec2_t)(src[0].s1, src[1].s1);
}

// The Frobenius norm of a matrix is the sqrt of the sum of the squares of the elements
inline float norm2x2 (Diderot_Mat2x2_t m)
{
    return sqrt(dot(m[0],m[0]) + dot(m[1],m[1]));
}


/********** 3x3 matrix functions **********/

inline void zero3x3 (Diderot_Mat3x3_t dst)
{
    dst[0] = VEC3(0.0, 0.0, 0.0);
    dst[1] = VEC3(0.0, 0.0, 0.0);
    dst[2] = VEC3(0.0, 0.0, 0.0);
}

inline void identity3x3 (Diderot_Mat3x3_t dst)
{
    dst[0] = VEC3(1.0, 0.0, 0.0);
    dst[1] = VEC3(0.0, 1.0, 0.0);
    dst[2] = VEC3(0.0, 0.0, 1.0);
}

inline void copy3x3 (Diderot_Mat3x3_t dst, Diderot_Mat3x3_t src)
{
    dst[0] = src[0];
    dst[1] = src[1];
    dst[2] = src[2];
}

inline void scale3x3 (Diderot_Mat3x3_t dst, float s, Diderot_Mat3x3_t src)
{
    Diderot_vec3_t scale = VEC3(s, s, s);
    dst[0] = scale * src[0];
    dst[1] = scale * src[1];
    dst[2] = scale * src[2];
}

inline void add3x3 (Diderot_Mat3x3_t dst, Diderot_Mat3x3_t a, Diderot_Mat3x3_t b)
{
    dst[0] = a[0] + b[0];
    dst[1] = a[1] + b[1];
    dst[2] = a[2] + b[2];
}

inline void sub3x3 (Diderot_Mat3x3_t dst, Diderot_Mat3x3_t a, Diderot_Mat3x3_t b)
{
    dst[0] = a[0] - b[0];
    dst[1] = a[1] - b[1];
    dst[2] = a[2] - b[2];
}

inline Diderot_vec3_t mulVec3Mat3x3 (Diderot_vec3_t v, Diderot_Mat3x3_t m)
{
    return VEC3(
	dot(v, VEC3(m[0].s0, m[1].s0, m[2].s0)),
	dot(v, VEC3(m[0].s1, m[1].s1, m[2].s1)),
	dot(v, VEC3(m[0].s2, m[1].s2, m[2].s2)));
}

inline Diderot_vec3_t mulMat3x3Vec3 (Diderot_Mat3x3_t m, Diderot_vec3_t v)
{
    return VEC3(dot(m[0], v), dot(m[1], v), dot(m[2], v));
}

inline void mulMat3x3Mat3x3 (Diderot_Mat3x3_t dst, Diderot_Mat3x3_t m1, Diderot_Mat3x3_t m2)
{
    dst[0] = VEC3(
	dot(m1[0], VEC3(m2[0].s0, m2[1].s0, m2[2].s0)),
	dot(m1[0], VEC3(m2[0].s1, m2[1].s1, m2[2].s1)),
	dot(m1[0], VEC3(m2[0].s2, m2[1].s2, m2[2].s2)));
    dst[1] = VEC3(
	dot(m1[1], VEC3(m2[0].s0, m2[1].s0, m2[2].s0)),
	dot(m1[1], VEC3(m2[0].s1, m2[1].s1, m2[2].s1)),
	dot(m1[1], VEC3(m2[0].s2, m2[1].s2, m2[2].s2)));
    dst[2] = VEC3(
	dot(m1[2], VEC3(m2[0].s0, m2[1].s0, m2[2].s0)),
	dot(m1[2], VEC3(m2[0].s1, m2[1].s1, m2[2].s1)),
	dot(m1[2], VEC3(m2[0].s2, m2[1].s2, m2[2].s2)));
} 
inline void transpose3x3 (Diderot_Mat3x3_t dst, Diderot_Mat3x3_t src)
{
    dst[0] = VEC3(src[0].s0, src[1].s0, src[2].s0);
    dst[1] = VEC3(src[0].s1, src[1].s1, src[2].s1);
    dst[2] = VEC3(src[0].s2, src[1].s2, src[2].s2);
}

// The Frobenius norm of a matrix is the sqrt of the sum of the squares of the elements
inline float norm3x3 (Diderot_Mat3x3_t m)
{
    return sqrt(dot(m[0],m[0]) + dot(m[1],m[1]) + dot(m[2],m[2]));
}


/********** 4x4 matrix functions **********/

inline void copy4x4 (Diderot_Mat4x4_t dst, Diderot_Mat4x4_t src)
{
    dst[0] = src[0];
    dst[1] = src[1];
    dst[2] = src[2];
    dst[3] = src[3];
}

inline void transpose4x4 (Diderot_Mat4x4_t dst, Diderot_Mat4x4_t src)
{
    dst[0] = (Diderot_vec4_t)(src[0].s0, src[1].s0, src[2].s0, src[3].s0);
    dst[1] = (Diderot_vec4_t)(src[0].s1, src[1].s1, src[2].s1, src[3].s1);
    dst[2] = (Diderot_vec4_t)(src[0].s2, src[1].s2, src[2].s2, src[3].s2);
    dst[3] = (Diderot_vec4_t)(src[0].s3, src[1].s3, src[2].s3, src[3].s3);
}

// The Frobenius norm of a matrix is the sqrt of the sum of the squares of the elements
inline float norm4x4 (Diderot_Mat4x4_t m)
{
    return sqrt(dot(m[0],m[0]) + dot(m[1],m[1]) + dot(m[2],m[2]) + dot(m[3],m[3]));
}

#endif /* !_DIDEROT_CL_INLINE_MATRIX_H_ */