portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-03 21:44:33 +00:00
Code Releases Activity

Cleaning up the recursion for traceIndex<n> after the changes the enable G++ to

Browse Source 
compile it again.
This commit is contained in:
Peter Boyle
2015-07-01 23:43:57 +01:00
parent 69a3d3203c
commit 84ec7c40cd
7 changed files with 331 additions and 329 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/Make.inc
View File

@@ -1,4 +1,4 @@
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h
CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

0
lib/tensors/Tensor_peek.h → lib/Old/Tensor_peek.h
View File

0
lib/tensors/Tensor_poke.h → lib/Old/Tensor_poke.h
View File

5
lib/Tensors.h
View File

@@ -8,11 +8,12 @@
#include <tensors/Tensor_outer.h>
#include <tensors/Tensor_transpose.h>
#include <tensors/Tensor_trace.h>
#include <tensors/Tensor_index.h>
#include <tensors/Tensor_Ta.h>
#include <tensors/Tensor_determinant.h>
#include <tensors/Tensor_exp.h>
#include <tensors/Tensor_peek.h>
#include <tensors/Tensor_poke.h>
//#include <tensors/Tensor_peek.h>
//#include <tensors/Tensor_poke.h>
#include <tensors/Tensor_reality.h>
#include <tensors/Tensor_unary.h>
#include <tensors/Tensor_extract_merge.h>

321
lib/tensors/Tensor_index.h Normal file
View File

@@ -0,0 +1,321 @@
#ifndef GRID_TENSOR_INDEX_H
#define GRID_TENSOR_INDEX_H
////////////////////////////////////////////////////////////////////////////////////////
// Recursion for trace, transpose, peek, poke a specific index
////////////////////////////////////////////////////////////////////////////////////////
// Allow trace to recurse if vector, but never terminate on a vector
// trace of a different index can distribute across the vector index in a replicated way
// but we do not trace a vector index.
namespace Grid {
/* Needed?
template<int Level> inline ComplexF traceIndex(const ComplexF arg) { return arg;}
template<int Level> inline ComplexD traceIndex(const ComplexD arg) { return arg;}
template<int Level> inline RealF traceIndex(const RealF arg) { return arg;}
template<int Level> inline RealD traceIndex(const RealD arg) { return arg;}
*/
template<int Level>
class TensorIndexRecursion {
public:
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))> ret;
ret._internal = TensorIndexRecursion<Level-1>::traceIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i][j]);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for peeking a specific index
////////////////////////////////////////////
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i);
return ret;
}
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i,int j) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i,j);
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii,int jj) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii,jj);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii);
}}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii,jj);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for poking a specific index
////////////////////////////////////////////
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0))> &arg, int i)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i);
}
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0,0))> &arg, int i,int j)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i,j);
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i);
}
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i,j);
}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i);
}}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i,j);
}}
}
////////////////////////////////////////////
// Recursion for transposing a specific index
////////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iVector<vtype,N> arg) -> iVector<vtype,N>
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i][j]);
}}
return ret;
}
};
////////////////////////////
// strip const & ref quali's
////////////////////////////
#define RemoveCRV(a) typename std::remove_const<typename std::remove_reference<decltype(a)>::type>::type
template<>
class TensorIndexRecursion<0> {
public:
/////////////////////////////////////////
// Ends recursion for trace (scalar/vector/matrix)
/////////////////////////////////////////
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<RemoveCRV(arg._internal)>
{
iScalar<RemoveCRV(arg._internal)> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0])>
{
iScalar<RemoveCRV(arg._internal[0])> ret;
ret._internal=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+ arg._internal[i];
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0][0])>
{
iScalar<RemoveCRV(arg._internal[0][0])> ret;
ret=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+arg._internal[i][i];
}
return ret;
}
/////////////////////////////////////////
// Ends recursion for transpose scalar/matrix ; no way to terminate on vector
/////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iMatrix<vtype,N> ret;
ret=zero;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = ret._internal[i][j]+arg._internal[i][j];
}}
return ret;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// End recursion for peeking a specific index; single index on vector, double index on matrix
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii];
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii][jj];
return ret;
}
// Vector poke, one index
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iScalar<vtype> &arg,int i)
{
ret._internal[i] = arg._internal;
}
// Matrix poke two indices
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iScalar<vtype> &arg,int i,int j)
{
ret._internal[i][j] = arg._internal;
}
};
////////////////////////////////////////////////////////////////////////////////////////////////////////
// External wrappers
////////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Level,class vtype> inline auto traceIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::traceIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto transposeIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::transposeIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i,int j) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i,j);
return ret;
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0)) &arg,int i)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i);
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0,0)) &arg,int i,int j)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i,j);
}
#undef RemoveCRV
}
#endif

330
lib/tensors/Tensor_trace.h
View File

@@ -1,8 +1,10 @@
#ifndef GRID_MATH_TRACE_H
#define GRID_MATH_TRACE_H
namespace Grid {
//////////////////////////////////////////////////////////////////
// Traces: both all indices and a specific index
// Traces: both all indices and a specific index. Indices must be
// either scalar or matrix
/////////////////////////////////////////////////////////////////
inline ComplexF trace( const ComplexF &arg){ return arg;}
@@ -10,11 +12,6 @@ inline ComplexD trace( const ComplexD &arg){ return arg;}
inline RealF trace( const RealF &arg){ return arg;}
inline RealD trace( const RealD &arg){ return arg;}
template<int Level> inline ComplexF traceIndex(const ComplexF arg) { return arg;}
template<int Level> inline ComplexD traceIndex(const ComplexD arg) { return arg;}
template<int Level> inline RealF traceIndex(const RealF arg) { return arg;}
template<int Level> inline RealD traceIndex(const RealD arg) { return arg;}
template<class vtype,int N>
inline auto trace(const iMatrix<vtype,N> &arg) -> iScalar<decltype(trace(arg._internal[0][0]))>
{
@@ -25,6 +22,7 @@ inline auto trace(const iMatrix<vtype,N> &arg) -> iScalar<decltype(trace(arg._in
}
return ret;
}
template<class vtype>
inline auto trace(const iScalar<vtype> &arg) -> iScalar<decltype(trace(arg._internal))>
{
@@ -33,326 +31,6 @@ inline auto trace(const iScalar<vtype> &arg) -> iScalar<decltype(trace(arg._inte
return ret;
}
/////////////////////////////////////
// Alternate implementation .. count
// indices from left. Annoying but
// I'm working around compiler bugs in gcc and earlier icpc
/////////////////////////////////////
// Allow to recurse if vector, but never terminate on a vector
// ltrace of a different index can distribute across the vector index in a replicated way
// but we do not ltrace a vector index.
template<int Level>
class TensorIndexRecursion {
public:
////////////////////////////////////////////
// Recursion for tracing a specific index
////////////////////////////////////////////
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))> ret;
ret._internal = TensorIndexRecursion<Level-1>::traceIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i][j]);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for peeking a specific index
////////////////////////////////////////////
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i);
return ret;
}
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i,int j) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i,j);
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii,int jj) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii,jj);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii);
}}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii,jj);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for poking a specific index
////////////////////////////////////////////
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0))> &arg, int i)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i);
}
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0,0))> &arg, int i,int j)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i,j);
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i);
}
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i,j);
}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i);
}}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i,j);
}}
}
////////////////////////////////////////////
// Recursion for transposing a specific index
////////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iVector<vtype,N> arg) -> iVector<vtype,N>
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i][j]);
}}
return ret;
}
};
////////////////////////////
// strip const & ref quali's
////////////////////////////
#define RemoveCRV(a) typename std::remove_const<typename std::remove_reference<decltype(a)>::type>::type
template<>
class TensorIndexRecursion<0> {
public:
/////////////////////////////////////////
// Ends recursion for trace (scalar/vector/matrix)
/////////////////////////////////////////
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<RemoveCRV(arg._internal)>
{
iScalar<RemoveCRV(arg._internal)> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0])>
{
iScalar<RemoveCRV(arg._internal[0])> ret;
ret._internal=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+ arg._internal[i];
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0][0])>
{
iScalar<RemoveCRV(arg._internal[0][0])> ret;
ret=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+arg._internal[i][i];
}
return ret;
}
/////////////////////////////////////////
// Ends recursion for transpose scalar/vector/matrix
/////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iVector<vtype,N> arg) -> iVector<vtype,N>
{
iScalar<vtype> ret;
ret._internal=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+ arg._internal[i];
}
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iScalar<vtype> ret;
ret=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+arg._internal[i][i];
}
return ret;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// End recursion for peeking a specific index; single index on vector, double index on matrix
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii];
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii][jj];
return ret;
}
// Vector poke, one index
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iScalar<vtype> &arg,int i)
{
ret._internal[i] = arg._internal;
}
// Matrix poke two indices
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iScalar<vtype> &arg,int i,int j)
{
ret._internal[i][j] = arg._internal;
}
};
////////////////////////////////////////////////////////////////////////////////////////////////////////
// External wrappers
////////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Level,class vtype> inline auto traceIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::traceIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto transposeIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::transposeIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i,int j) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i,j);
return ret;
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0)) &arg,int i)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i);
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0,0)) &arg,int i,int j)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i,j);
}
#undef RemoveCRV
}
#endif

2
lib/tensors/Tensor_transpose.h
View File

@@ -59,6 +59,7 @@ template<class vtype>
// Transpose a specific index; instructive to compare this style of recursion termination
// to that of adj; which is easiers?
////////////////////////////////////////////////////////////////////////////////////////////
#if 0
template<int Level,class vtype,int N> inline
typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,Level>::value, iMatrix<vtype,N> >::type
transposeIndex (const iMatrix<vtype,N> &arg)
@@ -96,6 +97,7 @@ transposeIndex (const iScalar<vtype> &arg)
{
return arg;
}
#endif
}
#endif