#ifndef _matrix_algebra_3x3_hpp_
#define _matrix_algebra_3x3_hpp_

//@HEADER
// ************************************************************************
//
// MiniFE: Simple Finite Element Assembly and Solve
// Copyright (2006-2013) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
//
// ************************************************************************
//@HEADER

#ifndef KERNEL_PREFIX
#define KERNEL_PREFIX
#endif

namespace miniFE {

template<typename Scalar>
#ifdef __CUDACC__
  __host__ __device__
#endif
KERNEL_PREFIX void fill(Scalar* begin, Scalar* end, const Scalar& val)
{
  while(begin != end) {*begin++ = val;}
}

template<typename Scalar>
KERNEL_PREFIX void inverse_and_determinant3x3(const Scalar* J, Scalar* invJ, Scalar& detJ)
{
  //hardwired "3x3" in function-name allows us to assume
  //that J and invJ have length 9:

  Scalar J00 = J[0];
  Scalar J01 = J[1];
  Scalar J02 = J[2];

  Scalar J10 = J[3];
  Scalar J11 = J[4];
  Scalar J12 = J[5];

  Scalar J20 = J[6];
  Scalar J21 = J[7];
  Scalar J22 = J[8];

  Scalar term0 = J22*J11 - J21*J12;
  Scalar term1 = J22*J01 - J21*J02;
  Scalar term2 = J12*J01 - J11*J02;

  detJ = J00*term0 - J10*term1 + J20*term2;

  Scalar inv_detJ = 1.0/detJ;

  invJ[0] =  term0*inv_detJ;
  invJ[1] = -term1*inv_detJ;
  invJ[2] =  term2*inv_detJ;

  invJ[3] = -(J22*J10 - J20*J12)*inv_detJ;
  invJ[4] =  (J22*J00 - J20*J02)*inv_detJ;
  invJ[5] = -(J12*J00 - J10*J02)*inv_detJ;

  invJ[6] =  (J21*J10 - J20*J11)*inv_detJ;
  invJ[7] = -(J21*J00 - J20*J01)*inv_detJ;
  invJ[8] =  (J11*J00 - J10*J01)*inv_detJ;
}

template<typename Scalar>
KERNEL_PREFIX void matmat3x3(const Scalar* A, const Scalar* B, Scalar* C)
{
  //hardwired "3x3" in function-name allows us to assume args have length 9:
  //A,B,C are all assumed to be ordered such that columns are contiguous.

  const Scalar zero = 0;
  miniFE::fill(C, C+9, zero);

  for(int i=0; i<3; ++i) {
    for(int j=0; j<3; ++j) {
      C[i+j*3] = A[i+0]*B[j*3+0]
               + A[i+3]*B[j*3+1]
               + A[i+6]*B[j*3+2];
    }
  }
}

template<typename Scalar>
KERNEL_PREFIX Scalar determinant3x3(const Scalar* J)
{
  //hardwired "3x3" in function-name allows us to assume that J has length 9:

  Scalar J00 = J[0];
  Scalar J01 = J[1];
  Scalar J02 = J[2];

  Scalar J10 = J[3];
  Scalar J11 = J[4];
  Scalar J12 = J[5];

  Scalar J20 = J[6];
  Scalar J21 = J[7];
  Scalar J22 = J[8];

  Scalar term0 = J22*J11 - J21*J12;
  Scalar term1 = J22*J01 - J21*J02;
  Scalar term2 = J12*J01 - J11*J02;

  Scalar detJ = J00*term0 - J10*term1 + J20*term2;

  return detJ;
}

template<typename Scalar>
KERNEL_PREFIX void matmat3x3_X_3xn(const Scalar* A, int n, const Scalar* B, Scalar* C)
{
  //A is 3x3, B is 3xn. So C is also 3xn.
  //A,B,C are all assumed to be ordered such that columns are contiguous.

  Scalar* Cj = C;
  const Scalar* Bj = B;
  for(int j=0; j<n; ++j) {
    Cj[0] = A[0]*Bj[0] + A[3]*Bj[1] + A[6]*Bj[2];
    Cj[1] = A[1]*Bj[0] + A[4]*Bj[1] + A[7]*Bj[2];
    Cj[2] = A[2]*Bj[0] + A[5]*Bj[1] + A[8]*Bj[2];
    Bj += 3;
    Cj += 3;
  }
}

template<typename Scalar>
KERNEL_PREFIX void matTransMat3x3_X_3xn(const Scalar* A, int n, const Scalar* B, Scalar* C)
{
  //A is 3x3, B is 3xn. So C is also 3xn.
  //A,B,C are all assumed to be ordered such that columns are contiguous.

  Scalar* Cj = C;
  const Scalar* Bj = B;
  for(int j=0; j<n; ++j) {
    Cj[0] = A[0]*Bj[0] + A[1]*Bj[1] + A[2]*Bj[2];
    Cj[1] = A[3]*Bj[0] + A[4]*Bj[1] + A[5]*Bj[2];
    Cj[2] = A[6]*Bj[0] + A[7]*Bj[1] + A[8]*Bj[2];
    Bj += 3;
    Cj += 3;
  }
}

}//namespace miniFE

#endif