diff --git a/lib/Make.inc b/lib/Make.inc index a660a775..dab78b37 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_poke.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_peek.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Threads.h ./Grid.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h +HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_poke.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_peek.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Threads.h ./Grid.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h CCFILES=./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/Simd.h b/lib/Simd.h index f93c3070..7da0d7ad 100644 --- a/lib/Simd.h +++ b/lib/Simd.h @@ -46,6 +46,18 @@ namespace Grid { inline ComplexF innerProduct(const ComplexF & l, const ComplexF & r) { return conjugate(l)*r; } inline RealD innerProduct(const RealD & l, const RealD & r) { return l*r; } inline RealF innerProduct(const RealF & l, const RealF & r) { return l*r; } + + inline ComplexD Reduce(const ComplexD& r){ return r; } + inline ComplexF Reduce(const ComplexF& r){ return r; } + inline RealD Reduce(const RealD& r){ return r; } + inline RealF Reduce(const RealF& r){ return r; } + + inline RealD toReal(const ComplexD& r){ return real(r); } + inline RealF toReal(const ComplexF& r){ return real(r); } + inline RealD toReal(const RealD& r){ return r; } + inline RealF toReal(const RealF& r){ return r; } + + //////////////////////////////////////////////////////////////////////////////// //Provide support functions for basic real and complex data types required by Grid diff --git a/lib/Tensors.h b/lib/Tensors.h index a2246f4a..90c9d7a8 100644 --- a/lib/Tensors.h +++ b/lib/Tensors.h @@ -9,6 +9,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/lib/lattice/Lattice_ET.h b/lib/lattice/Lattice_ET.h index caa45e69..82d11218 100644 --- a/lib/lattice/Lattice_ET.h +++ b/lib/lattice/Lattice_ET.h @@ -1,4 +1,3 @@ - #ifndef GRID_LATTICE_ET_H #define GRID_LATTICE_ET_H diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h index 91b39007..6d72ef31 100644 --- a/lib/lattice/Lattice_unary.h +++ b/lib/lattice/Lattice_unary.h @@ -24,5 +24,18 @@ PARALLEL_FOR_LOOP return ret; } + template Lattice expMat(const Lattice &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){ + Lattice ret(rhs._grid); + ret.checkerboard = rhs.checkerboard; + conformable(ret,rhs); +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss++){ + ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); + } + return ret; + } + + + } #endif diff --git a/lib/tensors/Tensor_Ta.h b/lib/tensors/Tensor_Ta.h index 0eea1e6d..e49ede35 100644 --- a/lib/tensors/Tensor_Ta.h +++ b/lib/tensors/Tensor_Ta.h @@ -1,5 +1,7 @@ #ifndef GRID_MATH_TA_H #define GRID_MATH_TA_H + + namespace Grid { /////////////////////////////////////////////// @@ -36,7 +38,8 @@ namespace Grid { /////////////////////////////////////////////// - // ProjectOnGroup function for scalar, vector, matrix + // ProjectOnGroup function for scalar, vector, matrix + // Projects on orthogonal, unitary group /////////////////////////////////////////////// @@ -59,22 +62,23 @@ namespace Grid { { // need a check for the group type? iMatrix ret(arg); - double nrm; + RealD nrm; + vtype inner; for(int c1=0;c1 inline auto Determinant(const iScalar&r) -> iScalar - { - iScalar ret; - ret._internal = Determinant(r._internal); - return ret; - } - - template::TensorLevel == 0 >::type * =nullptr> - inline iScalar Determinant(const iMatrix &arg) - { - iMatrix ret(arg); - iScalar 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; - } - - - - - /////////////////////////////////////////////// - // Exponentiate function for scalar, vector, matrix - /////////////////////////////////////////////// - - - template inline iScalar Exponentiate(const iScalar&r, double alpha, int Nexp) - { - iScalar ret; - ret._internal = Exponentiate(r._internal, alpha, Nexp); - return ret; - } - - - template::TensorLevel == 0 >::type * =nullptr> - inline iMatrix Exponentiate(const iMatrix &arg, double alpha, int Nexp) - { - iMatrix unit(1.0); - iMatrix temp(unit); - - for(int i=Nexp; i>=1;--i){ - temp *= alpha/double(i); - temp = unit + temp*arg; - } - return ProjectOnGroup(temp); - } - - } #endif diff --git a/lib/tensors/Tensor_determinant.h b/lib/tensors/Tensor_determinant.h new file mode 100644 index 00000000..7c8af9ed --- /dev/null +++ b/lib/tensors/Tensor_determinant.h @@ -0,0 +1,45 @@ +#ifndef GRID_MATH_DET_H +#define GRID_MATH_DET_H +namespace Grid { + /////////////////////////////////////////////// + // Determinant function for scalar, vector, matrix + /////////////////////////////////////////////// + inline ComplexF Determinant( const ComplexF &arg){ return arg;} + inline ComplexD Determinant( const ComplexD &arg){ return arg;} + inline RealF Determinant( const RealF &arg){ return arg;} + inline RealD Determinant( const RealD &arg){ return arg;} + + template inline auto Determinant(const iScalar&r) -> iScalar + { + iScalar ret; + ret._internal = Determinant(r._internal); + return ret; + } + + template::TensorLevel == 0 >::type * =nullptr> + inline iScalar Determinant(const iMatrix &arg) + { + iMatrix ret(arg); + iScalar 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; + } + + + +} +#endif diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h new file mode 100644 index 00000000..f8437116 --- /dev/null +++ b/lib/tensors/Tensor_exp.h @@ -0,0 +1,37 @@ +#ifndef GRID_MATH_EXP_H +#define GRID_MATH_EXP_H + +#define DEFAULT_MAT_EXP 12 + +namespace Grid { + + /////////////////////////////////////////////// + // Exponentiate function for scalar, vector, matrix + /////////////////////////////////////////////// + + + template inline iScalar Exponentiate(const iScalar&r, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP) + { + iScalar ret; + ret._internal = Exponentiate(r._internal, alpha, Nexp); + return ret; + } + + + template::TensorLevel == 0 >::type * =nullptr> + inline iMatrix Exponentiate(const iMatrix &arg, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP ) + { + iMatrix unit(1.0); + iMatrix temp(unit); + + for(int i=Nexp; i>=1;--i){ + temp *= alpha/ComplexD(i); + temp = unit + temp*arg; + } + return ProjectOnGroup(temp);//maybe not strictly necessary + } + + + +} +#endif diff --git a/tests/Test_main.cc b/tests/Test_main.cc index 9db38059..8eecf0c3 100644 --- a/tests/Test_main.cc +++ b/tests/Test_main.cc @@ -215,21 +215,21 @@ int main (int argc, char ** argv) random(SerialRNG, cm); std::cout << cm << std::endl; - //cm = Ta(cm); - //TComplex tracecm= trace(cm); - //std::cout << cm << std::endl; + cm = Ta(cm); + TComplex tracecm= trace(cm); + std::cout << cm << std::endl; - cm = ProjectOnGroup(cm); + + cm = Exponentiate(cm, 1.0, 12); std::cout << cm << " " << std::endl; - cm = ProjectOnGroup(cm); - std::cout << cm << " " << std::endl; - - - TComplex det = Determinant(cm); - + Complex det = Determinant(cm); std::cout << "determinant: " << det << std::endl; - cm = Exponentiate(cm, 1.0, 10); + + cm = ProjectOnGroup(cm); std::cout << cm << " " << std::endl; + cm = ProjectOnGroup(cm); + std::cout << cm << " " << std::endl; + det = Determinant(cm); std::cout << "determinant: " << det << std::endl; @@ -244,6 +244,9 @@ int main (int argc, char ** argv) LatticeComplex trscMat(&Fine); trscMat = trace(scMat); // Trace + // Exponentiate test + cMat = expMat(cMat, ComplexD(1.0, 0.0)); + // LatticeComplex trlcMat(&Fine); // trlcMat = trace(lcMat); // Trace involving iVector - now generates error