/*************************************************************************************

Grid physics library, www.github.com/paboyle/Grid 

Source file: lib/qcd/spin/Dirac.h

Copyright (C) 2015
Copyright (C) 2016

Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/*  END LEGAL */
#ifndef GRID_QCD_DIRAC_H
#define GRID_QCD_DIRAC_H

// code generated by the Mathematica notebook gamma-gen/gamma-gen.nb
#include <array>
#include <Grid/qcd/spin/Gamma.h>

namespace Grid {
  
inline QCD::Gamma adj (const QCD::Gamma & g)
{
   return QCD::Gamma (QCD::Gamma::adj[g.g]);
}

namespace QCD {

inline Gamma operator*(const Gamma & g1, const Gamma & g2)
{
   return Gamma (Gamma::mul[g1.g][g2.g]);
}

// FIXME
//
// Optimisation; switch over to a "multGammaX(ret._internal,arg._internal)" style early and
// note that doing so from the lattice operator will avoid copy back and case switch overhead, as
// was done for the tensor math operator to remove operator * notation early
//

//left multiply
template<class vtype> 
inline auto operator*(const Gamma &G, const iScalar<vtype> &arg)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
{
  iScalar<vtype> ret;
  ret._internal=G*arg._internal;
  return ret;
}

template<class vtype,int N>
inline auto operator*(const Gamma &G, const iVector<vtype, N> &arg)
->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type 
{
  iVector<vtype,N> ret;
  for(int i=0;i<N;i++){
    ret._internal[i]=G*arg._internal[i];
  }
  return ret;
}

template<class vtype, int N>
inline auto operator*(const Gamma &G, const iMatrix<vtype, N> &arg) 
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
{
  iMatrix<vtype,N> ret;
  for(int i=0;i<N;i++){
  for(int j=0;j<N;j++){
    ret._internal[i][j]=G*arg._internal[i][j];
  }}
  return ret;
}

//right multiply
template<class vtype>
inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
{
  iScalar<vtype> ret;
  ret._internal=arg._internal*G;
  return ret;
}

template<class vtype, int N>
inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G) 
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
{
  iMatrix<vtype,N> ret;
  for(int i=0;i<N;i++){
  for(int j=0;j<N;j++){
    ret._internal[i][j]=arg._internal[i][j]*G;
  }}
  return ret;
}

}}

#endif