mirror of
				https://github.com/paboyle/Grid.git
				synced 2025-11-04 05:54:32 +00:00 
			
		
		
		
	Adding several iMatrix utilities
This commit is contained in:
		@@ -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/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.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/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/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.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_extract_merge.h ./tensors/Tensor_inner.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=./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_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/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=./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
 | 
			
		||||
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
 | 
			
		||||
 
 | 
			
		||||
@@ -63,26 +63,70 @@ namespace Grid {
 | 
			
		||||
      for(int c1=0;c1<N;c1++){
 | 
			
		||||
	nrm = 0.0; 
 | 
			
		||||
	for(int c2=0;c2<N;c2++)
 | 
			
		||||
	  nrm = real(innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]));
 | 
			
		||||
	  nrm += real(innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]));
 | 
			
		||||
	nrm = 1.0/sqrt(nrm);
 | 
			
		||||
	std::cout << "norm : "<< nrm << "\n";
 | 
			
		||||
	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;
 | 
			
		||||
	  for(int c=0; c<N; ++c)
 | 
			
		||||
	    pr += ret._internal[c1][c]*ret._internal[b][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];
 | 
			
		||||
	  }
 | 
			
		||||
	}
 | 
			
		||||
	  
 | 
			
		||||
      }
 | 
			
		||||
      // assuming the determinant is ok
 | 
			
		||||
      return ret;
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  /////////////////////////////////////////////// 
 | 
			
		||||
  // 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 auto Determinant(const iMatrix<vtype,N> &arg)-> iScalar<decltype(Determinant(arg._internal[0][0]))>
 | 
			
		||||
    {
 | 
			
		||||
      iMatrix<vtype,N> ret(arg);
 | 
			
		||||
      iScalar<decltype(Determinant(arg._internal[0][0]))> det = 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
 | 
			
		||||
  /////////////////////////////////////////////// 
 | 
			
		||||
 
 | 
			
		||||
@@ -1,9 +1,5 @@
 | 
			
		||||
 | 
			
		||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_gamma Test_main Test_multishift_sqrt Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
 | 
			
		||||
Test_GaugeAction_LDADD=-lGrid
 | 
			
		||||
bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_gamma Test_GaugeAction Test_main Test_multishift_sqrt Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
Test_cayley_cg_SOURCES=Test_cayley_cg.cc
 | 
			
		||||
@@ -74,6 +70,10 @@ Test_gamma_SOURCES=Test_gamma.cc
 | 
			
		||||
Test_gamma_LDADD=-lGrid
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
 | 
			
		||||
Test_GaugeAction_LDADD=-lGrid
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
Test_main_SOURCES=Test_main.cc
 | 
			
		||||
Test_main_LDADD=-lGrid
 | 
			
		||||
 | 
			
		||||
 
 | 
			
		||||
@@ -216,16 +216,22 @@ int main (int argc, char ** argv)
 | 
			
		||||
    scm=transposeIndex<1>(scm);
 | 
			
		||||
    
 | 
			
		||||
   
 | 
			
		||||
    //random(SerialRNG, cm);
 | 
			
		||||
    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 << "  "<< tracecm << std::endl;
 | 
			
		||||
 | 
			
		||||
    cm = ProjectOnGroup(cm);
 | 
			
		||||
    std::cout << cm << "  " << std::endl;
 | 
			
		||||
 | 
			
		||||
    TComplex det = Determinant(cm);
 | 
			
		||||
    
 | 
			
		||||
    std::cout << "determinant: " << det <<  std::endl;
 | 
			
		||||
    cm = Exponentiate(cm, 1.0, 10);
 | 
			
		||||
    std::cout << cm << "  " << std::endl;
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
//    Foo = Foo+scalar; // LatticeColourMatrix+Scalar
 | 
			
		||||
//    Foo = Foo*scalar; // LatticeColourMatrix*Scalar
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user