Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /branches/lamont/src/include/Diderot/inline-tensor3.h
ViewVC logotype

View of /branches/lamont/src/include/Diderot/inline-tensor3.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2081 - (download) (as text) (annotate)
Mon Nov 5 23:26:06 2012 UTC (6 years, 11 months ago) by lamonts
File size: 9832 byte(s)
Creating new developmented branch based on vis12
/*! \file inline-tensor3.h
 *
 * \author John Reppy
 *
 * Operations on third-order tensors.
 */

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

#ifndef _DIDEROT_INLINE_TENSOR3_H_
#define _DIDEROT_INLINE_TENSOR3_H_

#ifndef _DIDEROT_TYPES_H_
#include "types.h"
#endif
#ifndef _DIDEROT_INLINE_VEC2_H_
#  include "inline-vec2.h"
#endif
#ifndef _DIDEROT_INLINE_VEC3_H_
#  include "inline-vec3.h"
#endif
#ifndef _DIDEROT_INLINE_VEC4_H_
#  include "inline-vec4.h"
#endif

/********** 2x2x2 tensor functions **********/

STATIC_INLINE void copyTen2x2x2 (Diderot_Ten2x2x2_t dst, Diderot_Ten2x2x2_t src)
{
    dst[0][0].v = src[0][0].v;
    dst[0][1].v = src[0][1].v;
    dst[1][0].v = src[1][0].v;
    dst[1][1].v = src[1][1].v;
}

// vector times tensor contraction: [v dot T]_jk = v_i * T_ijk
STATIC_INLINE void mulVec2Ten2x2x2 (Diderot_Mat2x2_t dst, Diderot_vec2_t v, Diderot_Ten2x2x2_t T)
{
    Diderot_union2_t u;
    u.v = v;
    dst[0].r[0] = u.r[0] * T[0][0].r[0] + u.r[1] * T[1][0].r[0];
    dst[0].r[1] = u.r[0] * T[0][0].r[1] + u.r[1] * T[1][0].r[1];
    dst[1].r[0] = u.r[0] * T[0][1].r[0] + u.r[1] * T[1][1].r[0];
    dst[1].r[1] = u.r[0] * T[0][1].r[1] + u.r[1] * T[1][1].r[1];
}

// tensor times vector contraction: [T dot v]_ij = T_ijk * v_k
STATIC_INLINE void mulTen2x2x2Vec2 (Diderot_Mat2x2_t dst, Diderot_Ten2x2x2_t T, Diderot_vec2_t v)
{
    dst[0].v = vec2(dot2(T[0][0].v, v), dot2(T[0][1].v, v));
    dst[1].v = vec2(dot2(T[1][0].v, v), dot2(T[1][1].v, v));
}

// The Frobenius norm of a tensor
STATIC_INLINE Diderot_real_t normTen2x2x2 (Diderot_Ten2x2x2_t T)
{
    Diderot_real_t sumSq =
	dot2(T[0][0].v,T[0][0].v) + dot2(T[0][1].v,T[0][1].v) +
	dot2(T[1][0].v,T[1][0].v) + dot2(T[1][1].v,T[1][1].v);
    return SQRT(sumSq);
}

// tensor times matrix ":" product
// (T : M)_i = T_ijk * M_jk
STATIC_INLINE Diderot_vec2_t colonProdTen2x2x2Mat2x2 (Diderot_Ten2x2x2_t T, Diderot_Mat2x2_t m)
{
#define V(i)	dot2(T[i][0].v, m[0].v) + dot2(T[i][1].v, m[1].v)

    return vec2(V(0), V(1));

#undef V
}

// matrix times tensor ":" product
// (M : T)_k = M_ij * T_ijk
STATIC_INLINE Diderot_vec2_t colonProdMat2x2Ten2x2x2 (Diderot_Mat2x2_t m, Diderot_Ten2x2x2_t T)
{
#define V(k)	\
	m[0].r[0]*T[0][0].r[k] + m[0].r[1]*T[0][1].r[k] + \
	m[1].r[0]*T[1][0].r[k] + m[1].r[1]*T[1][1].r[k]

    return vec2(V(0), V(1));

#undef V
}

// tensor times tensor ":" product
// (T1 : T2)_il = T1_ijk * T2_jkl
STATIC_INLINE void colonProdTen2x2x2Ten2x2x2 (Diderot_Mat2x2_t dst, Diderot_Ten2x2x2_t T1, Diderot_Ten2x2x2_t T2)
{
#define M(i,j)	\
	T1[i][0].r[0]*T2[0][0].r[j] + T1[i][0].r[1]*T2[0][1].r[j] + \
	T1[i][1].r[0]*T2[1][0].r[j] + T1[i][1].r[1]*T2[1][1].r[j]

    dst[0].v = vec2(M(0,0), M(0,1));
    dst[1].v = vec2(M(1,0), M(1,1));

#undef M
}

/********** 3x3x3 tensor functions **********/

STATIC_INLINE void copyTen3x3x3 (Diderot_Ten3x3x3_t dst, Diderot_Ten3x3x3_t src)
{
    dst[0][0].v = src[0][0].v;
    dst[0][1].v = src[0][1].v;
    dst[0][2].v = src[0][2].v;
    dst[1][0].v = src[1][0].v;
    dst[1][1].v = src[1][1].v;
    dst[1][2].v = src[1][2].v;
    dst[2][0].v = src[2][0].v;
    dst[2][1].v = src[2][1].v;
    dst[2][2].v = src[2][2].v;
}

// vector times tensor contraction
STATIC_INLINE void mulVec3Ten3x3x3 (Diderot_Mat3x3_t dst, Diderot_vec3_t v, Diderot_Ten3x3x3_t T)
{
    Diderot_union3_t u;
    u.v = v;
    dst[0].r[0] = u.r[0] * T[0][0].r[0] + u.r[1] * T[1][0].r[0] + u.r[2] * T[2][0].r[0];
    dst[0].r[1] = u.r[0] * T[0][0].r[1] + u.r[1] * T[1][0].r[1] + u.r[2] * T[2][0].r[1];
    dst[0].r[2] = u.r[0] * T[0][0].r[2] + u.r[1] * T[1][0].r[2] + u.r[2] * T[2][0].r[2];
    dst[1].r[0] = u.r[0] * T[0][1].r[0] + u.r[1] * T[1][1].r[0] + u.r[2] * T[2][1].r[0];
    dst[1].r[1] = u.r[0] * T[0][1].r[1] + u.r[1] * T[1][1].r[1] + u.r[2] * T[2][1].r[1];
    dst[1].r[2] = u.r[0] * T[0][1].r[2] + u.r[1] * T[1][1].r[2] + u.r[2] * T[2][1].r[2];
    dst[2].r[0] = u.r[0] * T[0][2].r[0] + u.r[1] * T[1][2].r[0] + u.r[2] * T[2][2].r[0];
    dst[2].r[1] = u.r[0] * T[0][2].r[1] + u.r[1] * T[1][2].r[1] + u.r[2] * T[2][2].r[1];
    dst[2].r[2] = u.r[0] * T[0][2].r[2] + u.r[1] * T[1][2].r[2] + u.r[2] * T[2][2].r[2];
}

// tensor times vector contraction
STATIC_INLINE void mulTen3x3x3Vec3 (Diderot_Mat3x3_t dst, Diderot_Ten3x3x3_t T, Diderot_vec3_t v)
{
    dst[0].v = vec3(dot3(T[0][0].v, v), dot3(T[0][1].v, v), dot3(T[0][2].v, v));
    dst[1].v = vec3(dot3(T[1][0].v, v), dot3(T[1][1].v, v), dot3(T[1][2].v, v));
    dst[2].v = vec3(dot3(T[2][0].v, v), dot3(T[2][1].v, v), dot3(T[2][2].v, v));
}

// The Frobenius norm of a tensor
STATIC_INLINE Diderot_real_t normTen3x3x3 (Diderot_Ten3x3x3_t T)
{
    Diderot_real_t sumSq =
	dot3(T[0][0].v,T[0][0].v) + dot3(T[0][1].v,T[0][1].v) + dot3(T[0][2].v,T[0][2].v) +
	dot3(T[1][0].v,T[1][0].v) + dot3(T[1][1].v,T[1][1].v) + dot3(T[1][2].v,T[1][2].v) +
	dot3(T[2][0].v,T[2][0].v) + dot3(T[2][1].v,T[2][1].v) + dot3(T[2][2].v,T[2][2].v);
    return SQRT(sumSq);
}

// tensor times matrix ":" product
// (T : M)_i = T_ijk * M_jk
STATIC_INLINE Diderot_vec3_t colonProdTen3x3x3Mat3x3 (Diderot_Ten3x3x3_t T, Diderot_Mat3x3_t m)
{
#define V(i)	dot3(T[i][0].v, m[0].v) + dot3(T[i][1].v, m[1].v) + dot3(T[i][2].v, m[2].v)

    return vec3(V(0), V(1), V(2));

#undef V
}

// matrix times tensor ":" product
// (M : T)_k = M_ij * T_ijk
STATIC_INLINE Diderot_vec3_t colonProdMat3x3Ten3x3x3 (Diderot_Mat3x3_t m, Diderot_Ten3x3x3_t T)
{
#define V(k)	\
	m[0].r[0]*T[0][0].r[k] + m[0].r[1]*T[0][1].r[k] + m[0].r[2]*T[0][2].r[k] + \
	m[1].r[0]*T[1][0].r[k] + m[1].r[1]*T[1][1].r[k] + m[1].r[2]*T[1][2].r[k] + \
	m[2].r[0]*T[0][0].r[k] + m[2].r[1]*T[2][1].r[k] + m[2].r[2]*T[2][2].r[k]

    return vec3(V(0), V(1), V(2));

#undef V
}

// tensor times tensor ":" product
// (T1 : T2)_il = T1_ijk * T2_jkl
STATIC_INLINE void colonProdTen3x3x3Ten3x3x3 (Diderot_Mat3x3_t dst, Diderot_Ten3x3x3_t T1, Diderot_Ten3x3x3_t T2)
{
#define M(i,j)	\
	T1[i][0].r[0]*T2[0][0].r[j] + T1[i][0].r[1]*T2[0][1].r[j] + T1[i][0].r[2]*T2[0][2].r[j] + \
	T1[i][1].r[0]*T2[1][0].r[j] + T1[i][1].r[1]*T2[1][1].r[j] + T1[i][1].r[2]*T2[1][2].r[j] + \
	T1[i][2].r[0]*T2[2][0].r[j] + T1[i][2].r[1]*T2[2][1].r[j] + T1[i][2].r[2]*T2[2][2].r[j]

    dst[0].v = vec3(M(0,0), M(0,1), M(0,2));
    dst[1].v = vec3(M(1,0), M(1,1), M(1,2));
    dst[2].v = vec3(M(2,0), M(2,1), M(2,2));

#undef M
}


/********** 4x4x4 tensor functions **********/

STATIC_INLINE void copyTen4x4x4 (Diderot_Ten4x4x4_t dst, Diderot_Ten4x4x4_t src)
{
    for (int i = 0;  i < 4;  i++) {
	dst[i][0].v = src[i][0].v;
	dst[i][1].v = src[i][1].v;
	dst[i][2].v = src[i][2].v;
	dst[i][3].v = src[i][3].v;
    }
}

// vector times tensor contraction
STATIC_INLINE void mulVec4Ten4x4x4 (Diderot_Mat4x4_t dst, Diderot_vec4_t v, Diderot_Ten4x4x4_t T)
{
    Diderot_union4_t u;
    u.v = v;
    dst[0].r[0] = u.r[0] * T[0][0].r[0] + u.r[1] * T[1][0].r[0] + u.r[2] * T[2][0].r[0] + u.r[3] * T[3][0].r[0];
    dst[0].r[1] = u.r[0] * T[0][0].r[1] + u.r[1] * T[1][0].r[1] + u.r[2] * T[2][0].r[1] + u.r[3] * T[3][0].r[1];
    dst[0].r[2] = u.r[0] * T[0][0].r[2] + u.r[1] * T[1][0].r[2] + u.r[2] * T[2][0].r[2] + u.r[3] * T[3][0].r[2];
    dst[0].r[3] = u.r[0] * T[0][0].r[3] + u.r[1] * T[1][0].r[3] + u.r[2] * T[2][0].r[3] + u.r[3] * T[3][0].r[3];
    dst[1].r[0] = u.r[0] * T[0][1].r[0] + u.r[1] * T[1][1].r[0] + u.r[2] * T[2][1].r[0] + u.r[3] * T[3][1].r[0];
    dst[1].r[1] = u.r[0] * T[0][1].r[1] + u.r[1] * T[1][1].r[1] + u.r[2] * T[2][1].r[1] + u.r[3] * T[3][1].r[1];
    dst[1].r[2] = u.r[0] * T[0][1].r[2] + u.r[1] * T[1][1].r[2] + u.r[2] * T[2][1].r[2] + u.r[3] * T[3][1].r[2];
    dst[1].r[3] = u.r[0] * T[0][1].r[3] + u.r[1] * T[1][1].r[3] + u.r[2] * T[2][1].r[3] + u.r[3] * T[3][1].r[3];
    dst[2].r[0] = u.r[0] * T[0][2].r[0] + u.r[1] * T[1][2].r[0] + u.r[2] * T[2][2].r[0] + u.r[3] * T[3][2].r[0];
    dst[2].r[1] = u.r[0] * T[0][2].r[1] + u.r[1] * T[1][2].r[1] + u.r[2] * T[2][2].r[1] + u.r[3] * T[3][2].r[1];
    dst[2].r[2] = u.r[0] * T[0][2].r[2] + u.r[1] * T[1][2].r[2] + u.r[2] * T[2][2].r[2] + u.r[3] * T[3][2].r[2];
    dst[2].r[3] = u.r[0] * T[0][2].r[3] + u.r[1] * T[1][2].r[3] + u.r[2] * T[2][2].r[3] + u.r[3] * T[3][2].r[3];
    dst[3].r[0] = u.r[0] * T[0][3].r[0] + u.r[1] * T[1][3].r[0] + u.r[2] * T[2][3].r[0] + u.r[3] * T[3][3].r[0];
    dst[3].r[1] = u.r[0] * T[0][3].r[1] + u.r[1] * T[1][3].r[1] + u.r[2] * T[2][3].r[1] + u.r[3] * T[3][3].r[1];
    dst[3].r[2] = u.r[0] * T[0][3].r[2] + u.r[1] * T[1][3].r[2] + u.r[2] * T[2][3].r[2] + u.r[3] * T[3][3].r[2];
    dst[3].r[3] = u.r[0] * T[0][3].r[3] + u.r[1] * T[1][3].r[3] + u.r[2] * T[2][3].r[3] + u.r[3] * T[3][3].r[3];
}

// tensor times vector contraction
STATIC_INLINE void mulTen4x4x4Vec4 (Diderot_Mat4x4_t dst, Diderot_Ten4x4x4_t T, Diderot_vec4_t v)
{
    dst[0].v = vec4(dot4(T[0][0].v, v), dot4(T[0][1].v, v), dot4(T[0][2].v, v), dot4(T[0][3].v, v));
    dst[1].v = vec4(dot4(T[1][0].v, v), dot4(T[1][1].v, v), dot4(T[1][2].v, v), dot4(T[1][3].v, v));
    dst[2].v = vec4(dot4(T[2][0].v, v), dot4(T[2][1].v, v), dot4(T[2][2].v, v), dot4(T[2][3].v, v));
    dst[3].v = vec4(dot4(T[3][0].v, v), dot4(T[3][1].v, v), dot4(T[3][2].v, v), dot4(T[3][3].v, v));
}

// The Frobenius norm of a tensor
STATIC_INLINE Diderot_real_t normTen4x4x4 (Diderot_Ten4x4x4_t T)
{
    Diderot_real_t sumSq =
	dot4(T[0][0].v,T[0][0].v) + dot4(T[0][1].v,T[0][1].v) + dot4(T[0][2].v,T[0][2].v) + dot4(T[0][3].v,T[0][3].v) +
	dot4(T[1][0].v,T[1][0].v) + dot4(T[1][1].v,T[1][1].v) + dot4(T[1][2].v,T[1][2].v) + dot4(T[1][3].v,T[1][3].v) +
	dot4(T[2][0].v,T[2][0].v) + dot4(T[2][1].v,T[2][1].v) + dot4(T[2][2].v,T[2][2].v) + dot4(T[2][3].v,T[2][3].v) +
	dot4(T[2][0].v,T[2][0].v) + dot4(T[2][1].v,T[2][1].v) + dot4(T[2][2].v,T[2][2].v) + dot4(T[3][3].v,T[3][3].v);
    return SQRT(sumSq);
}

#endif /*! _DIDEROT_INLINE_TENSOR_H_ */

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0