Tensor_determinant.h

Go to the documentation of this file.

/*************************************************************************************
 
    Grid physics library, www.github.com/paboyle/Grid 
 
    Source file: ./lib/tensors/Tensor_determinant.h
 
    Copyright (C) 2015
 
Author: neo <cossu@post.kek.jp>
 
    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_MATH_DET_H
#define GRID_MATH_DET_H
 
NAMESPACE_BEGIN(Grid);

// Determinant function for scalar, vector, matrix
accelerator_inline ComplexF Determinant( const ComplexF &arg){    return arg;}
accelerator_inline ComplexD Determinant( const ComplexD &arg){    return arg;}
accelerator_inline RealF Determinant( const RealF &arg){    return arg;}
accelerator_inline RealD Determinant( const RealD &arg){    return arg;}
 
template<class vtype> accelerator_inline auto Determinant(const iScalar<vtype>&r) -> iScalar<decltype(Determinant(r._internal))>
{
  iScalar<decltype(Determinant(r._internal))> ret;
  ret._internal = Determinant(r._internal);
  return ret;
}
 
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr> 
accelerator_inline iScalar<vtype> Determinant(const iMatrix<vtype,N> &arg)
{
  iMatrix<vtype,N> ret(arg);
  iScalar<vtype> det = vtype(1.0);
  /* Conversion of matrix to upper triangular */
  for(int i = 0; i < N; i++){
    for(int j = 0; j < N; j++){
      if(j>i){
    vtype ratio = ret._internal[j][i]/ret._internal[i][i];
    for(int k = 0; k < N; k++){
      ret._internal[j][k] -= ratio * ret._internal[i][k];
    }
      }
    }
  }      
 
  for(int i = 0; i < N; i++)
    det *= ret._internal[i][i];   
 
  return det;
}
 
NAMESPACE_END(Grid);
 
#endif