mirror of
https://github.com/paboyle/Grid.git
synced 2025-10-31 03:54:33 +00:00
Lattice matrix exponential ok
This commit is contained in:
neo
@@ -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
|
||||
|
||||
12
lib/Simd.h
12
lib/Simd.h
@@ -47,6 +47,18 @@ namespace Grid {
|
||||
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
|
||||
//Single and double precision versions. Should be able to template this once only.
|
||||
|
||||
@@ -9,6 +9,8 @@
|
||||
#include <tensors/Tensor_transpose.h>
|
||||
#include <tensors/Tensor_trace.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_reality.h>
|
||||
|
||||
@@ -1,4 +1,3 @@
|
||||
|
||||
#ifndef GRID_LATTICE_ET_H
|
||||
#define GRID_LATTICE_ET_H
|
||||
|
||||
|
||||
@@ -24,5 +24,18 @@ PARALLEL_FOR_LOOP
|
||||
return ret;
|
||||
}
|
||||
|
||||
template<class obj> Lattice<obj> expMat(const Lattice<obj> &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){
|
||||
Lattice<obj> ret(rhs._grid);
|
||||
ret.checkerboard = rhs.checkerboard;
|
||||
conformable(ret,rhs);
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int ss=0;ss<rhs._grid->oSites();ss++){
|
||||
ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp);
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
@@ -1,5 +1,7 @@
|
||||
#ifndef GRID_MATH_TA_H
|
||||
#define GRID_MATH_TA_H
|
||||
|
||||
|
||||
namespace Grid {
|
||||
|
||||
///////////////////////////////////////////////
|
||||
@@ -37,6 +39,7 @@ namespace Grid {
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// 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<vtype,N> ret(arg);
|
||||
double nrm;
|
||||
RealD nrm;
|
||||
vtype inner;
|
||||
for(int c1=0;c1<N;c1++){
|
||||
nrm = 0.0;
|
||||
zeroit(inner);
|
||||
for(int c2=0;c2<N;c2++)
|
||||
nrm += real(innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]));
|
||||
nrm = 1.0/sqrt(nrm);
|
||||
std::cout << "norm : "<< nrm << "\n";
|
||||
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
|
||||
|
||||
nrm = 1.0/sqrt(Reduce(toReal(inner)));
|
||||
for(int c2=0;c2<N;c2++)
|
||||
ret._internal[c1][c2]*= nrm;
|
||||
|
||||
for (int b=c1+1; b<N; ++b){
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr = 0.0;
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
|
||||
zeroit(pr);
|
||||
for(int c=0; c<N; ++c)
|
||||
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
|
||||
|
||||
std::cout << "pr : "<< pr << "\n";
|
||||
for(int c=0; c<N; ++c){
|
||||
ret._internal[b][c] -= pr * ret._internal[c1][c];
|
||||
}
|
||||
@@ -86,74 +90,6 @@ 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<class vtype> 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>
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// Exponentiate function for scalar, vector, matrix
|
||||
///////////////////////////////////////////////
|
||||
|
||||
|
||||
template<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, double alpha, int Nexp)
|
||||
{
|
||||
iScalar<vtype> ret;
|
||||
ret._internal = Exponentiate(r._internal, alpha, Nexp);
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
|
||||
inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &arg, double alpha, int Nexp)
|
||||
{
|
||||
iMatrix<vtype,N> unit(1.0);
|
||||
iMatrix<vtype,N> temp(unit);
|
||||
|
||||
for(int i=Nexp; i>=1;--i){
|
||||
temp *= alpha/double(i);
|
||||
temp = unit + temp*arg;
|
||||
}
|
||||
return ProjectOnGroup(temp);
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
45
lib/tensors/Tensor_determinant.h
Normal file
45
lib/tensors/Tensor_determinant.h
Normal file
@@ -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<class vtype> 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>
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
37
lib/tensors/Tensor_exp.h
Normal file
37
lib/tensors/Tensor_exp.h
Normal file
@@ -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<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP)
|
||||
{
|
||||
iScalar<vtype> ret;
|
||||
ret._internal = Exponentiate(r._internal, alpha, Nexp);
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
|
||||
inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &arg, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP )
|
||||
{
|
||||
iMatrix<vtype,N> unit(1.0);
|
||||
iMatrix<vtype,N> 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
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user