diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0bf7a033e..474b90296 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,15 +1,41 @@ -types: +stages: + - SCA + - style - build - test + - deploy -GQMCT: - type: build +Ising: + stage: build script: - - . configure.sh - - cd Libraries - - make - - cd .. - - cd Prog_8 - - make - + - cmake -E make_directory build + - cd build + - cmake -G "Unix Makefiles" -DCMAKE_BUILD_TYPE=Release .. + - cmake --build . --target qmct_Hub --config Release + +Hubbard: + stage: build + script: + - cmake -E make_directory build + - cd build + - cmake -G "Unix Makefiles" -DCMAKE_BUILD_TYPE=Release .. + - cmake --build . --target qmct_Hub --config Release + +SPT: + stage: build + script: + - cmake -E make_directory build + - cd build + - cmake -G "Unix Makefiles" -DCMAKE_BUILD_TYPE=Release .. + - cmake --build . --target qmct_SPT --config Release + +Analysis: + stage: build + script: + - cmake -E make_directory build + - cd build + - cmake -G "Unix Makefiles" -DCMAKE_BUILD_TYPE=Release .. + - cmake --build . --target cov_eq --config Release + - cmake --build . --target cov_tau --config Release + - cmake --build . --target jackv5 --config Release diff --git a/Analysis_7/Compile_cov b/Analysis_7/Compile_cov deleted file mode 100644 index d3643a78f..000000000 --- a/Analysis_7/Compile_cov +++ /dev/null @@ -1,14 +0,0 @@ -TARGET= cov_tau.out -OBJS= cov_tau.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Analysis_7/Compile_en b/Analysis_7/Compile_en deleted file mode 100644 index 8767339d4..000000000 --- a/Analysis_7/Compile_en +++ /dev/null @@ -1,14 +0,0 @@ -TARGET= jackv5.out -OBJS= jackv5.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Analysis_7/Compile_eq b/Analysis_7/Compile_eq deleted file mode 100644 index 9613a2271..000000000 --- a/Analysis_7/Compile_eq +++ /dev/null @@ -1,14 +0,0 @@ -TARGET= cov_eq.out -OBJS= cov_eq.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Analysis_7/Makefile b/Analysis_7/Makefile deleted file mode 100644 index c92c401c8..000000000 --- a/Analysis_7/Makefile +++ /dev/null @@ -1,19 +0,0 @@ -FC= $(mpif90) -FC= $(f90) -FLAGS= $(FL) -LIBS= $(Libs)/Modules/modules_90.a \ - $(Libs)/MyEis/libeis.a \ - $(Libs)/MyNag/libnag.a \ - $(Libs)/MyLin/liblin.a \ - $(LIB_BLAS_LAPACK) - -all: - cp $(Libs)/Modules/*.mod . ;\ - (make -f Compile_en FC="$(FC)" FLAGS="$(FLAGS)" LIBS="$(LIBS)" ) ;\ - (make -f Compile_cov FC="$(FC)" FLAGS="$(FLAGS)" LIBS="$(LIBS)" ) ;\ - (make -f Compile_eq FC="$(FC)" FLAGS="$(FLAGS)" LIBS="$(LIBS)" ) -clean: - (make -f Compile_eq clean );\ - (make -f Compile_cov clean );\ - (make -f Compile_en clean );\ - rm *.mod *~ \#* *.out diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 000000000..172c23aea --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,91 @@ +# CMake project file for General QMCT + +################################################## +# Define the project and the depencies that it has +################################################## + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) +PROJECT(General_QMCT Fortran) + +# Set the General QMCT version +SET(VERSION 1.0) + +# Add our local modlues to the module path +SET(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules/") + +# Uncomment if it is required that Fortran 90 is supported +IF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) + MESSAGE(FATAL_ERROR "Fortran compiler does not support F90") +ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) + +# Set some options the user may choose +# Uncomment the below if you want the user to choose a parallelization library +OPTION(USE_MPI "Use the MPI library for parallelization" OFF) +#OPTION(USE_OPENMP "Use OpenMP for parallelization" OFF) + +# This INCLUDE statement executes code that sets the compile flags for DEBUG, +# RELEASE, and TESTING. You should review this file and make sure the flags +# are to your liking. +INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake) +# Locate and set parallelization libraries. There are some CMake peculiarities +# taken care of here, such as the fact that the FindOpenMP routine doesn't know +# about Fortran. +INCLUDE(${CMAKE_MODULE_PATH}/SetParallelizationLibrary.cmake) +# Setup the LAPACK libraries. This also takes care of peculiarities, such as +# the fact the searching for MKL requires a C compiler, and that the results +# are not stored in the cache. +INCLUDE(${CMAKE_MODULE_PATH}/SetUpLAPACK.cmake) + +# There is an error in CMAKE with this flag for pgf90. Unset it +GET_FILENAME_COMPONENT(FCNAME ${CMAKE_Fortran_COMPILER} NAME) +IF(FCNAME STREQUAL "pgf90") + UNSET(CMAKE_SHARED_LIBRARY_LINK_Fortran_FLAGS) +ENDIF(FCNAME STREQUAL "pgf90") + +############################################################ +# Define the actual files and folders that make up the build +############################################################ + +# Define the executable name +SET(PROGEXE qmct) +SET(BLA_STATIC TRUE) + +#SET(CMAKE_AR gcc-ar) +#SET(CMAKE_RANLIB gcc-ranlib) + +# Define some directories +SET(SRC ${CMAKE_SOURCE_DIR}/src) +SET(LIB ${CMAKE_SOURCE_DIR}/lib) +SET(BIN ${CMAKE_SOURCE_DIR}/bin) + +# Define the library names +SET(MODULES modules) +SET(MYEIS myeis) +SET(MYNAG mynag) +SET(MYLIN mylin) + +SET(SRCMODULES ${SRC}/Modules) +SET(SRCEIS ${SRC}/MyEis) +SET(SRCNAG ${SRC}/MyNag) +SET(SRCLIN ${SRC}/MyLin) +SET(SRCPROG ${SRC}/Prog) +SET(SRCANALYSIS ${SRC}/Analysis) + +# Have the .mod files placed in the lib folder +SET(CMAKE_Fortran_MODULE_DIRECTORY ${LIB}) + +# The sources for the libraries and have it placed in the lib folder +ADD_SUBDIRECTORY(${SRCMODULES} ${LIB}/Modules) +ADD_SUBDIRECTORY(${SRCEIS} ${LIB}/MyEis) + +ADD_SUBDIRECTORY(${SRCNAG} ${LIB}/MyNAG) +ADD_SUBDIRECTORY(${SRCLIN} ${LIB}/MyLin) + +# The source for the FOO binary and have it placed in the bin folder +ADD_SUBDIRECTORY(${SRCPROG} ${BIN}) +ADD_SUBDIRECTORY(${SRCANALYSIS} ${BIN}/Analysis) + +# Add a distclean target to the Makefile +ADD_CUSTOM_TARGET(distclean + COMMAND ${CMAKE_COMMAND} -P ${CMAKE_SOURCE_DIR}/distclean.cmake +) diff --git a/Libraries/Makefile b/Libraries/Makefile deleted file mode 100644 index 4f76ba2c8..000000000 --- a/Libraries/Makefile +++ /dev/null @@ -1,16 +0,0 @@ -FC=$(f90) -FLAGS=$(FL) - -all: - (cd Modules;make FC=$(FC) FLAGS="$(FLAGS)") ;\ - (cd MyEis;make FC=$(FC) FLAGS="$(FLAGS)") ;\ - (cd MyEis;make FC=$(FC) FLAGS="$(FLAGS)") ;\ - (cd MyNag;make FC=$(FC) FLAGS="$(FLAGS)") ;\ - (cd MyLin;make FC=$(FC) FLAGS="$(FLAGS)") - -clean: - (cd Modules;make clean);\ - (cd MyEis;make clean);\ - (cd MyNag;make clean);\ - (cd MyLin;make clean) - diff --git a/Libraries/Modules/Compile b/Libraries/Modules/Compile deleted file mode 100644 index b602a560e..000000000 --- a/Libraries/Modules/Compile +++ /dev/null @@ -1,13 +0,0 @@ -OBJS=mat_mod.o Random_Wrap.o errors.o Files_mod.o maxent.o matrix.o maxent_stoch.o fourier.o \ - Histogram.o lattices_v3.o Natural_constants.o log_mesh.o precdef.mod.o \ - Histogram_v2.o - -$(LIB): $(OBJS) - ar -r $(LIB) $(OBJS) - -.SUFFIXES: .f90 -.f90.o: - $(FC) $(SUFFIX) $(FLAGS) $< - -clean: - rm $(OBJS) diff --git a/Libraries/Modules/Makefile b/Libraries/Modules/Makefile deleted file mode 100644 index da09025d1..000000000 --- a/Libraries/Modules/Makefile +++ /dev/null @@ -1,15 +0,0 @@ -#FC= $(f90) -#FC= mpxlf90 -#FLAGS= -c -q64 -O4 -#FLAGS= -c -O3 -fbounds-check -FLAGS= -c -O3 -SUFFIX= -qsuffix=f=f90 -LF= -LIB=modules_90.a - -all: - (make -f Compile FC="$(FC)" FLAGS="$(FLAGS)" LIB="$(LIB)" ) - -clean: - (make -f Compile clean ) ;\ - rm *.mod *~ \#* diff --git a/Libraries/Modules/Makefile_Juropa b/Libraries/Modules/Makefile_Juropa deleted file mode 100644 index fb198fdf5..000000000 --- a/Libraries/Modules/Makefile_Juropa +++ /dev/null @@ -1,15 +0,0 @@ -FC= ifort -#FC= mpxlf90 -#FLAGS= -c -q64 -O4 -#FLAGS= -c -O1 -pg -FLAGS= -c -O3 -SUFFIX= -qsuffix=f=f90 -LF= -warn all -LIB=modules_90.a - -all: - (make -f Compile FC="$(FC)" FLAGS="$(FLAGS)" LIB="$(LIB)" ) - -clean: - (make -f Compile clean ) ;\ - rm *.mod diff --git a/Libraries/Modules/Makefile_cl b/Libraries/Modules/Makefile_cl deleted file mode 100644 index 624f8fd2e..000000000 --- a/Libraries/Modules/Makefile_cl +++ /dev/null @@ -1,14 +0,0 @@ -FC= ifort -#FC= mpxlf90 -#FLAGS= -c -q64 -O4 -FLAGS= -c -O3 -SUFFIX= -qsuffix=f=f90 -LF= -LIB=modules_90.a - -all: - (make -f Compile FC="$(FC)" FLAGS="$(FLAGS)" LIB="$(LIB)" ) - -clean: - (make -f Compile clean ) ;\ - rm *.mod *~ \#* diff --git a/Libraries/Modules/pre1 b/Libraries/Modules/pre1 deleted file mode 100644 index 38a649688..000000000 --- a/Libraries/Modules/pre1 +++ /dev/null @@ -1,12 +0,0 @@ -PRE = cpp -PREF = -P -OBJ= maxent_stoch.f90 - -all: $(OBJ) - -.SUFFIXES: .G90 .f90 -.G90.f90: - $(PRE) $(PREF) $? $@ - -clean: - rm $(OBJ) diff --git a/Prog_7/Compile_Hub b/Prog_7/Compile_Hub deleted file mode 100644 index ace754b8c..000000000 --- a/Prog_7/Compile_Hub +++ /dev/null @@ -1,16 +0,0 @@ -TARGET= Hubb.out -OBJS= control_mod.o Operator.o print_bin_mod.o Hamiltonian_Hub.o Hop_mod.o UDV_WRAP.o tau_m.o main.o wrapul.o cgr1.o wrapgrup.o wrapur.o upgrade.o \ - nranf.o wrapgrdo.o outconfc.o inconfc.o cgr2.o cgr2_2.o cgr2_1.o - - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Prog_7/Compile_SPT b/Prog_7/Compile_SPT deleted file mode 100644 index e7b63f777..000000000 --- a/Prog_7/Compile_SPT +++ /dev/null @@ -1,20 +0,0 @@ -TARGET= SPT.out -OBJS= control_mod.o Operator.o print_bin_mod.o Hamiltonian_SPT.o Hop_mod.o UDV_WRAP.o tau_m.o main.o wrapul.o cgr1.o wrapgrup.o wrapur.o upgrade.o \ - nranf.o wrapgrdo.o outconfc.o inconfc.o cgr2.o cgr2_2.o cgr2_1.o - -#block.o block_obs.o Mol_Dyn.o Hubb.o inconfc.o outconfc.o npbc.o salph.o sli.o \ -# sthop.o wrapul.o wrapur.o cgr1.o \ -# wrapgrdo.o mmthr.o mmthrm1.o mmthl.o mmthlm1.o obser.o upgradeu.o wrapgrup.o \ -# preq.o cgr2.o propr.o proprm1.o obsert.o prtau.o tau_m.o set_Hopping.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Prog_7/Ham_hop b/Prog_7/Ham_hop deleted file mode 100644 index e69de29bb..000000000 diff --git a/Prog_7/Ham_obser.f90 b/Prog_7/Ham_obser.f90 deleted file mode 100644 index 8fa394ab1..000000000 --- a/Prog_7/Ham_obser.f90 +++ /dev/null @@ -1,9 +0,0 @@ - Module - - Implicit none - - Complex (Kind=8) :: Phase - - - end SUBROUTINE Ham_obser - diff --git a/Prog_7/Hamiltonian_SPT.f90 b/Prog_7/Hamiltonian_SPT.f90 deleted file mode 100644 index e16b5adf5..000000000 --- a/Prog_7/Hamiltonian_SPT.f90 +++ /dev/null @@ -1,538 +0,0 @@ - !Model Hamiltonian for interaction-induced topological reduction - Module Hamiltonian - - Use Operator_mod - Use Lattices_v3 - Use MyMats - Use Random_Wrap - Use Files_mod - Use Matrix - - - Type (Operator), dimension(:,:), allocatable :: Op_V - Type (Operator), dimension(:,:), allocatable :: Op_T - Integer, allocatable :: nsigma(:,:) - Integer :: Ndim, N_FL, N_SUN, Ltrot - - - - ! What is below is private - - Type (Lattice), private :: Latt - Integer, parameter, private :: Norb=16 - Integer, allocatable, private :: List(:,:), Invlist(:,:) - Integer, private :: L1, L2 - real (Kind=8), private :: Ham_T, Ham_Vint, Ham_Lam - real (Kind=8), private :: Dtau, Beta - Character (len=64), private :: Model, Lattice_type - Complex (Kind=8), private :: Gamma_M(4,4,5), Sigma_M(2,2,0:3) - - - ! Observables - Integer, private :: Nobs - Complex (Kind=8), allocatable, private :: obs_scal(:) - Complex (Kind=8), allocatable, private :: Den_eq(:,:,:), Den_eq0(:) - - ! For time displaced - Integer, private :: NobsT - Complex (Kind=8), private :: Phase_tau - Complex (Kind=8), allocatable, private :: Green_tau(:,:,:,:), Den_tau(:,:,:,:) - - contains - - Subroutine Ham_Set - - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - integer :: ierr - - NAMELIST /VAR_lattice/ L1, L2, Lattice_type, Model - - NAMELIST /VAR_SPT/ ham_T, Ham_Vint, Ham_Lam, Dtau, Beta - - -#ifdef MPI - Integer :: Isize, Irank - Integer :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) - If (Irank == 0 ) then -#endif - OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) - IF (ierr /= 0) THEN - WRITE(*,*) 'unable to open ',ierr - STOP - END IF - READ(5,NML=VAR_lattice) - CLOSE(5) -#ifdef MPI - endif - - CALL MPI_BCAST(L1 ,1 ,MPI_INTEGER, 0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(L2 ,1 ,MPI_INTEGER, 0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Model ,64 ,MPI_CHARACTER, 0,MPI_COMM_WORLD,IERR) - CALL MPI_BCAST(Lattice_type,64 ,MPI_CHARACTER, 0,MPI_COMM_WORLD,IERR) -#endif - Call Ham_latt - - N_FL = 1 - N_SUN = 1 - -#ifdef MPI - If (Irank == 0 ) then -#endif - OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) - READ(5,NML=VAR_SPT) - CLOSE(5) -#ifdef MPI - endif - - CALL MPI_BCAST(ham_T ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(ham_Vint ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(ham_Lam ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Dtau ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Beta ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) -#endif - - Call Ham_hop - Ltrot = nint(beta/dtau) -#ifdef MPI - If (Irank == 0) then -#endif - Open (Unit = 50,file="info",status="unknown",position="append") - Write(50,*) '=====================================' - Write(50,*) 'Model is : ', Model - Write(50,*) 'Lattice is : ', Lattice_type - Write(50,*) 'Beta : ', Beta - Write(50,*) 'dtau,Ltrot : ', dtau,Ltrot - Write(50,*) 't : ', Ham_T - Write(50,*) 'V : ', Ham_Vint - Write(50,*) 'Lambda : ', Ham_Lam - close(50) -#ifdef MPI - endif -#endif - call Ham_V - end Subroutine Ham_Set -!============================================================================= - - Subroutine Ham_Latt - Implicit none - !Set the lattice - Integer :: no, I, nc - Real (Kind=8) :: a1_p(2), a2_p(2), L1_p(2), L2_p(2) - If ( Lattice_type =="Square" ) then - a1_p(1) = 1.0 ; a1_p(2) = 0.d0 - a2_p(1) = 0.0 ; a2_p(2) = 1.d0 - L1_p = dble(L1)*a1_p - L2_p = dble(L2)*a2_p - Call Make_Lattice( L1_p, L2_p, a1_p, a2_p, Latt ) - !Write(6,*) 'Lattice: ', Ndim - else - Write(6,*) "Lattice not yet implemented!" - Stop - endif - - Ndim = Latt%N*Norb - Allocate (List(Ndim,Norb), Invlist(Latt%N,Norb)) - nc = 0 - Do I = 1,Latt%N - Do no = 1,Norb - nc = nc + 1 - List(nc,1) = I - List(nc,2) = no - Invlist(I,no) = nc - ! no = 1..4 Xi_1 - ! no = 5..8 Xi_2 - ! no = 9..12 Xi_3 - ! no = 13..16 Xi_4 - Enddo - Enddo - - end Subroutine Ham_Latt - -!=================================================================================== - Subroutine Ham_hop - Implicit none - - ! Setup the hopping - ! Per flavor, the hopping is given by: - ! e^{-dtau H_t} = Prod_{n=1}^{Ncheck} e^{-dtau_n H_{n,t}} - - Integer :: I, I1, I2,I3,no, no1, n, Ncheck, nc , nth - Real (Kind=8) :: X - Complex (Kind=8) :: Z - - - ! Setup Gamma matrices - Gamma_M = cmplx(0.d0,0.d0) - Sigma_M = cmplx(0.d0,0.d0) - Sigma_M(1,1,0) = cmplx( 1.d0, 0.d0) - Sigma_M(2,2,0) = cmplx( 1.d0, 0.d0) - Sigma_M(1,2,1) = cmplx( 1.d0, 0.d0) - Sigma_M(2,1,1) = cmplx( 1.d0, 0.d0) - Sigma_M(1,2,2) = cmplx( 0.d0,-1.d0) - Sigma_M(2,1,2) = cmplx( 0.d0, 1.d0) - Sigma_M(1,1,3) = cmplx( 1.d0, 0.d0) - Sigma_M(2,2,3) = cmplx(-1.d0, 0.d0) - Do no = 1,2 - Do no1 = 1,2 - Gamma_M(no+2,no1 ,1) = Sigma_M(no,no1,0) - Gamma_M(no ,no1+2,1) = Sigma_M(no,no1,0) - Gamma_M(no+2,no1 ,2) = cmplx( 0.d0,-1.d0)*Sigma_M(no,no1,0) - Gamma_M(no ,no1+2,2) = cmplx( 0.d0, 1.d0)*Sigma_M(no,no1,0) - Gamma_M(no ,no1 ,3) = cmplx( 1.d0, 0.d0)*Sigma_M(no,no1,1) - Gamma_M(no+2,no1+2,3) = cmplx(-1.d0, 0.d0)*Sigma_M(no,no1,1) - Gamma_M(no ,no1 ,4) = cmplx( 1.d0, 0.d0)*Sigma_M(no,no1,2) - Gamma_M(no+2,no1+2,4) = cmplx(-1.d0, 0.d0)*Sigma_M(no,no1,2) - Gamma_M(no ,no1 ,5) = cmplx( 1.d0, 0.d0)*Sigma_M(no,no1,3) - Gamma_M(no+2,no1+2,5) = cmplx(-1.d0, 0.d0)*Sigma_M(no,no1,3) - Enddo - Enddo - - Ncheck = 1 - allocate(Op_T(Ncheck,N_FL)) - do n = 1,N_FL - Do nc = 1,NCheck - Call Op_make(Op_T(nc,n),Ndim) - DO I = 1, Ndim - Op_T(nc,n)%P(I) = I - enddo - Do I = 1,Latt%N - do nth = 0,3 - do no = 1,4 - do no1 = 1,4 - Z = cmplx(1.d0*Ham_T,0.d0)*Gamma_M(no,no1,3) - Op_T(nc,n)%O( invlist(I ,no + 4*nth), invlist(I ,no1 + 4*nth ) ) = Z - enddo - enddo - I1 = Latt%nnlist(I,1,0) - do no = 1,4 - do no1 = 1,4 - Z = (cmplx(0.d0,Ham_T)*Gamma_M(no,no1,1) + cmplx(Ham_T,0.d0)*Gamma_M(no,no1,3))/cmplx(2.d0,0.d0) - Op_T(nc,n)%O( invlist(I ,no + 4*nth), invlist(I1,no1 + 4*nth ) ) = Z - Op_T(nc,n)%O( invlist(I1,no1 + 4*nth), invlist(I ,no + 4*nth ) ) = conjg(Z) - enddo - enddo - I2 = Latt%nnlist(I,0,1) - do no = 1,4 - do no1 = 1,4 - Z = (cmplx(0.d0,Ham_Lam)*Gamma_M(no,no1,2) + cmplx(Ham_T,0.d0)*Gamma_M(no,no1,3))/cmplx(2.d0,0.d0) - Op_T(nc,n)%O( invlist(I ,no + 4*nth), invlist(I2,no1 + 4*nth ) ) = Z - Op_T(nc,n)%O( invlist(I2,no1 + 4*nth), invlist(I ,no + 4*nth ) ) = conjg(Z) - enddo - enddo - enddo - enddo - Op_T(nc,n)%g=cmplx(-Dtau,0.d0) - Call Op_set(Op_T(nc,n)) - ! Just for tests - Do I = 1, Ndim - Write(6,*) Op_T(nc,n)%E(i) - enddo - enddo - enddo - - - end Subroutine Ham_hop -!=================================================================================== - Subroutine Ham_V - - Implicit none - - Integer :: nf, nth, n, n1, n2, n3, n4, I, I1, I2, J, Ix, Iy, nc, no,no1, ns, npm - Integer :: nxy - Real (Kind=8) :: X_p(2), X1_p(2), X2_p(2), X, XJ, Xpm - - Complex (Kind=8) :: Ps(4,4,2), Ps_G5(4,4,2), Tmp(4,4), Z - Complex (Kind=8) :: Sx(16,16,2,2), Sy(16,16,2,2) - - - Ps = cmplx(0.d0,0.d0) - Call mmult (Tmp,Gamma_M(:,:,3), Gamma_M(:,:,4) ) - do ns = 1,2 - if (ns == 1) X = 1.d0/2d0 - if (ns == 2) X = -1.d0/2.d0 - Do I = 1,4 - Do J = 1,4 - Z = cmplx(0.d0,0.d0) - if ( I == J ) Z = cmplx(1.d0/2.d0,0.d0) - Ps(I,J,ns) = Z + cmplx(0.d0,X) * tmp(I,J) - Enddo - Enddo - Enddo - - Do ns = 1,2 - Call mmult ( Ps_G5(:,:,ns), Ps(:,:,ns), Gamma_M(:,:,5) ) - enddo - - Sx = cmplx(0.d0,0.d0) - Sy = cmplx(0.d0,0.d0) - Do ns = 1,2 - Do npm = 1,2 - if (npm == 1) Xpm = 1.0 - if (npm == 2) Xpm = -1.0 - Do no = 1,4 - do no1 = 1,4 - Sx(no , no1 + 4 ,ns,npm) = cmplx(1.d0, 0.d0)*Ps_G5(no,no1,ns) - Sx(no +4 , no1 ,ns,npm) = cmplx(1.d0, 0.d0)*Ps_G5(no,no1,ns) - Sx(no +8 , no1 + 12,ns,npm) = cmplx(xpm, 0.d0)*Ps_G5(no,no1,ns) - Sx(no+12 , no1 + 8 ,ns,npm) = cmplx(xpm, 0.d0)*Ps_G5(no,no1,ns) - - Sy(no , no1 + 4 ,ns,npm) = cmplx(0.d0, -1.d0 )*Ps_G5(no,no1,ns) - Sy(no +4 , no1 ,ns,npm) = cmplx(0.d0, 1.d0 )*Ps_G5(no,no1,ns) - Sy(no +8 , no1 + 12,ns,npm) = cmplx(0.d0, 1.d0*xpm)*Ps_G5(no,no1,ns) - Sy(no+12 , no1 + 8 ,ns,npm) = cmplx(0.d0, -1.d0*xpm)*Ps_G5(no,no1,ns) - enddo - enddo - enddo - enddo - - - ! Number of opertors 8 per unit cell - Allocate( Op_V(8*Latt%N,N_FL) ) - do nf = 1,N_FL - do i = 1, 8*Latt%N - Call Op_make(Op_V(i,nf),Norb) - enddo - enddo - nc = 0 - Do nf = 1,N_FL - do nxy = 1,2 - do ns = 1,2 - do npm = 1,2 - Xpm = 1.d0 - if (npm == 2) Xpm = -1.d0 - Do i = 1,Latt%N - nc = nc + 1 - Do no = 1,Norb - Op_V(nc,nf)%P(no) = Invlist(I,no) - enddo - Do no = 1,Norb - Do no1 = 1,Norb - If (nxy == 1) Op_V(nc,nf)%O(no,no1) = Sx(no,no1,ns,npm) - If (nxy == 2) Op_V(nc,nf)%O(no,no1) = Sy(no,no1,ns,npm) - Enddo - Enddo - Op_V(nc,nf)%g = SQRT(CMPLX(-Xpm*DTAU*Ham_Vint/8.d0,0.D0)) - Op_V(nc,nf)%alpha = cmplx(0.d0,0.d0) - Op_V(nc,nf)%type = 2 - Call Op_set( Op_V(nc,nf) ) - ! The operator reads: - ! g*s*( c^{dagger} O c - alpha )) - ! with s the HS field. - Enddo - Enddo - Enddo - Enddo - Enddo - - end Subroutine Ham_V - -!=================================================================================== - Real (Kind=8) function S0(n,nt) - Implicit none - Integer, Intent(IN) :: n,nt - Integer :: i, nt1 - S0 = 1.d0 - end function S0 - -!=================================================================================== - Subroutine Alloc_obs(Ltau) - - Implicit none - Integer, Intent(In) :: Ltau - Integer :: I - Allocate ( Obs_scal(5) ) - Allocate ( Den_eq(Latt%N,Norb,Norb), Den_eq0(Norb) ) - If (Ltau == 1) then - Allocate ( Green_tau(Latt%N,Ltrot+1,Norb,Norb), Den_tau(Latt%N,Ltrot+1,Norb,Norb) ) - endif - - end Subroutine Alloc_obs - -!=================================================================================== - - Subroutine Init_obs(Ltau) - - Implicit none - Integer, Intent(In) :: Ltau - - Integer :: I,n - - Nobs = 0 - Obs_scal = cmplx(0.d0,0.d0) - Den_eq = cmplx(0.d0,0.d0) - Den_eq0 = cmplx(0.d0,0.d0) - - If (Ltau == 1) then - NobsT = 0 - Phase_tau = cmplx(0.d0,0.d0) - Green_tau = cmplx(0.d0,0.d0) - Den_tau = cmplx(0.d0,0.d0) - endif - - end Subroutine Init_obs - -!======================================================================== - Subroutine Obser(GR,Phase,Ntau) - - Implicit none - - Complex (Kind=8), INTENT(IN) :: GR(Ndim,Ndim,N_FL) - Complex (Kind=8), Intent(IN) :: PHASE - Integer, INTENT(IN) :: Ntau - - !Local - Complex (Kind=8) :: GRC(Ndim,Ndim,N_FL), ZK - Complex (Kind=8) :: Zrho, Zkin, ZPot, Z, ZP,ZS - Integer :: I,J, no,no1, n, n1, imj, nf, I1, I2, J1, J2 - - Real (Kind=8) :: G(4,4), X, FI, FJ - - Nobs = Nobs + 1 - ZP = PHASE/cmplx(Real(Phase,kind=8),0.d0) - ZS = cmplx(Real(Phase,kind=8)/Abs(Real(Phase,kind=8)), 0.d0) - - - Do nf = 1,N_FL - Do I = 1,Ndim - Do J = 1,Ndim - ZK = cmplx(0.d0,0.d0) - If ( I == J ) ZK = cmplx(1.d0,0.d0) - GRC(I,J,nf) = ZK - GR(J,I,nf) - Enddo - Enddo - Enddo - ! GRC(i,j,nf) = < c^{dagger}_{j,nf } c_{j,nf } > - ! Compute scalar observables. - - Zkin = cmplx(0.d0,0.d0) - Do nf = 1,N_FL - Do J = 1,Op_T(1,nf)%N - J1 = Op_T(1,nf)%P(J) - DO I = 1,Op_T(1,nf)%N - I1 = Op_T(1,nf)%P(I) - Zkin = Zkin + Op_T(1,nf)%O(i,j)*Grc(i1,j1,nf) - Enddo - ENddo - Enddo - Zkin = Zkin*cmplx( dble(N_SUN), 0.d0 ) - - Zrho = cmplx(0.d0,0.d0) - Do nf = 1,N_FL - Do I = 1,Ndim - Zrho = Zrho + Grc(i,i,nf) - enddo - enddo - Zrho = Zrho*cmplx( dble(N_SUN), 0.d0 ) - - ZPot = cmplx(0.d0,0.d0) - - Obs_scal(1) = Obs_scal(1) + zrho * ZP*ZS - Obs_scal(2) = Obs_scal(2) + zkin * ZP*ZS - Obs_scal(3) = Obs_scal(3) + Zpot * ZP*ZS - Obs_scal(4) = Obs_scal(4) + (zkin + Zpot)*ZP*ZS - Obs_scal(5) = Obs_scal(5) + ZS - ! You will have to allocate more space if you want to include more scalar observables. - DO I1 = 1,Ndim - I = List(I1,1) - no = List(I1,2) - DO J1 = 1, Ndim - J = List(J1,1) - no1 = list(J1,2) - imj = latt%imj(I,J) - - DEN_Eq (imj,no,no1) = DEN_Eq (imj,no,no1) + & - & ( GRC(I1,J1,1) * GR (I1,J1,1) + & - & GRC(I1,I1,1) * GRC(J1,J1,1) ) * ZP*ZS - - enddo - Den_eq0(no) = Den_eq0(no) + GRC(I1,I1,1)*ZP*ZS - enddo - - end Subroutine Obser -!========================================================== - - Subroutine Pr_obs(LTAU) - - Use Print_bin_mod - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - - Integer, Intent(In) :: Ltau - - Character (len=64) :: File_pr - Complex (Kind=8) :: Phase_bin -#ifdef MPI - Integer :: Isize, Irank, Ierr - Integer :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'In Pr_obs', LTAU -!!$#else -!!$ Write(6,*) 'In Pr_obs', LTAU -!!$#endif - - Phase_bin = Obs_scal(5)/cmplx(dble(Nobs),0.d0) - File_pr ="Den_eq" - Call Print_bin(Den_eq, Den_eq0, Latt, Nobs, Phase_bin, file_pr) - - File_pr ="ener" - Call Print_scal(Obs_scal, Nobs, file_pr) - If (Ltau == 1) then - Phase_tau = Phase_tau/cmplx(dble(NobsT),0.d0) - File_pr = "Green_tau" - Call Print_bin_tau(Green_tau,Latt,NobsT,Phase_tau, file_pr,dtau) - File_pr = "Den_tau" - Call Print_bin_tau(Den_tau,Latt,NobsT,Phase_tau, file_pr,dtau) - endif -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'out Pr_obs', LTAU -!!$#else -!!$ Write(6,*) 'out Pr_obs', LTAU -!!$#endif - end Subroutine Pr_obs -!========================================================== - - Subroutine OBSERT(NT, GT0,G0T,G00,GTT, PHASE) - Implicit none - - Integer , INTENT(IN) :: NT - Complex (Kind=8), INTENT(IN) :: GT0(Ndim,Ndim,N_FL),G0T(Ndim,Ndim,N_FL),G00(Ndim,Ndim,N_FL),GTT(Ndim,Ndim,N_FL) - Complex (Kind=8), INTENT(IN) :: Phase - - !Locals - Complex (Kind=8) :: Z, ZP, ZS - Integer :: IMJ, I, J - - ZP = PHASE/cmplx(Real(Phase,kind=8),0.d0) - ZS = cmplx(Real(Phase,kind=8)/Abs(Real(Phase,kind=8)), 0.d0) - If (NT == 0 ) then - Phase_tau = Phase_tau + ZS - NobsT = NobsT + 1 - endif - If ( N_FL == 1 ) then - Z = cmplx(dble(N_SUN),0.d0) - Do I = 1,Latt%N - Do J = 1,Latt%N - imj = latt%imj(I,J) - Green_tau(imj,nt+1,1,1) = green_tau(imj,nt+1,1,1) + Z * GT0(I,J,1) * ZP* ZS - Den_tau (imj,nt+1,1,1) = Den_tau (imj,nt+1,1,1) - Z * GT0(I,J,1)*G0T(J,I,1) * ZP* ZS - Enddo - Enddo - Endif - end Subroutine OBSERT - - - end Module Hamiltonian diff --git a/Prog_7/Makefile b/Prog_7/Makefile deleted file mode 100644 index 5db5d6c9a..000000000 --- a/Prog_7/Makefile +++ /dev/null @@ -1,22 +0,0 @@ -FC= $(mpif90) -FC= $(f90) -FLAGS= $(FL) -LF = $(Lflags) -LIBS= $(Libs)/Modules/modules_90.a \ - $(Libs)/MyEis/libeis.a \ - $(Libs)/MyNag/libnag.a \ - $(Libs)/MyLin/liblin.a \ - $(LIB_BLAS_LAPACK) - -Hub: - cp $(Libs)/Modules/*.mod . ;\ - (make -f Compile_Hub FC="$(FC)" FLAGS="$(FLAGS)" LIBS="$(LIBS)" ) - -SPT: - cp $(Libs)/Modules/*.mod . ;\ - (make -f Compile_SPT FC="$(FC)" FLAGS="$(FLAGS)" LIBS="$(LIBS)" LF="$(LF)" ) - -clean: - (make -f Compile_Hub clean );\ - (make -f Compile_SPT clean );\ - rm *.mod *~ \#* diff --git a/Prog_7/Operator.f90 b/Prog_7/Operator.f90 deleted file mode 100644 index 35256dcd1..000000000 --- a/Prog_7/Operator.f90 +++ /dev/null @@ -1,420 +0,0 @@ -Module Operator_mod - - Use MyMats - - Implicit none - - Real (Kind=8) :: Phi(-2:2,2), Gaml(-2:2,2) - Integer :: NFLIPL(-2:2,3) - - - ! What information should the operator contain - Type Operator - Integer :: N - complex (kind=8), pointer :: O(:,:), U (:,:) - Real (kind=8), pointer :: E(:) - Integer, pointer :: P(:) - complex (kind=8) :: g - complex (kind=8) :: alpha - Integer :: Type - ! P is an N X Ndim matrix such that P.T*O*P* = A - ! P has only one non-zero entry per column which is specified by P - ! All in all. g * Phi(s,type) * ( c^{dagger} A c + alpha ) - ! The variable Type allows you to define the type of HS. - end type Operator - - -Contains - - Subroutine Op_SetHS - Implicit none - Integer :: n - Phi = 0.d0 - do n = -2,2 - Phi(n,1) = real(n,kind=8) - enddo - Phi(-2,2) = - SQRT(2.D0 * ( 3.D0 + SQRT(6.D0) ) ) - Phi(-1,2) = - SQRT(2.D0 * ( 3.D0 - SQRT(6.D0) ) ) - Phi( 1,2) = SQRT(2.D0 * ( 3.D0 - SQRT(6.D0) ) ) - Phi( 2,2) = SQRT(2.D0 * ( 3.D0 + SQRT(6.D0) ) ) - - Do n = -2,2 - gaml(n,1) = 1.d0 - Enddo - GAML(-2,2) = 1.D0 - SQRT(6.D0)/3.D0 - GAML( 2,2) = 1.D0 - SQRT(6.D0)/3.D0 - GAML(-1,2) = 1.D0 + SQRT(6.D0)/3.D0 - GAML( 1,2) = 1.D0 + SQRT(6.D0)/3.D0 - - NFLIPL(-2,1) = -1 - NFLIPL(-2,2) = 1 - NFLIPL(-2,3) = 2 - - NFLIPL(-1,1) = 1 - NFLIPL(-1,2) = 2 - NFLIPL(-1,3) = -2 - - NFLIPL( 1,1) = 2 - NFLIPL( 1,2) = -2 - NFLIPL( 1,3) = -1 - - NFLIPL( 2,1) = -2 - NFLIPL( 2,2) = -1 - NFLIPL( 2,3) = 1 - - end Subroutine Op_SetHS - - Subroutine Op_phase(Phase,OP_V,Nsigma,N_SUN) ! This also goes in Operator (Input is nsigma, Op_V). - Implicit none - - Complex (Kind=8), Intent(Inout) :: Phase - Integer, Intent(IN) :: N_SUN - Integer, dimension(:,:), Intent(In) :: Nsigma - Type (Operator), dimension(:,:), Intent(In) :: Op_V - - Integer :: n, nf, nt - - do nf = 1,Size(Op_V,2) - do n = 1,size(Op_V,1) - do nt = 1,size(nsigma,2) - Phase = Phase*exp( Op_V(n,nf)%g * Op_V(n,nf)%alpha * Phi(nsigma(n,nt),Op_V(n,nf)%type) ) - enddo - enddo - enddo - Phase = Phase**dble(N_SUN) - - end Subroutine Op_phase - - - subroutine Op_make(Op,N) - Implicit none - Type (Operator), intent(INOUT) :: Op - Integer, Intent(IN) :: N - Allocate (Op%O(N,N), Op%U(N,N), Op%E(N), Op%P(N)) - Op%O = cmplx(0.d0,0.d0) - Op%U = cmplx(0.d0,0.d0) - Op%E = 0.d0 - Op%P = 0 - Op%N = N - Op%g = cmplx(0.d0,0.d0) - end subroutine Op_make - - subroutine Op_clear(Op,N) - Implicit none - Type (Operator), intent(INOUT) :: Op - Integer, Intent(IN) :: N - Deallocate (Op%O, Op%U, Op%E, Op%P) - end subroutine Op_clear - - subroutine Op_set(Op) - Implicit none - Type (Operator), intent(INOUT) :: Op - If (Op%N > 1) then - !Write(6,*) 'Calling diag', Op%O(1,2), Size(Op%O,1), Size(Op%U,1), Size(Op%E,1) - Call Diag(Op%O,Op%U,Op%E) - !Write(6,*) 'Calling diag 1' - else - Op%E(1) = Op%O(1,1) - Op%U(1,1) = cmplx(1.d0,0.d0) - endif - end subroutine Op_set - - - subroutine Op_exp(g,Op,Mat) - Implicit none - Type (Operator), Intent(IN) :: Op - Complex (Kind=8), Dimension(:,:), INTENT(OUT) :: Mat - Complex (Kind=8), INTENT(IN) :: g - Complex (Kind=8) :: Z, Z1 - - Integer :: n, i,j - - Mat = cmplx(0.d0,0.d0) - Do n = 1,Op%N - Z = exp(g*cmplx(Op%E(n),0.d0)) - do J = 1,Op%N - Z1 = Z*conjg(Op%U(J,n)) - Do I = 1,Op%N - Mat(I,J) = Mat(I,J) + Op%U(I,n)*Z1 - enddo - enddo - enddo - end subroutine Op_exp - - subroutine Op_mmultL(Mat,Op,spin,Ndim) - Implicit none - Integer :: Ndim - Type (Operator) , INTENT(IN ) :: Op - Complex (Kind=8), INTENT(INOUT) :: Mat (Ndim,Ndim) - Real (Kind=8), INTENT(IN ) :: spin - - - ! Local - Complex (Kind=8) :: VH(Ndim,Op%N), Z, Z1 - Integer :: n, i, m, m1 - - - ! In Mat - ! Out Mat = Mat*exp(spin*Op) - - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Z = exp(Op%g*cmplx(Op%E(n)*spin,0.d0)) - Do m = 1,Op%N - Z1 = Op%U(m,n)* Z - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z1 = conjg(Op%U(n,m)) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - - - end subroutine Op_mmultL - - subroutine Op_mmultR(Mat,Op,spin,Ndim) - Implicit none - Integer :: Ndim - Type (Operator) , INTENT(IN ) :: Op - Complex (Kind=8), INTENT(INOUT) :: Mat (Ndim,Ndim) - Real (Kind=8), INTENT(IN ) :: spin - - - ! Local - Complex (Kind=8) :: VH(Ndim,Op%N), Z, Z1 - Integer :: n, i, m, m1 - - ! In Mat - ! Out Mat = exp(spin*Op)*Mat - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Z1 = exp(Op%g*cmplx(Op%E(n)*spin,0.d0)) - Do m = 1,Op%N - Z = conjg(Op%U(m,n))* Z1 - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z* Mat(Op%P(m),I) - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z = Op%U(n,m) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z* Mat(Op%P(m),I) - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - - - end subroutine Op_mmultR - - Subroutine Op_Wrapup(Mat,Op,spin,Ndim,N_Type) - - Implicit none - - Integer :: Ndim - Type (Operator) , INTENT(IN ) :: Op - Complex (Kind=8), INTENT(INOUT) :: Mat (Ndim,Ndim) - Real (Kind=8), INTENT(IN ) :: spin - Integer, INTENT(IN) :: N_Type - - ! Local - Complex (Kind=8) :: VH(Ndim,Op%N), Z, Z1 - Integer :: n, i, m, m1 - - - - - !!!!! N_Type ==1 - ! exp(Op%g*spin*Op%E)*(Op%U^{dagger})*Mat*Op%U*exp(-Op%g*spin*Op%E) - ! - !!!!! - !!!!! N_Type == 2 - ! Op%U * Mat * (Op%U^{dagger}) - !!!!! - If (N_type == 1) then - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Z = exp(-Op%g*cmplx(Op%E(n)*spin,0.d0)) - Do m = 1,Op%N - Z1 = Op%U(m,n) * Z - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Z = exp(Op%g*cmplx(Op%E(n)*spin,0.d0)) - Do m = 1,Op%N - Z1 = Z * conjg(Op%U(m,n)) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z1* Mat(Op%P(m),I) - Enddo - enddo - enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - elseif (N_Type == 2) then - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z1 = conjg(Op%U(n,m)) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z1 = Op%U(n,m) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z1* Mat(Op%P(m),I) - Enddo - enddo - enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - endif - end Subroutine Op_Wrapup - - Subroutine Op_Wrapdo(Mat,Op,spin,Ndim,N_Type) - - Implicit none - - Integer :: Ndim - Type (Operator) , INTENT(IN ) :: Op - Complex (Kind=8), INTENT(INOUT) :: Mat (Ndim,Ndim) - Real (Kind=8), INTENT(IN ) :: spin - Integer, INTENT(IN) :: N_Type - - ! Local - Complex (Kind=8) :: VH(Ndim,Op%N), Z, Z1 - Integer :: n, i, m, m1 - - !!!!! N_Type == 1 - ! Op%U*exp(-Op%g*spin*Op%E)*Mat*exp(Op%g*spin*Op%E)*(Op%U^{dagger}) - ! - !!!!! - !!!!! N_Type == 2 - ! (Op%U^{dagger}) * Mat * Op%U - !!!!! - If (N_type == 1) then - VH = cmplx(0.d0,0.d0) - Do m = 1,Op%N - Z = exp(Op%g*cmplx(Op%E(m)*spin,0.d0)) - do n = 1,Op%N - Z1 = Z * conjg(Op%U(n,m)) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - VH = cmplx(0.d0,0.d0) - Do m = 1,Op%N - Z = exp(-Op%g*cmplx(Op%E(m)*spin,0.d0)) - do n = 1,Op%N - Z1 = Z * Op%U(n,m) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z1* Mat(Op%P(m),I) - Enddo - enddo - enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - elseif (N_Type == 2) then - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z1 = Op%U(m,n) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 - Enddo - enddo - Enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(I,Op%P(n)) = VH(I,n) - Enddo - Enddo - - VH = cmplx(0.d0,0.d0) - do n = 1,Op%N - Do m = 1,Op%N - Z1 = conjg(Op%U(m,n)) - DO I = 1,Ndim - VH(I,n) = VH(I,n) + Z1* Mat(Op%P(m),I) - Enddo - enddo - enddo - Do n = 1,Op%N - Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) - Enddo - Enddo - endif - - end Subroutine Op_Wrapdo - - -end Module Operator_mod diff --git a/Prog_7/print_bin_mod.f90 b/Prog_7/print_bin_mod.f90 deleted file mode 100644 index 9f5297d0e..000000000 --- a/Prog_7/print_bin_mod.f90 +++ /dev/null @@ -1,296 +0,0 @@ - Module Print_bin_mod - - Interface Print_bin - module procedure Print_bin_C, Print_bin_R - end Interface Print_bin - - Interface Print_bin_tau - module procedure Print_bin_tau_C - end Interface Print_bin_tau - - Contains - - Subroutine Print_bin_C(Dat_eq,Dat_eq0,Latt, Nobs, Phase_bin, file_pr) - Use Lattices_v3 - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - Complex (Kind=8), Dimension(:,:,:), Intent(inout):: Dat_eq - Complex (Kind=8), Dimension(:) , Intent(inout):: Dat_eq0 - Type (Lattice), Intent(In) :: Latt - Complex (Kind=8), Intent(Inout):: Phase_bin - Character (len=64), Intent(In) :: File_pr - Integer, Intent(In) :: Nobs - - ! Local - Integer :: Norb, I, no,no1 - Complex (Kind=8), allocatable :: Tmp(:,:,:), Tmp1(:) - Real (Kind=8) :: x_p(2) -#ifdef MPI - Complex (Kind=8):: Z - Integer :: Ierr, Isize, Irank - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - Norb = size(Dat_eq,3) - if ( .not. (Latt%N == Size(Dat_eq,1) ) ) then - Write(6,*) 'Error in Print_bin' - Stop - endif - Allocate (Tmp(Latt%N,Norb,Norb), Tmp1(Norb) ) - Dat_eq = Dat_eq/cmplx(dble(Nobs),0.d0) - Dat_eq0 = Dat_eq0/cmplx(dble(Nobs)*dble(Latt%N),0.d0) - -#ifdef MPI - I = Latt%N*Norb*Norb - Tmp = cmplx(0.d0,0.d0) - CALL MPI_REDUCE(Dat_eq,Tmp,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Dat_eq = Tmp/CMPLX(DBLE(ISIZE),0.D0) - I = 1 - CALL MPI_REDUCE(Phase_bin,Z,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Phase_bin= Z/CMPLX(DBLE(ISIZE),0.D0) - - I = Norb - CALL MPI_REDUCE(Dat_eq0,Tmp1,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Dat_eq0 = Tmp1/CMPLX(DBLE(ISIZE),0.D0) - - If (Irank == 0 ) then -#endif - do no = 1,Norb - do no1 = 1,Norb - Call Fourier_R_to_K(Dat_eq(:,no,no1), Tmp(:,no,no1), Latt) - enddo - enddo - Open (Unit=10,File=File_pr, status="unknown", position="append") - Write(10,*) dble(Phase_bin),Norb,Latt%N - do no = 1,Norb - Write(10,*) Dat_eq0(no) - enddo - do I = 1,Latt%N - x_p = dble(Latt%listk(i,1))*Latt%b1_p + dble(Latt%listk(i,2))*Latt%b2_p - Write(10,*) X_p(1), X_p(2) - do no = 1,Norb - do no1 = 1,Norb - Write(10,*) tmp(I,no,no1) - enddo - enddo - enddo - close(10) -#ifdef MPI - Endif -#endif - - deallocate (Tmp, tmp1 ) - - - End Subroutine Print_bin_C - - -!========================================================= - - Subroutine Print_bin_R(Dat_eq,Dat_eq0,Latt, Nobs, Phase_bin, file_pr) - Use Lattices_v3 - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - Real (Kind=8), Dimension(:,:,:), Intent(inout) :: Dat_eq - Real (Kind=8), Dimension(:) , Intent(inout) :: Dat_eq0 - Type (Lattice), Intent(In) :: Latt - Complex (Kind=8), Intent(Inout) :: Phase_bin - Character (len=64), Intent(In) :: File_pr - Integer, Intent(In) :: Nobs - - ! Local - Integer :: Norb, I, no,no1 - Real (Kind=8), allocatable :: Tmp(:,:,:), Tmp1(:) - Real (Kind=8) :: x_p(2) -#ifdef MPI - Integer :: Ierr, Isize, Irank - Complex (Kind=8) :: Z - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - Norb = size(Dat_eq,3) - if ( .not. (Latt%N == Size(Dat_eq,1) ) ) then - Write(6,*) 'Error in Print_bin' - Stop - endif - Allocate (Tmp(Latt%N,Norb,Norb), Tmp1(Norb) ) - Dat_eq = Dat_eq/dble(Nobs) - Dat_eq0 = Dat_eq0/(dble(Nobs)*dble(Latt%N)) -#ifdef MPI - I = Latt%N*Norb*Norb - Tmp = 0.d0 - CALL MPI_REDUCE(Dat_eq,Tmp,I,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Dat_eq = Tmp/DBLE(ISIZE) - I = 1 - CALL MPI_REDUCE(Phase_bin,Z,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Phase_bin= Z/CMPLX(DBLE(ISIZE),0.D0) - If (Irank == 0 ) then - - I = Norb - CALL MPI_REDUCE(Dat_eq0,Tmp1,I,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Dat_eq0 = Tmp1/CMPLX(DBLE(ISIZE),0.D0) - -#endif - - do no = 1,Norb - do no1 = 1,Norb - Call Fourier_R_to_K(Dat_eq(:,no,no1), Tmp(:,no,no1), Latt) - enddo - enddo - Open (Unit=10,File=File_pr, status="unknown", position="append") - Write(10,*) dble(Phase_bin),Norb,Latt%N - do no = 1,Norb - Write(10,*) Dat_eq0(no) - enddo - do I = 1,Latt%N - x_p = dble(Latt%listk(i,1))*Latt%b1_p + dble(Latt%listk(i,2))*Latt%b2_p - Write(10,*) X_p(1), X_p(2) - do no = 1,Norb - do no1 = 1,Norb - Write(10,*) tmp(I,no,no1) - enddo - enddo - enddo - close(10) -#ifdef MPI - endif -#endif - deallocate (Tmp ) - - End Subroutine Print_bin_R -!============================================================ - Subroutine Print_scal(Obs, Nobs, file_pr) - - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - Complex (Kind=8), Dimension(:), Intent(inout) :: Obs - Character (len=64), Intent(In) :: File_pr - Integer, Intent(In) :: Nobs - - ! Local - Integer :: Norb,I - Complex (Kind=8), allocatable :: Tmp(:) -#ifdef MPI - Integer :: Ierr, Isize, Irank - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - Norb = size(Obs,1) - Allocate ( Tmp(Norb) ) - Obs = Obs/cmplx(dble(Nobs),0.d0) -#ifdef MPI - Tmp = 0.d0 - CALL MPI_REDUCE(Obs,Tmp,Norb,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Obs = Tmp/cmplx(DBLE(ISIZE),0.d0) - if (Irank == 0 ) then -#endif - Open (Unit=10,File=File_pr, status="unknown", position="append") - WRITE(10,*) (Obs(I), I=1,size(Obs,1)) - close(10) -#ifdef MPI - endif -#endif - deallocate (Tmp ) - - End Subroutine Print_scal - -!============================================================== - Subroutine Print_bin_tau_C(Dat_tau,Latt, Nobs, Phase_bin, file_pr, dtau) - Use Lattices_v3 - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - Complex (Kind=8), Dimension(:,:,:,:), Intent(inout):: Dat_tau - Type (Lattice), Intent(In) :: Latt - Complex (Kind=8), Intent(In) :: Phase_bin - Character (len=64), Intent(In) :: File_pr - Integer, Intent(In) :: Nobs - Real (kind=8), Intent(In) :: dtau - - ! Local - Integer :: Norb, I, no,no1, LT, nt - Complex (Kind=8), allocatable :: Tmp(:,:,:,:) - Complex (Kind=8) :: Phase_mean - Real (Kind=8) :: x_p(2) -#ifdef MPI - Complex (Kind=8):: Z - Integer :: Ierr, Isize, Irank - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - Phase_mean = Phase_bin - Norb = size(Dat_tau,3) - if ( .not. (Latt%N == Size(Dat_tau,1) ) ) then - Write(6,*) 'Error in Print_bin' - Stop - endif - LT = Size(Dat_tau,2) - Allocate (Tmp(Latt%N,LT,Norb,Norb) ) - Dat_tau = Dat_tau/cmplx(dble(Nobs),0.d0) - -#ifdef MPI - I = Latt%N*Norb*Norb*LT - Tmp = cmplx(0.d0,0.d0) - CALL MPI_REDUCE(Dat_tau,Tmp,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Dat_tau = Tmp/CMPLX(DBLE(ISIZE),0.D0) - I = 1 - CALL MPI_REDUCE(Phase_mean,Z,I,MPI_COMPLEX16,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Phase_mean= Z/CMPLX(DBLE(ISIZE),0.D0) - If (Irank == 0 ) then -#endif - do nt = 1,LT - do no = 1,Norb - do no1 = 1,Norb - Call Fourier_R_to_K(Dat_tau(:,nt,no,no1), Tmp(:,nt,no,no1), Latt) - enddo - enddo - enddo - Open (Unit=10,File=File_pr, status="unknown", position="append") - Write(10,*) dble(Phase_mean),Norb,Latt%N, LT, dtau - do I = 1,Latt%N - x_p = dble(Latt%listk(i,1))*Latt%b1_p + dble(Latt%listk(i,2))*Latt%b2_p - Write(10,*) X_p(1), X_p(2) - Do nt = 1,LT - do no = 1,Norb - do no1 = 1,Norb - Write(10,*) tmp(I,nt,no,no1) - enddo - enddo - enddo - enddo - close(10) -#ifdef MPI - Endif -#endif - - deallocate (Tmp ) - - - End Subroutine Print_bin_tau_C - - - - end Module Print_bin_mod diff --git a/Prog_7/upgrade.f90 b/Prog_7/upgrade.f90 deleted file mode 100644 index 0dd7cac22..000000000 --- a/Prog_7/upgrade.f90 +++ /dev/null @@ -1,151 +0,0 @@ - Subroutine Upgrade(GR,N_op,NT,PHASE,Op_dim) - - Use Hamiltonian - Use Random_wrap - Use Control - Use Precdef - Implicit none - - Complex (Kind=double) :: GR(Ndim,Ndim, N_FL) - Integer, INTENT(IN) :: N_op, Nt, Op_dim - Complex (Kind=double) :: Phase - - ! Local :: - Complex (Kind=double) :: Mat(Op_dim,Op_Dim), Delta(Op_dim,N_FL) - Complex (Kind=double) :: Ratio(N_FL), Ratiotot, Z1 - Integer :: ns_new, ns_old, n,m,nf, i,j - Complex (Kind= double) :: ZK, Z, D_Mat - Integer, external :: nranf - - Real (Kind = double) :: Weight - Complex (Kind = double) :: u(Ndim,Op_dim), v(Ndim,Op_dim) - Complex (Kind = double) :: x_v(Ndim,Op_dim), y_v(Ndim,Op_dim), xp_v(Ndim,Op_dim) - Complex (Kind = double) :: s_xv(Op_dim), s_yu(Op_dim) - - Logical :: Log - - - if ( sqrt(dble(OP_V(n_op,1)%g*conjg(OP_V(n_op,1)%g))) < 1.D-6 ) return - - ! Compute the ratio - nf = 1 - ns_old = nsigma(n_op,nt) - If ( Op_V(n_op,nf)%type == 1) then - ns_new = -ns_old - else - ns_new = NFLIPL(Ns_old,nranf(3)) - endif - Do nf = 1,N_FL - Z1 = Op_V(n_op,nf)%g * cmplx( Phi(ns_new,Op_V(n_op,nf)%type) - Phi(ns_old,Op_V(n_op,nf)%type), 0.d0) - Do m = 1,Op_V(n_op,nf)%N - Z = exp( Z1* Op_V(n_op,nf)%E(m) ) - cmplx(1.d0,0.d0) - Delta(m,nf) = Z - do n = 1,Op_V(n_op,nf)%N - ZK = cmplx(0.d0,0.d0) - If (n == m ) ZK = cmplx(1.d0,0.d0) - Mat(n , m ) = ZK + ( ZK - GR( Op_V(n_op,nf)%P(n), Op_V(n_op,nf)%P(m),nf )) * Z - Enddo - Enddo - If (Size(Mat,1) == 1 ) then - D_mat = Mat(1,1) - elseif (Size(Mat,1) == 2 ) then - D_mat = Mat(1,1)*Mat(2,2) - Mat(2,1)*Mat(1,2) - else - D_mat = Det(Mat,Size(Mat,1)) - !Write(6,*) 'Not yet programed! ' - !Stop - endif - Ratio(nf) = D_Mat * exp( Z1*Op_V(n_op,nf)%alpha ) - Enddo - - Ratiotot = cmplx(1.d0,0.d0) - Do nf = 1,N_FL - Ratiotot = Ratiotot * Ratio(nf) - enddo - nf = 1 - Ratiotot = (Ratiotot**dble(N_SUN)) * cmplx(Gaml(ns_new, Op_V(n_op,nf)%type)/Gaml(ns_old, Op_V(n_op,nf)%type),0.d0) - Ratiotot = Ratiotot*cmplx(S0(n_op,nt),0.d0) - - - !Write(6,*) Ratiotot - - Weight = abs( real(Phase * Ratiotot, kind=double)/real(Phase,kind=double) ) - - Log = .false. - if ( Weight > ranf() ) Then - Log = .true. - Phase = Phase * Ratiotot/cmplx(weight,0.d0) - !Write(6,*) 'Accepted : ', Ratiotot - - Do nf = 1,N_FL - ! Setup u(i,n), v(n,i) - u = cmplx(0.d0,0.d0) - v = cmplx(0.d0,0.d0) - do n = 1,Op_V(n_op,nf)%N - u( Op_V(n_op,nf)%P(n), n) = Delta(n,nf) - do i = 1,Ndim - v(i,n) = - GR( Op_V(n_op,nf)%P(n), i, nf ) - enddo - v(Op_V(n_op,nf)%P(n), n) = cmplx(1.d0,0.d0) - GR( Op_V(n_op,nf)%P(n), Op_V(n_op,nf)%P(n), nf) - enddo - - - x_v = cmplx(0.d0,0.d0) - y_v = cmplx(0.d0,0.d0) - i = Op_V(n_op,nf)%P(1) - x_v(i,1) = u(i,1)/(cmplx(1.d0,0.d0) + v(i,1)*u(i,1) ) - Do i = 1,Ndim - y_v(i,1) = v(i,1) - enddo - do n = 2,Op_V(n_op,nf)%N - s_yu = cmplx(0.d0,0.d0) - s_xv = cmplx(0.d0,0.d0) - do m = 1,n-1 - do i = 1,Ndim - s_yu(m) = s_yu(m) + y_v(i,m)*u(i,n) - s_xv(m) = s_xv(m) + x_v(i,m)*v(i,n) - enddo - enddo - Do i = 1,Ndim - x_v(i,n) = u(i,n) - y_v(i,n) = v(i,n) - enddo - Z = cmplx(1.d0,0.d0) + u( Op_V(n_op,nf)%P(n), n)*v(Op_V(n_op,nf)%P(n),n) - do m = 1,n-1 - Z = Z - s_xv(m)*s_yu(m) - Do i = 1,Ndim - x_v(i,n) = x_v(i,n) - x_v(i,m)*s_yu(m) - y_v(i,n) = y_v(i,n) - y_v(i,m)*s_xv(m) - enddo - enddo - Do i = 1,Ndim - x_v(i,n) = x_v(i,n)/Z - Enddo - enddo - xp_v = cmplx(0.d0,0.d0) - do n = 1,Op_dim - do m = 1,Op_dim - j = Op_V(n_op,nf)%P(m) - do i = 1,Ndim - xp_v(i,n) = xp_v(i,n) + gr(i,j,nf)*x_v(j,n) - enddo - enddo - enddo - - do n = 1,Op_dim - do j = 1,Ndim - do i = 1,Ndim - gr(i,j,nf) = gr(i,j,nf) - xp_v(i,n)*y_v(j,n) - enddo - enddo - enddo - enddo - - ! Flip the spin - nsigma(n_op,nt) = ns_new - endif - - Call Control_upgrade(Log) - - - End Subroutine Upgrade diff --git a/Prog_7/wrapgrdo.f90 b/Prog_7/wrapgrdo.f90 deleted file mode 100644 index d167cab49..000000000 --- a/Prog_7/wrapgrdo.f90 +++ /dev/null @@ -1,82 +0,0 @@ - SUBROUTINE WRAPGRDO(GR,NTAU,PHASE) - - Use Hamiltonian - Use MyMats - Use Hop_mod - Implicit None - - Interface - Subroutine Upgrade(GR,N_op,NT,PHASE,Op_dim) - Use Hamiltonian - Implicit none - Complex (Kind=8) :: GR(Ndim,Ndim, N_FL) - Integer, INTENT(IN) :: N_op, Nt, Op_dim - Complex (Kind=8) :: Phase - End Subroutine Upgrade - End Interface - - ! Given GREEN at time NTAU => GREEN at time NTAU - 1, - ! Upgrade NTAU [LTROT:1] - - COMPLEX (Kind=8), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) - COMPLEX (Kind=8), INTENT(INOUT) :: PHASE - Integer :: NTAU - - ! Local - Complex (Kind=8) :: Mat_TMP(Ndim,Ndim) - Integer :: nf, N_Type, n, I,J - real (Kind=8) :: spin - - Do n = size(Op_V,1), 1, -1 - N_type = 2 - nf = 1 - spin = Phi(nsigma(n,ntau),Op_V(n,nf)%type) - do nf = 1,N_FL - Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type) - enddo - !Write(6,*) 'Upgrade : ', ntau,n - Call Upgrade(GR,n,ntau,PHASE,Op_V(n,1)%N) - ! The spin has changed after the upgrade! - nf = 1 - spin = Phi(nsigma(n,ntau),Op_V(n,nf)%type) - N_type = 1 - do nf = 1,N_FL - Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type ) - enddo - enddo - DO nf = 1,N_FL - Call Hop_mod_mmthl (GR(:,:,nf), MAT_TMP, nf) - Call Hop_mod_mmthr_m1(MAT_TMP, GR(:,:,nf), nf) - !CALL MMULT(MAT_TMP , GR(:,:,nf) , Exp_T(:,:,nf) ) - !CALL MMULT(GR(:,:,nf), Exp_T_M1(:,:,nf), MAT_TMP ) - enddo - -!!$ ! Test -!!$ Mat_TMP = cmplx(0.d0,0.d0) -!!$ DO I = 1,Ndim -!!$ Mat_TMP(I,I) = cmplx(1.d0,0.d0) -!!$ Enddo -!!$ Do n = size(Op_V,1), 1, -1 -!!$ N_type = 2 -!!$ nf = 1 -!!$ spin = Phi(nsigma(n,ntau),Op_V(n,nf)%type) -!!$ Write(6,*) n, spin -!!$ do nf = 1,N_FL -!!$ Call Op_Wrapdo( Mat_tmp, Op_V(n,nf), spin, Ndim, N_Type) -!!$ enddo -!!$ !Upgrade -!!$ N_type = 1 -!!$ do nf = 1,N_FL -!!$ Call Op_Wrapdo( Mat_tmp, Op_V(n,nf), spin, Ndim, N_Type ) -!!$ enddo -!!$ enddo -!!$ -!!$ DO I = 1,Ndim -!!$ Do J = 1,NDIM -!!$ WRITE(6,*) I,J, Mat_tmp(I,J) -!!$ ENDDO -!!$ ENDDO -!!$ -!!$ STOP - - END SUBROUTINE WRAPGRDO diff --git a/Prog_7/wrapgrup.f90 b/Prog_7/wrapgrup.f90 deleted file mode 100644 index 9743dd137..000000000 --- a/Prog_7/wrapgrup.f90 +++ /dev/null @@ -1,53 +0,0 @@ - SUBROUTINE WRAPGRUP(GR,NTAU,PHASE) - - Use Hamiltonian - Use Hop_mod - Implicit none - - Interface - Subroutine Upgrade(GR,N_op,NT,PHASE,Op_dim) - Use Hamiltonian - Implicit none - Complex (Kind=8) :: GR(Ndim,Ndim, N_FL) - Integer, INTENT(IN) :: N_op, Nt, Op_dim - Complex (Kind=8) :: Phase - End Subroutine Upgrade - End Interface - - ! Given GRUP at time NTAU => GRUP at time NTAU + 1. - ! Upgrade NTAU + 1 NTAU: [0:LTROT-1] - - ! Arguments - COMPLEX (Kind=8), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) - COMPLEX (Kind=8), INTENT(INOUT) :: PHASE - INTEGER, INTENT(IN) :: NTAU - - !Local - Integer :: nf, N_Type, NTAU1,n - Complex (Kind=8) :: Mat_TMP(Ndim,Ndim) - Real (Kind=8) :: X - - ! Wrap up, upgrade ntau1. with B^{1}(tau1) - NTAU1 = NTAU + 1 - Do nf = 1,N_FL - CALL HOP_MOD_mmthr( GR(:,:,nf), MAT_TMP,nf) - CALL HOP_MOD_mmthl_m1(MAT_TMP,GR(:,:,nf), nf ) - !CALL MMULT ( MAT_TMP, Exp_T(:,:,nf), GR(:,:,nf) ) - !CALL MMULT ( GR(:,:,nf), MAT_TMP , Exp_T_M1(:,:,nf) ) - Enddo - Do n = 1,Size(Op_V,1) - Do nf = 1, N_FL - X = Phi(nsigma(n,ntau1),Op_V(n,nf)%type) - N_type = 1 - Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),X,Ndim,N_Type) - enddo - nf = 1 - !Write(6,*) 'Upgrade: ', ntau1,n - Call Upgrade(GR,N,ntau1,PHASE,Op_V(n,nf)%N) - do nf = 1,N_FL - N_type = 2 - Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),X,Ndim,N_Type) - enddo - Enddo - - END SUBROUTINE WRAPGRUP diff --git a/Prog_8/Compile_Hub b/Prog_8/Compile_Hub deleted file mode 100644 index ace754b8c..000000000 --- a/Prog_8/Compile_Hub +++ /dev/null @@ -1,16 +0,0 @@ -TARGET= Hubb.out -OBJS= control_mod.o Operator.o print_bin_mod.o Hamiltonian_Hub.o Hop_mod.o UDV_WRAP.o tau_m.o main.o wrapul.o cgr1.o wrapgrup.o wrapur.o upgrade.o \ - nranf.o wrapgrdo.o outconfc.o inconfc.o cgr2.o cgr2_2.o cgr2_1.o - - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Prog_8/Compile_Ising b/Prog_8/Compile_Ising deleted file mode 100644 index c0400f6e8..000000000 --- a/Prog_8/Compile_Ising +++ /dev/null @@ -1,15 +0,0 @@ -TARGET= Ising.out -OBJS= control_mod.o Operator.o print_bin_mod.o Hamiltonian_Ising.o Hop_mod.o UDV_WRAP.o tau_m.o main.o wrapul.o cgr1.o wrapgrup.o wrapur.o upgrade.o \ - nranf.o wrapgrdo.o outconfc.o inconfc.o cgr2.o cgr2_2.o cgr2_1.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Prog_8/Compile_SPT b/Prog_8/Compile_SPT deleted file mode 100644 index e7b63f777..000000000 --- a/Prog_8/Compile_SPT +++ /dev/null @@ -1,20 +0,0 @@ -TARGET= SPT.out -OBJS= control_mod.o Operator.o print_bin_mod.o Hamiltonian_SPT.o Hop_mod.o UDV_WRAP.o tau_m.o main.o wrapul.o cgr1.o wrapgrup.o wrapur.o upgrade.o \ - nranf.o wrapgrdo.o outconfc.o inconfc.o cgr2.o cgr2_2.o cgr2_1.o - -#block.o block_obs.o Mol_Dyn.o Hubb.o inconfc.o outconfc.o npbc.o salph.o sli.o \ -# sthop.o wrapul.o wrapur.o cgr1.o \ -# wrapgrdo.o mmthr.o mmthrm1.o mmthl.o mmthlm1.o obser.o upgradeu.o wrapgrup.o \ -# preq.o cgr2.o propr.o proprm1.o obsert.o prtau.o tau_m.o set_Hopping.o - - -.SUFFIXES: .f90 .f -.f.o .f90.o: - $(FC) -c -cpp -o $@ $(FLAGS) $< - -$(TARGET): $(OBJS) - $(FC) $(LF) -o $(TARGET) $(OBJS) $(LIBS) - -clean: - rm $(OBJS) - diff --git a/Prog_8/Hamiltonian_Hub.f90 b/Prog_8/Hamiltonian_Hub.f90 deleted file mode 100644 index 72a1d5f17..000000000 --- a/Prog_8/Hamiltonian_Hub.f90 +++ /dev/null @@ -1,539 +0,0 @@ - Module Hamiltonian - - Use Operator_mod - Use Lattices_v3 - Use MyMats - Use Random_Wrap - Use Files_mod - Use Matrix - Use Print_bin_mod - - - Type (Operator), dimension(:,:), allocatable :: Op_V - Type (Operator), dimension(:,:), allocatable :: Op_T - Integer, allocatable :: nsigma(:,:) - Integer :: Ndim, N_FL, N_SUN, Ltrot - !Complex (Kind=8), dimension(:,:,:), allocatable :: Exp_T(:,:,:), Exp_T_M1(:,:,:) - - ! ToDo. Public and private subroutines. - - ! What is below is private - Type (Lattice), private :: Latt - Integer, private :: L1, L2 - real (Kind=8), private :: ham_T , ham_U, Ham_chem - real (Kind=8), private :: Dtau, Beta - Character (len=64), private :: Model, Lattice_type - Logical, private :: One_dimensional - Integer, private :: N_coord - - - ! Observables - Integer, private :: Nobs, Norb - Complex (Kind=8), allocatable, private :: obs_scal(:) - Complex (Kind=8), allocatable, private :: Green_eq (:,:,:), SpinZ_eq (:,:,:), SpinXY_eq (:,:,:), & - & Den_eq(:,:,:) - Complex (Kind=8), allocatable, private :: Green_eq0 (:), SpinZ_eq0(:), SpinXY_eq0(:), & - & Den_eq0(:) - - ! For time displaced - Integer, private :: NobsT - Complex (Kind=8), private :: Phase_tau - Complex (Kind=8), allocatable, private :: Green_tau(:,:,:,:), Den_tau(:,:,:,:) - - contains - - - Subroutine Ham_Set - - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - integer :: ierr - - - NAMELIST /VAR_lattice/ L1, L2, Lattice_type, Model - - NAMELIST /VAR_Hubbard/ ham_T, ham_chem, ham_U, Dtau, Beta - - -#ifdef MPI - Integer :: Isize, Irank - Integer :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - - ! NAMELIST /VAR_Model/ N_FL, N_SUN, ham_T , ham_xi, ham_h, ham_J, ham_U, Ham_Vint, & - ! & Dtau, Beta - - -#ifdef MPI - If (Irank == 0 ) then -#endif - OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) - IF (ierr /= 0) THEN - WRITE(*,*) 'unable to open ',ierr - STOP - END IF - READ(5,NML=VAR_lattice) - CLOSE(5) - -#ifdef MPI - Endif - CALL MPI_BCAST(L1 ,1 ,MPI_INTEGER, 0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(L2 ,1 ,MPI_INTEGER, 0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Model ,64 ,MPI_CHARACTER, 0,MPI_COMM_WORLD,IERR) - CALL MPI_BCAST(Lattice_type,64 ,MPI_CHARACTER, 0,MPI_COMM_WORLD,IERR) -#endif - Call Ham_latt - - If ( Model == "Hubbard_Mz") then - N_FL = 2 - N_SUN = 1 - elseif ( Model == "Hubbard_SU2" ) then - N_FL = 1 - N_SUN = 2 - else - Write(6,*) "Model not yet implemented!" - Stop - endif -#ifdef MPI - If (Irank == 0 ) then -#endif - OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) - READ(5,NML=VAR_Hubbard) - CLOSE(5) -#ifdef MPI - endif - CALL MPI_BCAST(ham_T ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(ham_chem ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(ham_U ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Dtau ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Beta ,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr) -#endif - Call Ham_hop - - Ltrot = nint(beta/dtau) -#ifdef MPI - If (Irank == 0) then -#endif - Open (Unit = 50,file="info",status="unknown",position="append") - Write(50,*) '=====================================' - Write(50,*) 'Model is : ', Model - Write(50,*) 'Beta : ', Beta - Write(50,*) 'dtau,Ltrot : ', dtau,Ltrot - Write(50,*) 'N_SUN : ', N_SUN - Write(50,*) 'N_FL : ', N_FL - Write(50,*) 't : ', Ham_T - Write(50,*) 'Ham_U : ', Ham_U - Write(50,*) 'Ham_chem : ', Ham_chem - close(50) -#ifdef MPI - endif -#endif - call Ham_V - end Subroutine Ham_Set -!============================================================================= - Subroutine Ham_Latt - Implicit none - !Set the lattice - Real (Kind=8) :: a1_p(2), a2_p(2), L1_p(2), L2_p(2) - If ( Lattice_type =="Square" ) then - a1_p(1) = 1.0 ; a1_p(2) = 0.d0 - a2_p(1) = 0.0 ; a2_p(2) = 1.d0 - L1_p = dble(L1)*a1_p - L2_p = dble(L2)*a2_p - Call Make_Lattice( L1_p, L2_p, a1_p, a2_p, Latt ) - Ndim = Latt%N - !Write(6,*) 'Lattice: ', Ndim - One_dimensional = .false. - N_coord = 2 - If ( L1 == 1 .or. L2 == 1 ) then - One_dimensional = .true. - N_coord = 1 - If (L1 == 1 ) then - Write(6,*) ' For one dimensional systems set L2 = 1 ' - Stop - endif - endif - else - Write(6,*) "Lattice not yet implemented!" - Stop - endif - end Subroutine Ham_Latt - -!=================================================================================== - Subroutine Ham_hop - Implicit none - - !Setup the hopping - !Per flavor, the hopping is given by: - ! e^{-dtau H_t} = Prod_{n=1}^{Ncheck} e^{-dtau_n H_{n,t}} - - - Integer :: I, I1, I2, n, Ncheck,nc - Real (Kind=8) :: X - - Ncheck = 1 - allocate(Op_T(Ncheck,N_FL)) - do n = 1,N_FL - Do nc = 1,Ncheck - Call Op_make(Op_T(nc,n),Ndim) - If (One_dimensional ) then - DO I = 1, Latt%N - I1 = Latt%nnlist(I,1,0) - Op_T(nc,n)%O(I,I1) = cmplx(-Ham_T,0.d0) - Op_T(nc,n)%O(I1,I) = cmplx(-Ham_T,0.d0) - Op_T(nc,n)%O(I ,I) = cmplx(-Ham_chem,0.d0) - ENDDO - else - DO I = 1, Latt%N - I1 = Latt%nnlist(I,1,0) - I2 = Latt%nnlist(I,0,1) - Op_T(nc,n)%O(I,I1) = cmplx(-Ham_T, 0.d0) - Op_T(nc,n)%O(I1,I) = cmplx(-Ham_T, 0.d0) - Op_T(nc,n)%O(I,I2) = cmplx(-Ham_T, 0.d0) - Op_T(nc,n)%O(I2,I) = cmplx(-Ham_T, 0.d0) - Op_T(nc,n)%O(I ,I) = cmplx(-Ham_chem,0.d0) - ENDDO - endif - - Do I = 1,Latt%N - Op_T(nc,n)%P(i) = i - Enddo - if ( abs(Ham_T) < 1.E-6 .and. abs(Ham_chem) < 1.E-6 ) then - Op_T(nc,n)%g=cmplx(0.d0 ,0.d0) - else - Op_T(nc,n)%g=cmplx(-Dtau,0.d0) - endif - Op_T(nc,n)%alpha=cmplx(0.d0,0.d0) - !Write(6,*) 'In Ham_hop', Ham_T - Call Op_set(Op_T(nc,n)) - !Write(6,*) 'In Ham_hop 1' - !Do I = 1,Latt%N - ! Write(6,*) Op_T(n)%E(i) - !enddo - !Call Op_exp( cmplx(-Dtau,0.d0), Op_T(n), Exp_T (:,:,n) ) - !Call Op_exp( cmplx( Dtau,0.d0), Op_T(n), Exp_T_M1(:,:,n) ) - enddo - enddo - end Subroutine Ham_hop - -!=================================================================================== - - Subroutine Ham_V - - Implicit none - - Integer :: nf, nth, n, n1, n2, n3, n4, I, I1, I2, J, Ix, Iy, nc - Real (Kind=8) :: X_p(2), X1_p(2), X2_p(2), X - - - If (Model == "Hubbard_SU2") then - !Write(50,*) 'Model is ', Model - Allocate(Op_V(Latt%N,N_FL)) - do nf = 1,N_FL - do i = 1, Latt%N - Call Op_make(Op_V(i,nf),1) - enddo - enddo - Do nf = 1,N_FL - nc = 0 - Do i = 1,Latt%N - nc = nc + 1 - Op_V(nc,nf)%P(1) = I - Op_V(nc,nf)%O(1,1) = cmplx(1.d0 ,0.d0) - Op_V(nc,nf)%g = SQRT(CMPLX(-DTAU*ham_U/(DBLE(N_SUN)),0.D0)) - Op_V(nc,nf)%alpha = cmplx(-0.5d0,0.d0) - Op_V(nc,nf)%type = 2 - Call Op_set( Op_V(nc,nf) ) - ! The operator reads: - ! g*s*( c^{dagger} O c + alpha )) - ! with s the HS field. - Enddo - Enddo - Elseif (Model == "Hubbard_Mz") then - !Write(50,*) 'Model is ', Model - Allocate(Op_V(Latt%N,N_FL)) - do nf = 1,N_FL - do i = 1, Latt%N - Call Op_make(Op_V(i,nf),1) - enddo - enddo - Do nf = 1,N_FL - nc = 0 - X = 1.d0 - if (nf == 2) X = -1.d0 - Do i = 1,Latt%N - nc = nc + 1 - Op_V(nc,nf)%P(1) = I - Op_V(nc,nf)%O(1,1) = cmplx(1.d0 ,0.d0) - Op_V(nc,nf)%g = X*SQRT(CMPLX(DTAU*ham_U/2.d0 ,0.D0)) - Op_V(nc,nf)%alpha = cmplx(0.d0,0.d0) - Op_V(nc,nf)%type = 2 - Call Op_set( Op_V(nc,nf) ) - ! The operator reads: - ! g*s*( c^{dagger} O c - alpha )) - ! with s the HS field. - ! Write(6,*) nc,nf, Op_V(nc,nf)%g - Enddo - Enddo - Endif - end Subroutine Ham_V - -!=================================================================================== - Real (Kind=8) function S0(n,nt) - Implicit none - Integer, Intent(IN) :: n,nt - Integer :: i, nt1 - S0 = 1.d0 - - end function S0 - -!=================================================================================== - Subroutine Alloc_obs(Ltau) - - Implicit none - Integer, Intent(In) :: Ltau - Integer :: I - - Allocate ( Obs_scal(5) ) - Allocate ( Green_eq(Latt%N,1,1), SpinZ_eq(Latt%N,1,1), SpinXY_eq(Latt%N,1,1), & - & Den_eq(Latt%N,1,1) ) - Allocate ( Green_eq0(1), SpinZ_eq0(1), SpinXY_eq0(1), Den_eq0(1) ) - - - If (Ltau == 1) then - Allocate ( Green_tau(Latt%N,Ltrot+1,1,1), Den_tau(Latt%N,Ltrot+1,1,1) ) - endif - - end Subroutine Alloc_obs - -!=================================================================================== - - Subroutine Init_obs(Ltau) - - Implicit none - Integer, Intent(In) :: Ltau - - Integer :: I,n - - Nobs = 0 - Obs_scal = cmplx(0.d0,0.d0) - Green_eq = cmplx(0.d0,0.d0) - SpinZ_eq = cmplx(0.d0,0.d0) - SpinXY_eq = cmplx(0.d0,0.d0) - Den_eq = cmplx(0.d0,0.d0) - Green_eq0 = cmplx(0.d0,0.d0) - SpinZ_eq0 = cmplx(0.d0,0.d0) - SpinXY_eq0= cmplx(0.d0,0.d0) - Den_eq0 = cmplx(0.d0,0.d0) - - - If (Ltau == 1) then - NobsT = 0 - Phase_tau = cmplx(0.d0,0.d0) - Green_tau = cmplx(0.d0,0.d0) - Den_tau = cmplx(0.d0,0.d0) - endif - - end Subroutine Init_obs - -!======================================================================== - Subroutine Obser(GR,Phase,Ntau) - - Implicit none - - Complex (Kind=8), INTENT(IN) :: GR(Ndim,Ndim,N_FL) - Complex (Kind=8), Intent(IN) :: PHASE - Integer, INTENT(IN) :: Ntau - - !Local - Complex (Kind=8) :: GRC(Ndim,Ndim,N_FL), ZK - Complex (Kind=8) :: Zrho, Zkin, ZPot, Z, ZP,ZS - Integer :: I,J, no,no1, n, n1, imj, nf, I1, I2, J1, J2 - - Real (Kind=8) :: G(4,4), X, FI, FJ - - Nobs = Nobs + 1 - ZP = PHASE/cmplx(Real(Phase,kind=8),0.d0) - ZS = cmplx(Real(Phase,kind=8)/Abs(Real(Phase,kind=8)), 0.d0) - - - Do nf = 1,N_FL - Do I = 1,Ndim - Do J = 1,Ndim - ZK = cmplx(0.d0,0.d0) - If ( I == J ) ZK = cmplx(1.d0,0.d0) - GRC(I,J,nf) = ZK - GR(J,I,nf) - Enddo - Enddo - Enddo - ! GRC(i,j,nf) = < c^{dagger}_{j,nf } c_{j,nf } > - ! Compute scalar observables. - Zkin = cmplx(0.d0,0.d0) - Do nf = 1,N_FL - Do J = 1,Ndim - DO I = 1,Ndim - Zkin = Zkin + Op_T(1,nf)%O(i,j)*Grc(i,j,nf) - Enddo - ENddo - Enddo - Zkin = Zkin*cmplx( dble(N_SUN), 0.d0 ) - - Zrho = cmplx(0.d0,0.d0) - Do nf = 1,N_FL - Do I = 1,Ndim - Zrho = Zrho + Grc(i,i,nf) - enddo - enddo - Zrho = Zrho*cmplx( dble(N_SUN), 0.d0 ) - - ZPot = cmplx(0.d0,0.d0) - If ( Model == "Hubbard_SU2" ) then - Do I = 1,Ndim - ZPot = ZPot + Grc(i,i,1) * Grc(i,i,1) - Enddo - Zpot = Zpot*cmplx(ham_U,0.d0) - elseif ( Model == "Hubbard_Mz" ) then - Do I = 1,Ndim - ZPot = ZPot + Grc(i,i,1) * Grc(i,i,2) - Enddo - Zpot = Zpot*cmplx(ham_U,0.d0) - endif - - Obs_scal(1) = Obs_scal(1) + zrho * ZP*ZS - Obs_scal(2) = Obs_scal(2) + zkin * ZP*ZS - Obs_scal(3) = Obs_scal(3) + Zpot * ZP*ZS - Obs_scal(4) = Obs_scal(4) + (zkin + Zpot)*ZP*ZS - Obs_scal(5) = Obs_scal(5) + ZS - ! You will have to allocate more space if you want to include more scalar observables. - - ! Compute spin-spin, Green, and den-den correlation functions ! This is general N_SUN, and N_FL = 1 - If ( Model == "Hubbard_SU2" ) then - Z = cmplx(dble(N_SUN),0.d0) - Do I = 1,Latt%N - Do J = 1,Latt%N - imj = latt%imj(I,J) - GREEN_EQ (imj,1,1) = GREEN_EQ (imj,1,1) + Z * GRC(I,J,1) * ZP*ZS - SPINXY_Eq (imj,1,1) = SPINXY_Eq (imj,1,1) + Z * GRC(I,J,1) * GR(I,J,1) * ZP*ZS - SPINZ_Eq (imj,1,1) = SPINZ_Eq (imj,1,1) + Z * GRC(I,J,1) * GR(I,J,1) * ZP*ZS - DEN_Eq (imj,1,1) = DEN_Eq (imj,1,1) + ( & - & GRC(I,I,1) * GRC(J,J,1) *Z + & - & GRC(I,J,1) * GR(I,J,1) & - & ) * Z* ZP*ZS - ENDDO - Den_eq0(1) = Den_eq0(1) + Z * GRC(I,I,1) * ZP * ZS - ENDDO - elseif (Model == "Hubbard_Mz" ) Then - DO I = 1,Latt%N - DO J = 1, Latt%N - imj = latt%imj(I,J) - SPINZ_Eq (imj,1,1) = SPINZ_Eq (imj,1,1) + & - & ( GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2) + & - & (GRC(I,I,2) - GRC(I,I,1))*(GRC(J,J,2) - GRC(J,J,1)) ) * ZP*ZS - ! c^d_(i,u) c_(i,d) c^d_(j,d) c_(j,u) + c^d_(i,d) c_(i,u) c^d_(j,u) c_(j,d) - SPINXY_Eq (imj,1,1) = SPINXY_Eq (imj,1,1) + & - & ( GRC(I,J,1) * GR(I,J,2) + GRC(I,J,2) * GR(I,J,1) ) * ZP*ZS - - DEN_Eq (imj,1,1) = DEN_Eq (imj,1,1) + & - & ( GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2) + & - & (GRC(I,I,2) + GRC(I,I,1))*(GRC(J,J,2) + GRC(J,J,1)) ) * ZP*ZS - enddo - Den_eq0(1) = Den_eq0(1) + (GRC(I,I,2) + GRC(I,I,1)) * ZP*ZS - enddo - Endif - - - end Subroutine Obser -!========================================================== - Subroutine Pr_obs(LTAU) - - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - - Integer, Intent(In) :: Ltau - - Character (len=64) :: File_pr - Complex (Kind=8) :: Phase_bin -#ifdef MPI - Integer :: Isize, Irank, Ierr - Integer :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'In Pr_obs', LTAU -!!$#else -!!$ Write(6,*) 'In Pr_obs', LTAU -!!$#endif - - Phase_bin = Obs_scal(5)/cmplx(dble(Nobs),0.d0) - File_pr ="SpinZ_eq" - Call Print_bin(SpinZ_eq ,SpinZ_eq0, Latt, Nobs, Phase_bin, file_pr) - File_pr ="SpinXY_eq" - Call Print_bin(Spinxy_eq, Spinxy_eq0,Latt, Nobs, Phase_bin, file_pr) - File_pr ="Den_eq" - Call Print_bin(Den_eq, Den_eq0, Latt, Nobs, Phase_bin, file_pr) - File_pr ="Green_eq" - Call Print_bin(Green_eq , Green_eq0 ,Latt, Nobs, Phase_bin, file_pr) - - File_pr ="ener" - Call Print_scal(Obs_scal, Nobs, file_pr) - - If (Ltau == 1) then - Phase_tau = Phase_tau/cmplx(dble(NobsT),0.d0) - File_pr = "Green_tau" - Call Print_bin_tau(Green_tau,Latt,NobsT,Phase_tau, file_pr,dtau) - File_pr = "Den_tau" - Call Print_bin_tau(Den_tau,Latt,NobsT,Phase_tau, file_pr,dtau) - endif - -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'out Pr_obs', LTAU -!!$#else -!!$ Write(6,*) 'out Pr_obs', LTAU -!!$#endif - end Subroutine Pr_obs -!========================================================== - - Subroutine OBSERT(NT, GT0,G0T,G00,GTT, PHASE) - Implicit none - - Integer , INTENT(IN) :: NT - Complex (Kind=8), INTENT(IN) :: GT0(Ndim,Ndim,N_FL),G0T(Ndim,Ndim,N_FL),G00(Ndim,Ndim,N_FL),GTT(Ndim,Ndim,N_FL) - Complex (Kind=8), INTENT(IN) :: Phase - - !Locals - Complex (Kind=8) :: Z, ZP, ZS - Integer :: IMJ, I, J - - ZP = PHASE/cmplx(Real(Phase,kind=8),0.d0) - ZS = cmplx(Real(Phase,kind=8)/Abs(Real(Phase,kind=8)), 0.d0) - If (NT == 0 ) then - Phase_tau = Phase_tau + ZS - NobsT = NobsT + 1 - endif - If ( N_FL == 1 ) then - Z = cmplx(dble(N_SUN),0.d0) - Do I = 1,Latt%N - Do J = 1,Latt%N - imj = latt%imj(I,J) - Green_tau(imj,nt+1,1,1) = green_tau(imj,nt+1,1,1) + Z * GT0(I,J,1) * ZP* ZS - Den_tau (imj,nt+1,1,1) = Den_tau (imj,nt+1,1,1) - Z * GT0(I,J,1)*G0T(J,I,1) * ZP* ZS - Enddo - Enddo - Endif - end Subroutine OBSERT - - - end Module Hamiltonian diff --git a/Prog_8/Hop_mod.f90 b/Prog_8/Hop_mod.f90 deleted file mode 100644 index 0d6028fb4..000000000 --- a/Prog_8/Hop_mod.f90 +++ /dev/null @@ -1,217 +0,0 @@ -! This is for the Kondo project with tarun. - Module Hop_mod - - - Use Hamiltonian - Use Random_wrap - Use MyMats - - ! Private variables - Complex (Kind=8), allocatable, private :: Exp_T(:,:,:,:), Exp_T_M1(:,:,:,:) - Complex (Kind=8), allocatable, private :: U_HLP(:,:), U_HLP1(:,:), V_HLP(:,:), V_HLP1(:,:) - Integer, private, save :: Ncheck, Ndim_hop - Real (Kind=8), private, save :: Zero - - Contains - - subroutine Hop_mod_init - - Implicit none - - Integer :: nc, nf - Complex (Kind=8) :: g - - Ncheck = size(Op_T,1) - If ( size(Op_T,2) /= N_FL ) then - Write(6,*) 'Error in the number of flavors.' - Stop - Endif - Ndim_hop = Op_T(1,1)%N - Write(6,*) 'In Hop_mod: ', Ndim, Ndim_hop, Ncheck - Do nc = 1, Ncheck - do nf = 1,N_FL - if ( Ndim_hop /= Op_T(nc,nf)%N ) Then - Write(6,*) 'Different size of Hoppings not implemented ' - Stop - endif - enddo - enddo - - Allocate ( Exp_T (Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( Exp_T_M1(Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( V_Hlp(Ndim_hop,Ndim) ) - Allocate ( V_Hlp1(Ndim_hop,Ndim) ) - Allocate ( U_Hlp (Ndim, Ndim_hop) ) - Allocate ( U_Hlp1(Ndim, Ndim_hop) ) - - Exp_T = cmplx(0.d0,0.d0) - Exp_T_M1 = cmplx(0.d0,0.d0) - do nf = 1,N_FL - do nc = 1,Ncheck - g = Op_T(nc,nf)%g - Call Op_exp(g,Op_T(nc,nf),Exp_T(:,:,nc,nf)) - g = -Op_T(nc,nf)%g - Call Op_exp(g,Op_T(nc,nf),Exp_T_M1(:,:,nc,nf)) - enddo - enddo - - Zero = 1.E-12 - - end subroutine Hop_mod_init - -!============================================================================ - Subroutine Hop_mod_mmthr(In, Out,nf) - - - ! In: IN - ! Out: OUT = e^{ -dtau T }.IN - Implicit none - - Complex (Kind=8), intent(IN) :: IN(Ndim,Ndim) - Complex (Kind=8), intent(INOUT) :: OUT(Ndim,Ndim) - Integer :: nf - - !Local - Integer :: nc, I, n - - Out = In - do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - do I = 1,Ndim - do n = 1,Ndim_hop - V_Hlp(n,I) = Out(Op_T(nc,nf)%P(n),I) - enddo - enddo - Call mmult(V_HLP1,Exp_T(:,:,nc,nf),V_Hlp) - DO I = 1,Ndim - do n = 1,Ndim_hop - OUT(OP_T(nc,nf)%P(n),I) = V_hlp1(n,I) - Enddo - Enddo - Endif - Enddo - - end Subroutine Hop_mod_mmthr - -!============================================================================ - Subroutine Hop_mod_mmthr_m1(In, Out,nf) - - - ! In: IN - ! Out: OUT = e^{ dtau T }.IN - Implicit none - - Complex (Kind=8), intent(IN) :: IN(Ndim,Ndim) - Complex (Kind=8), intent(INOUT) :: OUT(Ndim,Ndim) - Integer :: nf - - !Local - Integer :: nc, I, n - - - Out = In - do nc = 1,Ncheck - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - do I = 1,Ndim - do n = 1,Ndim_hop - V_Hlp(n,I) = Out(Op_T(nc,nf)%P(n),I) - enddo - enddo - Call mmult(V_HLP1,Exp_T_m1(:,:,nc,nf),V_Hlp) - DO I = 1,Ndim - do n = 1,Ndim_hop - OUT(OP_T(nc,nf)%P(n),I) = V_hlp1(n,I) - Enddo - Enddo - Endif - Enddo - - end Subroutine Hop_mod_mmthr_m1 - -!============================================================================ - Subroutine Hop_mod_mmthl (In, Out,nf) - - - ! In: IN - ! Out: OUT = IN * e^{ -dtau T } - Implicit none - - Complex (Kind=8), intent(IN) :: IN(Ndim,Ndim) - Complex (Kind=8), intent(INOUT) :: OUT(Ndim,Ndim) - Integer :: nf - - !Local - Integer :: nc, I, n - - Out = In - do nc = 1, Ncheck - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - do n = 1,Ndim_hop - do I = 1,Ndim - U_Hlp(I,n) = Out(I,Op_T(nc,nf)%P(n)) - enddo - enddo - Call mmult(U_Hlp1,U_Hlp,Exp_T(:,:,nc,nf)) - do n = 1,Ndim_hop - DO I = 1,Ndim - OUT(I,OP_T(nc,nf)%P(n)) = U_hlp1(I,n) - Enddo - Enddo - Endif - Enddo - - end Subroutine Hop_mod_mmthl -!============================================================================ - Subroutine Hop_mod_mmthl_m1 (In, Out,nf) - - - ! In: IN - ! Out: OUT = IN * e^{ dtau T } - Implicit none - - Complex (Kind=8), intent(IN) :: IN(Ndim,Ndim) - Complex (Kind=8), intent(INOUT) :: OUT(Ndim,Ndim) - Integer :: nf - - !Local - Integer :: nc, I, n - - Out = In - do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - do n = 1,Ndim_hop - do I = 1,Ndim - U_Hlp(I,n) = Out(I,Op_T(nc,nf)%P(n)) - enddo - enddo - Call mmult(U_Hlp1,U_Hlp,Exp_T_M1(:,:,nc,nf)) - do n = 1,Ndim_hop - DO I = 1,Ndim - OUT(I,OP_T(nc,nf)%P(n)) = U_hlp1(I,n) - Enddo - Enddo - Endif - Enddo - - end Subroutine Hop_mod_mmthl_m1 - -!============================================================================ -!!$ Subroutine Hop_mod_test -!!$ -!!$ Implicit none -!!$ -!!$ Complex (Kind=8) :: IN(Ndim,Ndim),Out(Ndim,Ndim) -!!$ Complex (Kind=8) :: Test(Ndim,Ndim) -!!$ -!!$ Integer :: I,J -!!$ -!!$ DO I = 1,Ndim -!!$ DO J = 1,Ndim -!!$ IN(J,I) = cmplx(Ranf(),Ranf()) -!!$ ENDDO -!!$ ENDDO -!!$ -!!$ !Write(6,*) IN -!!$ end Subroutine Hop_mod_test - - end Module Hop_mod diff --git a/Prog_8/UDV_WRAP.f90 b/Prog_8/UDV_WRAP.f90 deleted file mode 100644 index 4fb367afb..000000000 --- a/Prog_8/UDV_WRAP.f90 +++ /dev/null @@ -1,135 +0,0 @@ - Module UDV_Wrap_mod - Use MyMats - Use Files_mod - - Contains - -!*************************************************************** - Subroutine UDV_Wrap_Pivot(A,U,D,V,NCON,N1,N2) - - Implicit NONE - COMPLEX (KIND=8), INTENT(IN), DIMENSION(:,:) :: A - COMPLEX (KIND=8), INTENT(INOUT), DIMENSION(:,:) :: U,V - COMPLEX (KIND=8), INTENT(INOUT), DIMENSION(:) :: D - INTEGER, INTENT(IN) :: NCON - INTEGER, INTENT(IN) :: N1,N2 - - ! Locals - REAL (Kind=8) :: VHELP(N2), XNORM(N2), XMAX, XMEAN - INTEGER :: IVPT(N2), IVPTM1(N2), I, J, K, IMAX - COMPLEX (KIND=8) :: A1(N1,N2), A2(N1,N2) - - DO I = 1,N2 - XNORM(I) = 0.D0 - DO J = 1,N1 - XNORM(I) = XNORM(I) + DBLE( A(J,I) * CONJG( A(J,I) ) ) - ENDDO - ENDDO - DO I = 1,N2 - VHELP(I) = XNORM(I) - ENDDO - - DO I = 1,N2 - XMAX = 0.D0 - DO J = 1,N2 - IF (VHELP(J).GT.XMAX) THEN - IMAX = J - XMAX = VHELP(J) - ENDIF - ENDDO - VHELP(IMAX) = -1.D0 - IVPTM1(IMAX)= I - IVPT(I) = IMAX - ENDDO - DO I = 1,N2 - K = IVPT(I) - DO J = 1,N1 - A1(J,I) = A(J,K) - ENDDO - ENDDO - - CALL UDV_Wrap(A1,U,D,V,NCON) - - A1 = V - DO I = 1,N2 - K = IVPTM1(I) - DO J = 1,N1 - V(J,I) = A1(J,K) - ENDDO - ENDDO - - - IF (NCON == 1) THEN - !Check the result A = U D V - DO J = 1,N2 - DO I = 1,N1 - A1(I,J) = D(I)*V(I,J) - ENDDO - ENDDO - Call MMULT (A2,U,A1) - CALL COMPARE(A,A2,XMAX,XMEAN) - Write (6,*) 'Check afer Pivoting', XMAX - ENDIF - - - - End Subroutine UDV_Wrap_Pivot -!*************************************************************** - Subroutine UDV_Wrap(A,U,D,V,NCON) - -#include "machine" - - Implicit None -#ifdef MPI - INCLUDE 'mpif.h' -#endif - COMPLEX (KIND=8), INTENT(IN), DIMENSION(:,:) :: A - COMPLEX (KIND=8), INTENT(INOUT), DIMENSION(:,:) :: U,V - COMPLEX (KIND=8), INTENT(INOUT), DIMENSION(:) :: D - INTEGER, INTENT(IN) :: NCON - - !Local - Complex (Kind=8), Allocatable :: A1(:,:),U1(:,:) - Integer :: I,J, N - character (len=64) :: file_sr, File -#ifdef MPI - INTEGER :: STATUS(MPI_STATUS_SIZE) - INTEGER :: Isize, Irank,Ierr - - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - File_sr = "SDV" -#ifdef MPI - File = File_i(File_sr, Irank) -#else - File = File_sr -#endif - !Open (Unit = 78,File=File, Status='UNKNOWN', action="write", position="append") - !Write(78,*) 'Call QR' - !Close(78) - CALL QR(A,U,V,NCON) - !Open (Unit = 78,File=File, Status='UNKNOWN', action="write", position="append") - !Write(78,*) 'End call QR' - !Close(78) - N = Size(V,1) - Allocate (A1(N,N),U1(N,N)) - A1 = V - !Open (Unit = 78,File=File, Status='UNKNOWN') - !Write(78,*) 'Call SVD' - !DO I = 1,N - ! Write(78,*) Real(V(I,I)) - !ENDDO - !Close(78) - CALL SVD(A1,U1,D,V,NCON) - !Open (Unit = 78,File=File, Status='UNKNOWN', action="write", position="append") - !Write(78,*) 'End call SVD' - !Close(78) - Call MMULT(A1,U,U1) - U = A1 - - End Subroutine UDV_Wrap - - End Module UDV_Wrap_mod - diff --git a/Prog_8/cgr1.f90 b/Prog_8/cgr1.f90 deleted file mode 100644 index b7fc92ca1..000000000 --- a/Prog_8/cgr1.f90 +++ /dev/null @@ -1,110 +0,0 @@ - SUBROUTINE CGR(PHASE,NVAR, GRUP, URUP,DRUP,VRUP, ULUP,DLUP,VLUP) - - Use UDV_Wrap_mod - - Implicit None - - !!! GRUP = (1 + UR*DR*VR*VL*DL*UL)^-1 - !!! NVAR = 1 Big scales are in DL - !!! NVAR = 2 Big scales are in DR - - !Arguments. - COMPLEX(Kind=8), Dimension(:,:), Intent(IN) :: URUP, VRUP, ULUP, VLUP - COMPLEX(Kind=8), Dimension(:), Intent(In) :: DLUP, DRUP - COMPLEX(Kind=8), Dimension(:,:), Intent(INOUT) :: GRUP - COMPLEX(Kind=8) :: PHASE - INTEGER :: NVAR - - !Local - COMPLEX (Kind=8), Dimension(:,:), Allocatable :: UUP, VUP, TPUP, TPUP1, & - & TPUPM1,TPUP1M1, UUPM1, VUP1 - COMPLEX (Kind=8), Dimension(:) , Allocatable :: DUP - COMPLEX (Kind=8) :: ZDUP1, ZDDO1, ZDUP2, ZDDO2, Z1, ZUP, ZDO, Z - Integer :: I,J, N_size, NCON, NR, NT, N - Real (Kind=8) :: X, Xmax - - N_size = SIZE(DLUP,1) - NCON = 0 - - Allocate( UUP(N_size,N_size), VUP(N_size,N_size), TPUP(N_size,N_size), TPUP1(N_size,N_size), & - & TPUPM1(N_size,N_size),TPUP1M1(N_size,N_size), UUPM1(N_size,N_size), VUP1(N_size,N_size), DUP(N_size) ) - - !Write(6,*) 'In CGR', N_size - CALL MMULT(VUP,VRUP,VLUP) - DO J = 1,N_size - DO I = 1,N_size - TPUP(I,J) = DRUP(I)*VUP(I,J)*DLUP(J) - ENDDO - ENDDO - CALL MMULT(UUP,ULUP,URUP) - DO J = 1,N_size - DO I = 1,N_size - UUPM1(I,J) = CONJG(UUP(J,I)) - ENDDO - ENDDO - DO J = 1,N_size - DO I = 1,N_size - TPUP(I,J) = TPUP(I,J) + UUPM1(I,J) - ENDDO - ENDDO - IF (NVAR.EQ.1) THEN - !WRITE(6,*) 'UDV of U + DR * V * DL' - CALL UDV_WRAP(TPUP,UUP,DUP,VUP,NCON) - !CALL UDV(TPUP,UUP,DUP,VUP,NCON) - CALL MMULT(TPUP,VUP,ULUP) - !Do I = 1,N_size - ! Write(6,*) DLUP(I) - !enddo - CALL INV(TPUP,TPUPM1,ZDUP1) - !WRITE(6,*) 'End called Inv' - CALL MMULT(TPUP1,URUP,UUP) - CALL INV(TPUP1,TPUP1M1,ZDUP2) - Z1 = ZDUP1*ZDUP2 - ELSEIF (NVAR.EQ.2) THEN - !WRITE(6,*) 'UDV of (U + DR * V * DL)^{*}' - DO J = 1,N_size - DO I = 1,N_size - TPUP1(I,J) = CONJG( TPUP(J,I) ) - ENDDO - ENDDO - CALL UDV_WRAP(TPUP1,UUP,DUP,VUP,NCON) - !CALL UDV(TPUP1,UUP,DUP,VUP,NCON) - DO J = 1,N_size - DO I = 1,N_size - TPUP(I,J) = CONJG( ULUP(J,I) ) - ENDDO - ENDDO - CALL MMULT(TPUPM1,TPUP,UUP) - DO J = 1,N_size - DO I = 1,N_size - VUP1(I,J) = CONJG( VUP(J,I) ) - ENDDO - ENDDO - CALL MMULT(TPUP1,URUP,VUP1) - CALL INV(TPUP1,TPUP1M1,ZDUP2) - CALL INV(TPUPM1, TPUP, ZDUP1) - Z1 = ZDUP2/ZDUP1 - ENDIF - DO I = 1,N_size - Z = DUP(I) - if (I == 1) Xmax = real(SQRT( Z* conjg(Z)),kind=8) - if ( real(SQRT( Z* conjg(Z)),kind=8) < Xmax ) Xmax = & - & real(SQRT( Z* conjg(Z)),kind=8) - ENDDO - !Write(6,*) 'Cgr1, Cutoff: ', Xmax - - - DO J = 1,N_size - DO I = 1,N_size - ZUP = CMPLX(0.D0,0.D0) - DO NR = 1,N_size - ZUP = ZUP + TPUPM1(I,NR)*TPUP1M1(NR,J)/DUP(NR) - ENDDO - GRUP(I,J) = ZUP - ENDDO - ENDDO - PHASE = Z1/SQRT( Z1* CONJG(Z1) ) - - Deallocate(UUP, VUP, TPUP,TPUP1,TPUPM1, TPUP1M1, UUPM1, VUP1, DUP ) - - END SUBROUTINE CGR diff --git a/Prog_8/cgr2.f90 b/Prog_8/cgr2.f90 deleted file mode 100644 index 213a01b25..000000000 --- a/Prog_8/cgr2.f90 +++ /dev/null @@ -1,122 +0,0 @@ - SUBROUTINE CGR2(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ) - - ! B2 = U2*D2*V2 - ! B1 = V1*D1*U1 - !Calc: ( 1 B1 )^-1 i.e. 2*LQ \times 2*LQ matrix - ! (-B2 1 ) - - - Use Precdef - Use UDV_WRAP_mod - Use MyMats - - Implicit none - - ! Arguments - Integer :: LQ - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - - - ! Local:: - Complex (Kind=double) :: U3B(2*LQ,2*LQ), V3B(2*LQ,2*LQ), HLPB1(2*LQ,2*LQ), HLPB2(2*LQ,2*LQ), & - & V2INV(LQ,LQ), V1INV(LQ,LQ), HLP2(LQ,LQ) - Complex (Kind=double) :: D3B(2*LQ) - Complex (Kind=double) :: Z - - Integer :: LQ2, I,J, M, ILQ, JLQ, NCON, I1, J1 - - LQ2 = LQ*2 - - HLPB1 = cmplx(0.D0,0.d0,double) - DO I = 1,LQ - HLPB1(I , I + LQ ) = D1(I) - HLPB1(I+LQ, I ) = -D2(I) - ENDDO - CALL INV(V2,V2INV,Z) - CALL INV(V1,V1INV,Z) - CALL MMULT(HLP2,V1INV,V2INV) - DO J = 1,LQ - DO I = 1,LQ - HLPB1(I,J) = HLP2(I,J) - ENDDO - ENDDO - CALL MMULT(HLP2,U1,U2) - DO I = 1,LQ - ILQ = I+LQ - DO J = 1,LQ - JLQ = J + LQ - HLPB1(ILQ,JLQ) = conjg( HLP2(J,I) ) ! = (U1*U2)^T - ENDDO - ENDDO - NCON = 0 - CALL UDV_wrap(HLPB1,U3B,D3B,V3B,NCON) - - - ! Multiplication: - ! U3B^T * ( V1INV 0 ) = U3B - ! ( 0 U2^T ) - - DO I = 1,LQ2 - DO J = 1,LQ2 - HLPB1(I,J) = conjg(U3B(J,I)) - ENDDO - ENDDO - HLPB2 = cmplx(0.D0,0.d0,double) - DO I = 1,LQ - DO J = 1,LQ - HLPB2(I,J) = V1INV(I,J) - ENDDO - ENDDO - DO I = 1,LQ - ILQ = I + LQ - DO J = 1,LQ - JLQ = J + LQ - HLPB2(ILQ,JLQ) = conjg(U2(J,I)) - ENDDO - ENDDO - CALL MMULT(U3B,HLPB1,HLPB2) - - - ! Multiplication: - ! ( V2INV 0 )*(V3B)^{-1} = V3B - ! ( 0 U1^T ) - - CALL INV(V3B,HLPB1,Z) - HLPB2 = cmplx(0.d0,0.d0,double) - DO I = 1,LQ - DO J = 1,LQ - HLPB2(I,J) = V2INV(I,J) - ENDDO - ENDDO - DO I = 1,LQ - ILQ = I + LQ - DO J = 1, LQ - JLQ = J + LQ - HLPB2(ILQ,JLQ) = conjg(U1(J,I)) - ENDDO - ENDDO - CALL MMULT(V3B,HLPB2,HLPB1) - - - ! G = V3B * D3B^{-1}* U3B - DO M = 1,LQ2 - Z = cone/D3B(M) - DO J = 1,LQ2 - U3B(M,J) = Z * U3B(M,J) - ENDDO - ENDDO - CALL MMULT(HLPB2, V3B, U3B) - DO I = 1,LQ - I1 = I+LQ - DO J = 1,LQ - J1 = J + LQ - GR00(I,J) = HLPB2(I ,J ) - GRTT(I,J) = HLPB2(I1,J1) - GRT0(I,J) = HLPB2(I1,J ) - GR0T(I,J) = HLPB2(I,J1 ) - ENDDO - ENDDO - - END SUBROUTINE CGR2 diff --git a/Prog_8/cgr2_1.f90 b/Prog_8/cgr2_1.f90 deleted file mode 100644 index 78297bf8d..000000000 --- a/Prog_8/cgr2_1.f90 +++ /dev/null @@ -1,539 +0,0 @@ - SUBROUTINE CGR2_1(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ, NVAR) - - ! B2 = U2*D2*V2 is right (i.e. from time slice 0 to tau) propagation to time tau - ! B1 = V1*D1*U1 is left (i.e. from time slice Ltrot to tau) propagation to time tau - !Calc: ( 1 B1 )^-1 ( G00 G0T ) - ! (-B2 1 ) == ( GT0 GTT ) - ! - ! G00 = (1 + B1*B2)^-1 G0T = -(1 - G00 )*B2^-1 - ! GT0 = B2 * G00 GTT = (1 + B2*B1)^-1 - - ! Here you want to compute G00, G0T, GT0 and GTT just by involving LQ x LQ matrix operations. - ! If NVAR == 1 then the large scales are in D1 - ! If NVAR == 2 then the large scales are in D2 - Use Precdef - Use MyMats - USe UDV_Wrap_mod - - Implicit none - - Interface - SUBROUTINE CGR(Z,NVAR, GRUP, URUP,DRUP,VRUP, ULUP,DLUP,VLUP) - COMPLEX(Kind=8), Dimension(:,:), Intent(In) :: URUP, VRUP, ULUP, VLUP - COMPLEX(Kind=8), Dimension(:), Intent(In) :: DLUP, DRUP - COMPLEX(Kind=8), Dimension(:,:), Intent(INOUT) :: GRUP - - COMPLEX(Kind=8) :: Z - END SUBROUTINE CGR - end Interface - - - ! Arguments - Integer, intent(in) :: LQ, NVAR - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - - - ! Local:: - Complex (Kind=double) :: HLP1(LQ,LQ), HLP2(LQ,LQ), U(LQ,LQ), D(LQ), V(LQ,LQ) - Complex (Kind=double) :: Z, Z1, Z2 - Real (Kind=double) :: X, Xmax, Xmin, X1, X2, Xmax1, Xmax2, Xmean - Integer :: I, J, M, NCON, NVAR1 - - Complex (Kind=double) :: V2inv(LQ,LQ), V1inv(LQ,LQ) - - - NCON = 0 - - Call INV( V2, V2inv, Z2) - CALL INV( V1, V1inv, Z1) - - - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = CONJG( U1(J,I) ) - ENDDO - ENDDO - CALL MMULT(HLP2,HLP1,U1) - HLP1 = cmplx(0.d0,0.d0,kind=8) - DO I = 1,LQ - HLP1(I,I) = cmplx(1.d0,0.d0,kind=8) - ENDDO - Xmax = 0.d0 - CALL COMPARE(HLP1, HLP2, XMAX, XMEAN) - - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = CONJG( U2(J,I) ) - ENDDO - ENDDO - CALL MMULT(HLP2,HLP1,U2) - HLP1 = cmplx(0.d0,0.d0,kind=8) - DO I = 1,LQ - HLP1(I,I) = cmplx(1.d0,0.d0,kind=8) - ENDDO - Xmax1 = 0.d0 - CALL COMPARE(HLP1, HLP2, XMAX1, XMEAN) - Write(77,*) "Cgr2_1 V1inv V2inv : ", Xmax, Xmax1 - -!!$ Xmax = 0.d0 -!!$ do I = 1,LQ -!!$ do j = 1,LQ -!!$ X = sqrt(dble(V1(i,j)*conjg(V1(i,j)))) -!!$ if (X > Xmax) Xmax = X -!!$ enddo -!!$ enddo -!!$ Write(77,*) 'In cgr2_1 Xmax V1: ', Xmax, Z2 -!!$ do I = 1,LQ -!!$ do j = 1,LQ -!!$ X = sqrt(dble(V2(i,j)*conjg(V2(i,j)))) -!!$ if (X > Xmax) Xmax = X -!!$ enddo -!!$ enddo -!!$ Write(77,*) 'In cgr2_1 Xmax V2: ', Xmax, Z1 - - ! Compute G00 - ! G00 = (1 + B1*B2)^-1 = (1 + V1 D1 U1 U2 D2 V2 )^-1 = - ! = ( V1 ( V1^-1 V2^-1 + D1 U1 U2 D2 ) V2 )^-1 = - ! = V2^-1 ( (V2 V1)^-1 + D1 U1 U2 D2 )^-1 V1^-1 - Call MMULT(HLP1,V1inv,V2inv) - Call MMULT(HLP2,U1,U2) - DO J = 1,LQ - DO I = 1,LQ - HLP2(I,J) = D1(I)*HLP2(I,J)*D2(J) + HLP1(I,J) - ENDDO - ENDDO - Xmax1 = dble( D1(1) ) - Xmax2 = dble( D2(1) ) - DO I = 2,LQ - If ( dble( D1(I) ) > Xmax1 ) Xmax1 = dble( D1(I) ) - If ( dble( D2(I) ) > Xmax2 ) Xmax2 = dble( D2(I) ) - Enddo - Nvar1 = 1 - If ( Xmax2 > Xmax1) Nvar1 = 2 - If (Nvar1 == 1) then - ! V2^-1 (UDV )^-1 V1^-1 = V2^-1 V^-1 D^-1 U^-1 V1^-1 - Call UDV_WRAP(HLP2, U, D, V, Ncon) - CALL INV (V,HLP2 ,Z ) - CALL MMULT(V,V2inv,HLP2) - DO J = 1,LQ - DO I = 1,LQ - V(I,J) = V(I,J)/D(J) - ENDDO - ENDDO - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = Conjg(U(J,I)) - ENDDO - ENDDO - CALL MMULT( HLP2, HLP1,V1inv) - CALL MMULT (GR00, V, HLP2) - else - ! V2^-1 (UDV )^(-1,*) V1^-1 = V2^-1 U D^-1 V^(-1,*) V1^-1 - DO J = 1,LQ - DO I = 1,LQ - HLP1(I,J) = conjg(HLP2(J,I)) - ENDDO - ENDDO - Call UDV_WRAP(HLP1, U, D, V, Ncon) - Call MMULT(HLP1, V2inv, U) - DO J = 1,LQ - DO I = 1,LQ - HLP1(I,J) = HLP1(I,J)/D(J) - ENDDO - ENDDO - CALL INV (V, HLP2, Z) - DO J = 1,LQ - DO I = 1,LQ - V(I,J) = CONJG(HLP2(J,I)) - ENDDO - ENDDO - CALL MMULT(HLP2,V,V1inv) - CALL MMULT(GR00,HLP1,HLP2) - endif - - ! Compute G0T - ! G00 = (1 + B1*B2)^-1 G0T = -(1 - (1 + B1*B2)^-1 )*B2^-1 = - ! = -( 1 + B1*B2 - 1) (1 + B1*B2)^-1 * B2^-1 = - ! = - B1 * B2 (1 + B1*B2)^-1 * B2^-1 = - ! = -( B2 ( 1+ B1*B2) * B2^-1 B1^-1)^-1 = - ! = -( B2 ( 1+ B1*B2) * (B1 B2)^-1)^-1 = - ! = -( B2 ( (B1 B2)^-1 + 1 ) )^-1 = - ! = -( B1^-1 + B2)^-1 = - ! -G0T*= ( B1*^-1 + B2*)^-1 = - ! = ( V1*^-1 D1*^-1 U1 + V2* D2* U2*)^-1 = - ! = ( V1*^-1 ( D1*^-1 U1 U2 + V1* V2* D2* ) U2* )^-1 = - ! = U2 ( D1*^-1 (U1 U2) + ( V2 V1)* D2* )^-1 V1* - ! = U2 ( D1*^-1 (U1 U2) + ( V2 V1)* D2* )^-1 V1* - ! B2 = U2*D2*V2 - ! B1 = V1*D1*U1 - Xmax2 = dble(cmplx(1.d0,0.d0)/D1(1)) - Xmax1 = dble(D2(1)) - Do I = 2,LQ - X2 = dble(cmplx(1.d0,0.d0)/D1(I)) - X1 = dble(D2(I)) - If ( X2 > Xmax2 ) Xmax2 = X2 - If ( X1 > Xmax1 ) Xmax1 = X1 - ENDDO - NVAR1 = 1 - If (Xmax2 > Xmax1) Nvar1 = 2 - Call MMULT(HLP1,U1,U2) - DO J = 1,LQ - DO I =1,LQ - HLP1(I,J) = HLP1(I,J)/conjg(D1(I)) - ENDDO - ENDDO - Call MMULT(V,V2,V1) - DO J = 1,LQ - DO I = 1,LQ - HLP2(I,J) = Conjg(V(J,I)) - ENDDO - ENDDO - DO J = 1,LQ - DO I =1,LQ - HLP2(I,J) = HLP1(I,J) + HLP2(I,J)*conjg(D2(J)) - ENDDO - ENDDO - NCON = 0 - IF ( NVAR1 == 1 ) Then - ! UDV of HLP2 - ! -G0T*= U2 V^-1 D^-1 U* V1* - CALL UDV_WRAP(HLP2,U,D,V,NCON) - CALL MMULT (HLP1, V1, U) - DO I = 1,LQ - DO J = 1,LQ - U(I,J) = conjg(HLP1(J,I)) - ENDDO - ENDDO - CALL INV(V,HLP2,Z) - Call MMULT(HLP1,U2,HLP2) - DO J = 1,LQ - Z = cmplx(1.d0,0.d0,kind=8)/D(J) - DO I = 1,LQ - HLP1(I,J) = HLP1(I,J)*Z - ENDDO - ENDDO - Call MMULT (HLP2,HLP1,U) - DO I = 1,LQ - DO J = 1,LQ - GR0T(I,J) = -conjg(HLP2(J,I)) - ENDDO - ENDDO - ELSE - ! UDV of HLP2* - ! -G0T*= U2 (U D V)*^-1 V1* = U2 U D*^-1 V*^-1 V1* - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = conjg(HLP2(J,I)) - ENDDO - ENDDO - CALL UDV_WRAP(HLP1,U,D,V,NCON) - CALL MMULT (HLP1, U2, U) - DO J = 1,LQ - Z = cmplx(1.d0,0.d0,kind=8)/D(J) - DO I = 1,LQ - HLP1(I,J) = HLP1(I,J)*Z - ENDDO - ENDDO - CALL INV(V,HLP2,Z) - Call MMULT(V,V1,HLP2) - DO I = 1,LQ - DO J = 1,LQ - HLP2(I,J) = conjg(V(J,I)) - ENDDO - ENDDO - Call MMULT (V,HLP1,HLP2) - DO I = 1,LQ - DO J = 1,LQ - GR0T(I,J) = -conjg(V(J,I)) - ENDDO - ENDDO - ENDIF - - - - - ! Compute GT0 - ! GT0 = B2 * G00 = ( ( 1 + B1* B2) * B2^-1 )^-1 = ( B2^-1 + B1)^-1 = - ! = (V2^-1 D2^-1 U2^-1 + V1 D1 U1)^-1 = - ! = ( (V2^-1 D2^-1 U2^-1 U1^-1 + V1 D1 ) U1 )^-1 = - ! = U1^-1 ( ( D2^-1 (U1 U2)^-1 + V2*V1 D1 ) )^-1 V2 - Xmax2 = dble(cmplx(1.d0,0.d0)/D2(1)) - Xmax1 = dble(D1(1)) - Do I = 2,LQ - X2 = dble(cmplx(1.d0,0.d0)/D2(I)) - X1 = dble(D1(I)) - If ( X2 > Xmax2 ) Xmax2 = X2 - If ( X1 > Xmax1 ) Xmax1 = X1 - ENDDO - NVAR1 = 1 - If (Xmax2 > Xmax1 ) NVAR1 = 2 - !Write(6,*) "CGR2_1: NVAR,NVAR1 ", NVAR, NVAR1 - Call MMULT(HLP2,U1,U2) - DO J = 1,LQ - DO I = 1,LQ - HLP1(I,J) = Conjg(HLP2(J,I)) - ENDDO - ENDDO - DO J = 1,LQ - DO I =1,LQ - HLP1(I,J) = HLP1(I,J)/D2(I) - ENDDO - ENDDO - Call MMULT(HLP2,V2,V1) - DO J = 1,LQ - DO I =1,LQ - HLP2(I,J) = HLP1(I,J) + HLP2(I,J)*D1(J) - ENDDO - ENDDO - NCON = 0 - IF ( NVAR1 == 1 ) Then - ! UDV of HLP2 - CALL UDV_WRAP(HLP2,U,D,V,NCON) - CALL MMULT (HLP1, V, U1) - CALL INV(HLP1,HLP2,Z) - DO J = 1,LQ - Z = cmplx(1.d0,0.d0)/D(J) - DO I = 1,LQ - HLP2(I,J) = HLP2(I,J)*Z - ENDDO - ENDDO - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = Conjg(U(J,I)) - ENDDO - ENDDO - CALL MMULT(U,HLP1,V2) - Call MMULT (GRT0, HLP2,U) - ELSE - !UDV of HLP2^* - DO J = 1,LQ - DO I =1,LQ - HLP1(I,J) = Conjg(HLP2(J,I)) - ENDDO - ENDDO - CALL UDV_WRAP(HLP1,U,D,V,NCON) - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = conjg(U1(J,I)) - ENDDO - ENDDO - CALL MMULT( HLP2, HLP1,U) - DO J = 1,LQ - DO I = 1,LQ - HLP2(I,J) = HLP2(I,J)/Conjg(D(J)) - ENDDO - ENDDO - DO I = 1,LQ - DO J = 1,LQ - HLP1(I,J) = conjg(V(J,I)) - ENDDO - ENDDO - CALL INV(HLP1,V,Z) - CALL MMULT(U,V,V2) - Call MMULT (GRT0, HLP2,U) - ENDIF - Xmin = abs(dble(D(1))) - DO I = 1,LQ - if (abs(dble(D(I))) < Xmin ) Xmin = abs(dble(D(I))) - ENDDO - Write(6,*) 'Cgr2_1 T0, Xmin: ', Xmin - - - !Compute GRTT - Z = cmplx(1.d0,0.d0,kind=8) - Z1 = cmplx(1.d0,0.d0,kind=8) - CALL CGR(Z,NVAR,GRTT, U2,D2,V2, U1,D1,V1) - - - END SUBROUTINE CGR2_1 - - -!!$ ! Compute G0T -!!$ ! G00 = (1 + B1*B2)^-1 G0T = -(1 - (1 + B1*B2)^-1 )*B2^-1 = -!!$ ! = -( 1 + B1*B2 - 1) (1 + B1*B2)^-1 * B2^-1 = -!!$ ! = - B1 * B2 (1 + B1*B2)^-1 * B2^-1 = -!!$ ! = -( B2 ( 1+ B1*B2) * B2^-1 B1^-1)^-1 = -!!$ ! = -( B2 ( 1+ B1*B2) * (B1 B2)^-1)^-1 = -!!$ ! = -( B2 ( (B1 B2)^-1 + 1 ) )^-1 = -!!$ ! = -( B1^-1 + B2)^-1 = -!!$ ! = -( U1^-1 D1^-1 V1^-1 + U2 D2 V2)^-1 = -!!$ ! = -( ( U1^-1 D1^-1 V1^-1 V2^-1 + U2 D2 ) V2 )^-1 = -!!$ ! = -( U1^-1( D1^-1 (V2 V1)^-1 + U1 U2 D2) V2 )^-1 = -!!$ ! = - V2^-1( D1^-1 (V2 V1)^-1 + U1 U2 D2)^-1 U1 -!!$ ! B2 = U2*D2*V2 -!!$ ! B1 = V1*D1*U1 -!!$ Call MMULT (HLP2, V1inv,V2inv) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = HLP2(I,J)/D1(I) -!!$ ENDDO -!!$ ENDDO -!!$ Call MMULT (HLP1, U1,U2) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = HLP2(I,J) + HLP1(I,J)*D2(J) -!!$ ENDDO -!!$ ENDDO -!!$ Xmax2 = dble(cmplx(1.d0,0.d0,Kind=8)/D1(1)) -!!$ Xmax1 = dble(D2(1)) -!!$ Do I = 2,LQ -!!$ X2 = dble(cmplx(1.d0,0.d0,Kind=8)/D1(I)) -!!$ X1 = dble(D2(I)) -!!$ If ( X2 > Xmax2 ) Xmax2 = X2 -!!$ If ( X1 > Xmax1 ) Xmax1 = X1 -!!$ ENDDO -!!$ NVAR1 = 1 -!!$ If (Xmax2 > Xmax1 ) NVAR1 = 2 -!!$ IF (NVAR1 == 1) Then -!!$ ! UDV of HLP2 -!!$ != - V2^-1( U D V )^-1 ) U1 = -!!$ != - V2^-1 V^-1 D^-1 U^-1 U1 = - (V V2)^-1 D^-1 U^-1 U1 -!!$ CALL UDV(HLP2,U,D,V,NCON) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP1(I,J) = conjg(U(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL MMULT( U, HLP1, U1 ) -!!$ CALL MMULT(HLP1,V, V2) -!!$ CALL INV (HLP1,V ,Z) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP1(I,J) = - V(I,J)/D(J) -!!$ ENDDO -!!$ ENDDO -!!$ CALL MMULT(GR0T, HLP1, U) -!!$ Else -!!$ ! UDV of HLP2^* -!!$ != - V2^-1( U D V)^*,-1 ) U1 = -!!$ != - V2^-1 U D^-1 V^*,-1 U1 -!!$ DO I = 1,LQ -!!$ DO J = 1,LQ -!!$ HLP1(J,I) = Conjg(HLP2(I,J)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL UDV(HLP1,U,D,V,NCON) -!!$ CALL INV(V2,HLP1,Z) -!!$ CALL MMULT(HLP2,HLP1,U) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP1(I,J) = -HLP2(I,J)/conjg(D(J)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL INV(V,HLP2,Z) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ V(I,J) = Conjg(HLP2(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL MMULT(HLP2,V,U1) -!!$ CALL MMULT(GR0T, HLP1,HLP2) -!!$ endif -!!$ Xmin = abs(dble(D(1))) -!!$ DO I = 1,LQ -!!$ if (abs(dble(D(I))) < Xmin ) Xmin = abs(dble(D(I))) -!!$ ENDDO -!!$ Write(6,*) 'Cgr2_1 0T, Xmin: ', Xmin - - - - -!!$ ! Compute G0T -!!$ ! G00 = (1 + B1*B2)^-1 G0T = -(1 - (1 + B1*B2)^-1 )*B2^-1 = -!!$ ! = -( 1 + B1*B2 - 1) (1 + B1*B2)^-1 * B2^-1 = -!!$ ! = - B1 * B2 (1 + B1*B2)^-1 * B2^-1 = -!!$ ! = -( B2 ( 1+ B1*B2) * B2^-1 B1^-1)^-1 = -!!$ ! = -( B2 ( 1+ B1*B2) * (B1 B2)^-1)^-1 = -!!$ ! = -( B2 ( (B1 B2)^-1 + 1 ) )^-1 = -!!$ ! = -( B1^-1 + B2)^-1 = -!!$ ! = -( U1^-1 D1^-1 V1^-1 + U2 D2 V2)^-1 = -!!$ ! = -(U2 (U2^-1 U1^-1 D1^-1 + D2 V2 V1 ) V1^-1 )^-1 = -!!$ ! = - V1 ( (U1 U2)^-1 D1^-1 + D2 V2 V1 )^-1 U2^-1 -!!$ ! B2 = U2*D2*V2 -!!$ ! B1 = V1*D1*U1 -!!$ Call MMULT (HLP1, U1,U2) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = Conjg(HLP1(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = HLP2(I,J) / D1(J) -!!$ ENDDO -!!$ ENDDO -!!$ -!!$ Call MMULT (HLP1, V2,V1) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = HLP2(I,J) + D2(I)*HLP1(I,J) -!!$ ENDDO -!!$ ENDDO -!!$ Xmax2 = dble(cmplx(1.d0,0.d0)/D1(1)) -!!$ Xmax1 = dble(D2(1)) -!!$ Do I = 2,LQ -!!$ X2 = dble(cmplx(1.d0,0.d0)/D1(I)) -!!$ X1 = dble(D2(I)) -!!$ If ( X2 > Xmax2 ) Xmax2 = X2 -!!$ If ( X1 > Xmax1 ) Xmax1 = X1 -!!$ ENDDO -!!$ NVAR1 = 1 -!!$ If (Xmax1 > Xmax2 ) NVAR1 = 2 -!!$ IF (NVAR1 == 1) Then -!!$ ! UDV of HLP2 -!!$ != - V1 ( U D V)^-1 U2^-1 -!!$ != - V1 V^-1 D^-1 U^-1 U2^-1 = - V1 V^-1 D^-1 (U2 U)^-1 -!!$ CALL UDV(HLP2,U,D,V,NCON) -!!$ CALL MMULT( HLP2, U2, U ) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP1(I,J) = conjg(HLP2(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL INV (V, HLP2 ,Z) -!!$ CALL MMULT(V, V1, HLP2) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = - V(I,J)/D(J) -!!$ ENDDO -!!$ ENDDO -!!$ CALL MMULT(GR0T, HLP2, HLP1) -!!$ Else -!!$ ! UDV of HLP2^* -!!$ != - V1 ( U D V)^(*,-1) U2^-1 -!!$ != - V1 U D^(*,-1) V^(*,-1) U2^-1 -!!$ DO I = 1,LQ -!!$ DO J = 1,LQ -!!$ HLP1(J,I) = Conjg(HLP2(I,J)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL UDV(HLP1,U,D,V,NCON) -!!$ CALL MMULT(HLP2,V1,U) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ HLP2(I,J) = -HLP2(I,J)/conjg(D(J)) -!!$ ENDDO -!!$ ENDDO -!!$ -!!$ CALL INV(V,HLP1,Z) -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ V(I,J) = Conjg(HLP1(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ DO J = 1,LQ -!!$ DO I = 1,LQ -!!$ U(I,J) = Conjg(U2(J,I)) -!!$ ENDDO -!!$ ENDDO -!!$ CALL MMULT(HLP1,V,U) -!!$ -!!$ CALL MMULT(GR0T, HLP2,HLP1) -!!$ endif -!!$ Xmin = abs(dble(D(1))) -!!$ DO I = 1,LQ -!!$ if (abs(dble(D(I))) < Xmin ) Xmin = abs(dble(D(I))) -!!$ ENDDO -!!$ Write(6,*) 'Cgr2_1 0T, Xmin: ', Xmin, NVAR1 diff --git a/Prog_8/cgr2_2.f90 b/Prog_8/cgr2_2.f90 deleted file mode 100644 index 87e8a462f..000000000 --- a/Prog_8/cgr2_2.f90 +++ /dev/null @@ -1,176 +0,0 @@ - SUBROUTINE CGR2_2(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ) - - - ! B2 = U2*D2*V2 is right (i.e. from time slice 0 to tau) propagation to time tau - ! B1 = V1*D1*U1 is left (i.e. from time slice Ltrot to tau) propagation to time tau - !Calc: ( 1 B1 )^-1 ( G00 G0T ) - ! (-B2 1 ) == ( GT0 GTT ) - ! - ! G00 = (1 + B1*B2)^-1 G0T = -(1 - G00)*B2^-1 - ! GT0 = B2 * G00 GTT = (1 + B2*B1)^-1 - - !( 1 V1*D1*U1 )^-1 ( ( V1 0 ) ( V1^-1 D1*U1 ) )^-1 - !(-U2*D2*V2 1 ) == ( ( 0 U2 ) * (-D2*V2 U2^-1 ) ) == I - ! You should transpose before carrying out the singular value decomposition - ! - ! - ! ( ( V1 0 ) ( V1^-1 D1*U1 )^*^* )^-1 ( V1^-1 0 ) - ! I == ( ( 0 U2 ) * (-D2*V2 U2^-1 ) ) = (UDV^*)^(-1) * ( 0 U2^-1) = - ! - ! ( V1^-1 0 ) - ! == U * D^(*,-1) * V^(*,-1) * ( 0 U2^-1) - - ! Let's see if this could work. - Use Precdef - Use MyMats - Use UDV_WRAP_mod - Implicit none - - ! Arguments - Integer, intent(in) :: LQ - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - - - ! Local:: - Complex (Kind=double) :: U3B(2*LQ,2*LQ), V3B(2*LQ,2*LQ), HLPB1(2*LQ,2*LQ), HLPB2(2*LQ,2*LQ), & - & V2INV(LQ,LQ), V1INV(LQ,LQ), HLP2(LQ,LQ) - Complex (Kind=double) :: D3B(2*LQ) - Complex (Kind=double) :: Z - Real (Kind=double) :: X, Xmax - - Integer :: LQ2, I,J, M, ILQ, JLQ, NCON, I1, J1,N - - LQ2 = LQ*2 - NCON = 0 - - If (dble(D1(1)) > dble(D2(1)) ) Then - - !Write(6,*) "D1(1) > D2(1)", dble(D1(1)), dble(D2(1)) - - HLPB2 = cmplx(0.D0,0.d0,double) - CALL INV(V1,V1INV,Z) - DO J = 1,LQ - DO I = 1,LQ - HLPB2(I , J ) = V1INV(I,J) - HLPB2(I , J+LQ ) = D1(I)*U1(I,J) - HLPB2(I+LQ, J+LQ ) = Conjg(U2(J,I)) - HLPB2(I+LQ, J ) = -D2(I)*V2(I,J) - ENDDO - ENDDO - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB1(I,J) = Conjg(HLPB2(J,I)) - ENDDO - ENDDO - - !CALL UDV_wrap(HLPB1,U3B,D3B,V3B,NCON) - CALL UDV_wrap_Pivot(HLPB1,U3B,D3B,V3B,NCON,LQ2,LQ2) - -!!$!!!!!!!!!!!!! Tests -!!$ Xmax = 0.d0 -!!$ DO I = 1,LQ2 -!!$ DO J = 1,LQ2 -!!$ Z = cmplx(0.d0,0.d0) -!!$ DO N = 1,LQ2 -!!$ Z = Z + U3B(I,N) *conjg(U3B(J,N)) -!!$ ENDDO -!!$ if (I == J) Z = Z - cmplx(1.d0,0.d0) -!!$ X = real(SQRT( Z* conjg(Z)),kind=8) -!!$ if (X > Xmax) Xmax = X -!!$ ENDDO -!!$ ENDDO -!!$ !Write(6,*) 'Cgr2_2, ortho: ', Xmax -!!$ DO I = 1,LQ2 -!!$ Z = D3B(I) -!!$ if (I == 1) Xmax = real(SQRT( Z* conjg(Z)),kind=8) -!!$ if ( real(SQRT( Z* conjg(Z)),kind=8) < Xmax ) Xmax = & -!!$ & real(SQRT( Z* conjg(Z)),kind=8) -!!$ ENDDO -!!$ !Write(6,*) 'Cgr2_2, Cutoff: ', Xmax -!!$!!!!!!!!!!!!! End Tests - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB2(I,J) = Conjg(V3B(J,I)) - ENDDO - ENDDO - CALL INV(HLPB2,V3B,Z) - HLPB1 = cmplx(0.d0,0.d0,double) - DO I = 1,LQ - DO J = 1,LQ - HLPB1(I , J ) = V1INV(I,J) - HLPB1(I+LQ, J+LQ ) = Conjg(U2(J,I)) - ENDDO - ENDDO - CALL MMULT(HLPB2,V3B,HLPB1) - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB1(I,J) = Conjg(cmplx(1.d0,0.d0,double)/D3B(I))*HLPB2(I,J) - ENDDO - ENDDO - CALL MMULT(HLPB2,U3B,HLPB1) - DO I = 1,LQ - I1 = I+LQ - DO J = 1,LQ - J1 = J + LQ - GR00(I,J) = HLPB2(I ,J ) - GRTT(I,J) = HLPB2(I1,J1) - GRT0(I,J) = HLPB2(I1,J ) - GR0T(I,J) = HLPB2(I,J1 ) - ENDDO - ENDDO - Else - !Write(6,*) "D1(1) < D2(1)", dble(D1(1)), dble(D2(1)) - HLPB2 = cmplx(0.D0,0.d0,double) - CALL INV(V1,V1INV,Z) - DO J = 1,LQ - DO I = 1,LQ - HLPB2(I , J ) = Conjg(U2(J,I)) - HLPB2(I , J+LQ ) = -D2(I)*V2(I,J) - HLPB2(I+LQ, J+LQ ) = V1INV(I,J) - HLPB2(I+LQ, J ) = D1(I)*U1(I,J) - ENDDO - ENDDO - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB1(I,J) = Conjg(HLPB2(J,I)) - ENDDO - ENDDO - - !CALL UDV_wrap(HLPB1,U3B,D3B,V3B,NCON) - CALL UDV_wrap_Pivot(HLPB1,U3B,D3B,V3B,NCON,LQ2,LQ2) - - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB2(I,J) = Conjg(V3B(J,I)) - ENDDO - ENDDO - CALL INV(HLPB2,V3B,Z) - HLPB1 = cmplx(0.d0,0.d0,double) - DO I = 1,LQ - DO J = 1,LQ - HLPB1(I , J ) = Conjg(U2(J,I)) - HLPB1(I+LQ, J+LQ ) = V1INV(I,J) - ENDDO - ENDDO - CALL MMULT(HLPB2,V3B,HLPB1) - DO J = 1,LQ2 - DO I = 1,LQ2 - HLPB1(I,J) = Conjg(cmplx(1.d0,0.d0,double)/D3B(I))*HLPB2(I,J) - ENDDO - ENDDO - CALL MMULT(HLPB2,U3B,HLPB1) - DO I = 1,LQ - I1 = I+LQ - DO J = 1,LQ - J1 = J + LQ - GRTT(I,J) = HLPB2(I ,J ) - GR00(I,J) = HLPB2(I1,J1) - GR0T(I,J) = HLPB2(I1,J ) - GRT0(I,J) = HLPB2(I,J1 ) - ENDDO - ENDDO - Endif - - END SUBROUTINE CGR2_2 diff --git a/Prog_8/control_mod.f90 b/Prog_8/control_mod.f90 deleted file mode 100644 index ba1b18a04..000000000 --- a/Prog_8/control_mod.f90 +++ /dev/null @@ -1,142 +0,0 @@ - module Control - - Use MyMats - Implicit none - - real (Kind=8) , private, save :: XMEANG, XMAXG, XMAXP, CPU_time_st, CPU_time_en, Xmean_tau, Xmax_tau - Integer , private, save :: NCG, NCG_tau - Integer (Kind=8), private, save :: NC_up, ACC_up - - Contains - - subroutine control_init - Implicit none - XMEANG = 0.d0 - XMEAN_tau = 0.d0 - XMAXG = 0.d0 - XMAX_tau = 0.d0 - NCG = 0 - NCG_tau = 0 - NC_up = 0 - ACC_up = 0 - Call CPU_TIME(CPU_time_st) - end subroutine control_init - - Subroutine Control_upgrade(Log) - Implicit none - Logical :: Log - NC_up = NC_up + 1 - if (Log) ACC_up = ACC_up + 1 - end Subroutine Control_upgrade - - Subroutine Control_PrecisionG(A,B,Ndim) - Implicit none - - Integer :: Ndim - Complex (Kind=8) :: A(Ndim,Ndim), B(Ndim,Ndim) - Real (Kind=8) :: XMAX, XMEAN - - !Local - NCG = NCG + 1 - XMEAN = 0.d0 - XMAX = 0.d0 - CALL COMPARE(A, B, XMAX, XMEAN) - IF (XMAX > XMAXG) XMAXG = XMAX - XMEANG = XMEANG + XMEAN - !Write(6,*) 'Control', XMEAN, XMAX - End Subroutine Control_PrecisionG - - Subroutine Control_Precision_tau(A,B,Ndim) - Implicit none - - Integer :: Ndim - Complex (Kind=8) :: A(Ndim,Ndim), B(Ndim,Ndim) - Real (Kind=8) :: XMAX, XMEAN - - !Local - NCG_tau = NCG_tau + 1 - XMEAN = 0.d0 - XMAX = 0.d0 - CALL COMPARE(A, B, XMAX, XMEAN) - IF (XMAX > XMAX_tau) XMAX_tau = XMAX - XMEAN_tau = XMEAN_tau + XMEAN - !Write(6,*) 'Control_tau', XMEAN, XMAX - End Subroutine Control_Precision_tau - - - Subroutine Control_PrecisionP(Z,Z1) - Implicit none - Complex (Kind=8), INTENT(IN) :: Z,Z1 - Real (Kind=8) :: X - X = sqrt(dble((Z-Z1)*conjg(Z-Z1))) - if ( X > XMAXP ) XMAXP = X - End Subroutine Control_PrecisionP - - - Subroutine control_Print - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - Real (Kind=8) :: Time, Acc -#ifdef MPI - REAL (KIND=8) :: X - Integer :: Ierr, Isize, Irank - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - ACC = 0.d0 - IF (NC_up > 0 ) ACC = dble(ACC_up)/dble(NC_up) - Call CPU_TIME(CPU_time_en) - Time = CPU_time_en - CPU_time_st -#ifdef MPI - X = 0.d0 - CALL MPI_REDUCE(XMEANG,X,1,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - XMEANG = X/dble(Isize) - X = 0.d0 - CALL MPI_REDUCE(XMEAN_tau,X,1,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - XMEAN_tau = X/dble(Isize) - X = 0.d0 - CALL MPI_REDUCE(ACC,X,1,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - ACC = X/dble(Isize) - - X = 0.d0 - CALL MPI_REDUCE(Time,X,1,MPI_REAL8,MPI_SUM, 0,MPI_COMM_WORLD,IERR) - Time = X/dble(Isize) - - - CALL MPI_REDUCE(XMAXG,X,1,MPI_REAL8,MPI_MAX, 0,MPI_COMM_WORLD,IERR) - XMAXG = X - CALL MPI_REDUCE(XMAX_tau,X,1,MPI_REAL8,MPI_MAX, 0,MPI_COMM_WORLD,IERR) - XMAX_tau= X - - - CALL MPI_REDUCE(XMAXP,X,1,MPI_REAL8,MPI_MAX, 0,MPI_COMM_WORLD,IERR) - XMAXP = X - If (Irank == 0 ) then -#endif - - Open (Unit=50,file="info", status="unknown", position="append") - If (NCG > 0 ) then - XMEANG = XMEANG/dble(NCG) - Write(50,*) ' Precision Green Mean, Max : ', XMEANG, XMAXG - Write(50,*) ' Precision Phase, Max : ', XMAXP - endif - If ( NCG_tau > 0 ) then - XMEAN_tau = XMEAN_tau/dble(NCG_tau) - Write(50,*) ' Precision tau Mean, Max : ', XMEAN_tau, XMAX_tau - endif - Write(50,*) ' Acceptance : ', ACC - Write(50,*) ' CPU Time : ', Time - Close(50) -#ifdef MPI - endif -#endif - end Subroutine Control_Print - - end module control - - diff --git a/Prog_8/gperp.f90 b/Prog_8/gperp.f90 deleted file mode 100644 index acb61d3c2..000000000 --- a/Prog_8/gperp.f90 +++ /dev/null @@ -1,98 +0,0 @@ - Subroutine Gperp_sub( G, Gperp, Ndim,Irank) - - Use Precdef - Use MyMats - Implicit none - - ! Arguments - Integer, Intent(In) :: Ndim, Irank - Complex (kind=double), Intent(In) :: G(ndim,ndim) - Complex (kind=double), Intent(InOut) :: Gperp(ndim,ndim) - - ! Local space - Complex (Kind=double) :: A(ndim,ndim), W(ndim), VL(Ndim,ndim), VR(Ndim,ndim) - Character (len=1) :: JOBVL, JOBVR - Integer :: INFO, LDA, LDVL, LDVR, N, lp, LWORK, N_c,m, i, j, NCon - Complex (Kind=double) :: WORK(2*Ndim), U(Ndim,Ndim/2), Vec(Ndim),Z - Real (Kind=double) :: RWORK(2*ndim), X, Xmax, Xmean - Complex (Kind=double) :: U1(Ndim,Ndim/2), V(Ndim/2,Ndim/2), D(Ndim/2) - - - A = G - JOBVL = "N" - JOBVR = "V" - LDA = Ndim - LWORK = 2*Ndim - LDVL = Ndim - LDVR = Ndim - N = Ndim - - Call ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) - - !lp = 70 + Irank - !Write(lp,*) "Info: ", INFO - N_c = 0 - do n = 1,Ndim - !Write(lp,*) n, W(n) - if ( abs( dble(W(n)) ) < 0.00001 ) then - N_c = N_c + 1 - do i = 1,Ndim - U1(i,N_c) = VR(i,n) - enddo - endif - enddo - !Write(6,*) "N_c ", N_c - NCON = 0 - Call UDV (U1,U,D,V,NCON) - - ! Setpup G_perp - gperp = cmplx(0.d0,0.d0) - Do i = 1,Ndim - do j = 1,Ndim - do n = 1,Ndim/2 - Gperp(i,j) = Gperp(i,j) + U(i,n) * conjg( U(j,n) ) - enddo - enddo - enddo - -#ifdef Test_gperp - X = 0.05 - A = cmplx(1.d0-X,0.d0) * G + cmplx(x,0.d0)*Gperp - Call Inv(A,VR,Z) - Write(lp,*) "Det is ", Z - Call MMult(VL,A,VR) - VR = cmplx(0.d0,0.d0) - do i = 1,Ndim - VR(I,I) = cmplx(1.d0,0.d0) - enddo - Call Compare(VL,VR,Xmax,Xmean) - Write(lp,*) 'Compare: ', Xmax, Xmean - - ! This is for testing - do n = 1,N_c - Vec = cmplx(0.d0,0.d0) - do i = 1,Ndim - do j = 1,Ndim - Vec(i) = Vec(i) + G(i,j) * U(j,n) - enddo - enddo - X = 0.d0 - do i = 1,Ndim - X = X + dble( Vec(i) * conjg(Vec(i))) - enddo - X = sqrt(x) - Write(lp,*) 'n, G*v = ', n, X - enddo - - do n = 1,N_c - do m = n,N_c - Z = cmplx(0.d0,0.d0) - do j = 1,Ndim - Z = Z + Conjg(U(j,m)) * U(j,n) - enddo - Write(lp,*) "n,m,z ", n,m,z - enddo - enddo -#endif - - end Subroutine Gperp_sub diff --git a/Prog_8/inconfc.f90 b/Prog_8/inconfc.f90 deleted file mode 100644 index 97dffc123..000000000 --- a/Prog_8/inconfc.f90 +++ /dev/null @@ -1,126 +0,0 @@ - SUBROUTINE confin - - Use Hamiltonian - - Implicit none - -#include "machine" - - - -#ifdef MPI - INCLUDE 'mpif.h' - ! Local -#endif - - Integer :: I, IERR, ISIZE, IRANK, seed_in, K, iseed, Nt - Integer, dimension(:), allocatable :: Seed_vec - Real (Kind=8) :: X - Logical :: lconf - character (len=64) :: file_sr, File_tg - -#ifdef MPI - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - Allocate (Nsigma(Size(Op_V,1),Ltrot)) - -#ifdef MPI - INQUIRE (FILE='confin_0', EXIST=lconf) - If (lconf) Then - file_sr = "confin" - Call Get_seed_Len(K) - Allocate(Seed_vec(K)) - file_tg = File_i(file_sr,IRANK) - Open (Unit = 10, File=File_tg, status='old', ACTION='read') - Read(10,*) Seed_vec - Call Ranset(Seed_vec) - do NT = 1,LTROT - do I = 1,Size(Op_V,1) - Read(10,*) NSIGMA(I,NT) - enddo - enddo - close(10) - Deallocate(Seed_vec) - else - If (Irank == 0) then - Write(6,*) 'No initial configuration' - OPEN(UNIT=5,FILE='seeds',STATUS='old',ACTION='read',IOSTAT=ierr) - IF (ierr /= 0) THEN - WRITE(*,*) 'unable to open ',ierr - STOP - END IF - DO I = Isize-1,1,-1 - Read (5,*) Seed_in - CALL MPI_SEND(Seed_in,1,MPI_INTEGER, I, I+1024,MPI_COMM_WORLD,IERR) - enddo - Read(5,*) Seed_in - CLOSE(5) - else - CALL MPI_RECV(Seed_in, 1, MPI_INTEGER,0, IRANK + 1024, MPI_COMM_WORLD,STATUS,IERR) - endif - Call Get_seed_Len(K) - !Write(6,*) K - Allocate(Seed_vec(K)) - Do I = 1,K - X = Ranf_Imada(Seed_in) - Seed_vec(I) = Seed_in - enddo - Call Ranset(Seed_vec) - Deallocate(Seed_vec) - do NT = 1,LTROT - do I = 1,Size(Op_V,1) - X = RANF() - NSIGMA(I,NT) = 1 - IF (X.GT.0.5) NSIGMA(I,NT) = -1 - enddo - enddo - endif - -#else - INQUIRE (FILE='confin_0', EXIST=lconf) - If (lconf) Then - file_tg = "confin_0" - Call Get_seed_Len(K) - Allocate(Seed_vec(K)) - Open (Unit = 10, File=File_tg, status='old', ACTION='read') - Read(10,*) Seed_vec - Call Ranset(Seed_vec) - do NT = 1,LTROT - do I = 1,Size(Op_V,1) - Read(10,*) NSIGMA(I,NT) - enddo - enddo - close(10) - Deallocate(Seed_vec) - else - Write(6,*) 'No initial configuration' - OPEN(UNIT=5,FILE='seeds',STATUS='old',ACTION='read',IOSTAT=ierr) - IF (ierr /= 0) THEN - WRITE(*,*) 'unable to open ',ierr - STOP - END IF - Read (5,*) Seed_in - CLOSE(5) - Call Get_seed_Len(K) - !Write(6,*) K - Allocate(Seed_vec(K)) - Do I = 1,K - X = Ranf_Imada(Seed_in) - Seed_vec(I) = Seed_in - enddo - Call Ranset(Seed_vec) - Deallocate(Seed_vec) - do NT = 1,LTROT - do I = 1,Size(Op_V,1) - X = RANF() - NSIGMA(I,NT) = 1 - IF (X.GT.0.5) NSIGMA(I,NT) = -1 - enddo - enddo - endif -#endif - - END SUBROUTINE CONFIN diff --git a/Prog_8/machine b/Prog_8/machine deleted file mode 100644 index 2e1fc39e2..000000000 --- a/Prog_8/machine +++ /dev/null @@ -1 +0,0 @@ -#define noMPI diff --git a/Prog_8/main.f90 b/Prog_8/main.f90 deleted file mode 100644 index bd810ccd1..000000000 --- a/Prog_8/main.f90 +++ /dev/null @@ -1,449 +0,0 @@ -Program Main - - Use Operator_mod - Use Lattices_v3 - Use MyMats - Use Hamiltonian - Use Control - Use Tau_m_mod - Use Hop_mod - Implicit none -#include "machine" -#ifdef MPI - include 'mpif.h' -#endif - - - Interface - SUBROUTINE WRAPUL(NTAU1, NTAU, UL ,DL, VL) - Use Hamiltonian - Implicit none - COMPLEX (KIND=8) :: UL(Ndim,Ndim,N_FL), VL(Ndim,Ndim,N_FL) - COMPLEX (KIND=8) :: DL(Ndim,N_FL) - Integer :: NTAU1, NTAU - END SUBROUTINE WRAPUL - SUBROUTINE CGR(PHASE,NVAR, GRUP, URUP,DRUP,VRUP, ULUP,DLUP,VLUP) - Use UDV_Wrap_mod - Implicit None - COMPLEX(Kind=8), Dimension(:,:), Intent(In) :: URUP, VRUP, ULUP, VLUP - COMPLEX(Kind=8), Dimension(:), Intent(In) :: DLUP, DRUP - COMPLEX(Kind=8), Dimension(:,:), Intent(Inout) :: GRUP - COMPLEX(Kind=8) :: PHASE - INTEGER :: NVAR - END SUBROUTINE CGR - SUBROUTINE WRAPGRUP(GR,NTAU,PHASE) - Use Hamiltonian - Implicit none - COMPLEX (Kind=8), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) - COMPLEX (Kind=8), INTENT(INOUT) :: PHASE - INTEGER, INTENT(IN) :: NTAU - END SUBROUTINE WRAPGRUP - SUBROUTINE WRAPGRDO(GR,NTAU,PHASE) - Use Hamiltonian - Implicit None - COMPLEX (Kind=8), INTENT(INOUT) :: GR(NDIM,NDIM,N_FL) - COMPLEX (Kind=8), INTENT(INOUT) :: PHASE - Integer :: NTAU - end SUBROUTINE WRAPGRDO - SUBROUTINE WRAPUR(NTAU, NTAU1, UR, DR, VR) - Use Hamiltonian - Use UDV_Wrap_mod - Implicit None - COMPLEX (KIND=8) :: UR(Ndim,Ndim,N_FL), VR(Ndim,Ndim,N_FL) - COMPLEX (KIND=8) :: DR(Ndim,N_FL) - Integer :: NTAU1, NTAU - END SUBROUTINE WRAPUR - - end Interface - - COMPLEX (Kind=8), Dimension(:) , Allocatable :: D - COMPLEX (KIND=8), Dimension(:,:) , Allocatable :: TEST, A, U, V - - COMPLEX (Kind=8), Dimension(:,:) , Allocatable :: DL, DR - COMPLEX (Kind=8), Dimension(:,:,:), Allocatable :: UL, VL, UR, VR - COMPLEX (Kind=8), Dimension(:,:,:), Allocatable :: GR - - - Integer :: Nwrap, NSweep, NBin, Ltau, NSTM, NT, NT1, NVAR, LOBS_EN, LOBS_ST, NBC, NSW - Integer :: NTAU, NTAU1 - - NAMELIST /VAR_QMC/ Nwrap, NSweep, NBin, Ltau, LOBS_EN, LOBS_ST - - Integer :: Ierr, I,J,nf, nst, n - Complex (Kind=8) :: Z_ONE = cmplx(1.d0,0.d0), Phase, Z, Z1 - - ! Space for storage. - COMPLEX (Kind=8), Dimension(:,:,:) , Allocatable :: DST - COMPLEX (Kind=8), Dimension(:,:,:,:), Allocatable :: UST, VST - - ! For tests - Integer, external :: nranf - Real (kind=8) :: Weight - Integer :: nr,nth - Logical :: Log -#ifdef MPI - Integer :: Isize, Irank - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_INIT(ierr) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) -#endif - - ! Write(6,*) 'Call Ham_set' - Call Ham_set - ! Write(6,*) 'End Call Ham_set' - Call confin - Call Hop_mod_init - !Call Hop_mod_test - !stop - -#ifdef MPI - If ( Irank == 0 ) then -#endif - OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) - IF (ierr /= 0) THEN - WRITE(*,*) 'unable to open ',ierr - STOP - END IF - READ(5,NML=VAR_QMC) - CLOSE(5) -#ifdef MPI - Endif - CALL MPI_BCAST(Nwrap ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(NSweep ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(NBin ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(Ltau ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(LOBS_EN ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - CALL MPI_BCAST(LOBS_ST ,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) -#endif - - - Call control_init - Call Alloc_obs(Ltau) - Call Op_SetHS - - -!!$#ifdef Ising_test -!!$ ! Test Ising -!!$ DO NBC = 1, NBIN -!!$ Call Init_obs -!!$ DO NSW = 1, NSWEEP -!!$ do nth = 1,Ltrot*2*Latt%N -!!$ Nt = nranf(Ltrot) -!!$ Nr = nranf(2*Latt%N) -!!$ Weight = S0(nr,nt) -!!$ log =.false. -!!$ if (Weight > ranf()) then -!!$ nsigma(nr,nt) = - nsigma(nr,nt) -!!$ log =.true. -!!$ endif -!!$ Call Control_upgrade(log) -!!$ enddo -!!$ Call Obser -!!$ Enddo -!!$ Call Preq -!!$ Enddo -!!$ Call Ham_confout -!!$ Call control_Print -!!$ Stop -!!$ ! End Test Ising -!!$#endif - - Allocate( DL(NDIM,N_FL), DR(NDIM,N_FL) ) - Allocate( UL(NDIM,NDIM,N_FL), VL(NDIM,NDIM,N_FL), & - & UR(NDIM,NDIM,N_FL), VR(NDIM,NDIM,N_FL), GR(NDIM,NDIM,N_FL ) ) - NSTM = LTROT/NWRAP -#ifdef MPI - if ( Irank == 0 ) then -#endif - Open (Unit = 50,file="info",status="unknown",position="append") - Write(50,*) 'Sweeps : ', Nsweep - Write(50,*) 'Bin : ', NBin - Write(50,*) 'Measure Int. : ', LOBS_ST, LOBS_EN - Write(50,*) 'Stabilization,Wrap : ', Nwrap - Write(50,*) 'Nstm : ', NSTM - Write(50,*) 'Ltau : ', Ltau - close(50) -#ifdef MPI - endif -#endif - - Allocate ( UST(NDIM,NDIM,NSTM,N_FL), VST(NDIM,NDIM,NSTM,N_FL), DST(NDIM,NSTM,N_FL) ) - Allocate ( Test(Ndim,Ndim) ) - - NST = NINT( DBLE(LTROT)/DBLE(NWRAP) ) - !Write(6,*) "Write UL ", NST - Do nf = 1,N_FL - CALL INITD(UL(:,:,Nf),Z_ONE) - do I = 1,Ndim - DL(I,Nf) = Z_ONE - enddo - CALL INITD(VL(:,:,nf),Z_ONE) - DO I = 1,NDim - DO J = 1,NDim - UST(I,J,NST,nf) = UL(I,J,nf) - VST(I,J,NST,nf) = VL(I,J,nf) - ENDDO - ENDDO - DO I = 1,NDim - DST(I,NST,nf) = DL(I,nf) - ENDDO - - CALL INITD(UR(:,:,nf),Z_ONE) - CALL INITD(VR(:,:,nf),Z_ONE) - Do I = 1,Ndim - DR(I,nf) = Z_ONE - Enddo - Enddo - - DO NT = LTROT-NWRAP,NWRAP,-1 - IF ( MOD(NT,NWRAP) == 0 ) THEN - NT1 = NT + NWRAP - !Write(6,*) 'Calling Wrapul:', NT1,NT - CALL WRAPUL(NT1,NT,UL,DL, VL) - NST = NINT( DBLE(NT)/DBLE(NWRAP) ) - !Write(6,*) "Write UL ", NST - Do nf = 1,N_FL - DO I = 1,Ndim - DO J = 1,Ndim - UST(I,J,NST,nf) = UL(I,J,nf) - VST(I,J,NST,nf) = VL(I,J,nf) - ENDDO - ENDDO - DO I = 1,Ndim - DST(I,NST,nf) = DL(I,nf) - ENDDO - ENDDO - ENDIF - ENDDO - CALL WRAPUL(NWRAP,0, UL ,DL, VL) - - !WRITE(6,*) 'Filling up storage' - !Write(6,*) 'Done wrapping' - NVAR = 1 - Phase = cmplx(1.d0,0.d0) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf) ) - Phase = Phase*Z - Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) -#ifdef MPI - WRITE(6,*) 'Phase is: ', Irank, PHASE, GR(1,1,1) -#else - WRITE(6,*) 'Phase is: ', PHASE -!!$ if (N_FL == 1) then -!!$ Do n = 1,Ndim -!!$ Write(6,*) GR(1,n,1) -!!$ enddo -!!$ else -!!$ Do n = 1,Ndim -!!$ Write(6,*) GR(1,n,1), GR(1,n,2) -!!$ enddo -!!$ endif -#endif - - Call Control_init - - DO NBC = 1, NBIN - ! Here, you have the green functions on time slice 1. - ! Set bin observables to zero. - - Call Init_obs(Ltau) - DO NSW = 1, NSWEEP - - !Propagation from 1 to Ltrot - !Set the right storage to 1 - - do nf = 1,N_FL - CALL INITD(UR(:,:,nf),Z_ONE) - CALL INITD(VR(:,:,nf),Z_ONE) - do n = 1,Ndim - DR(n,nf)= Z_ONE - Enddo - Enddo - - DO NTAU = 0, LTROT-1 - NTAU1 = NTAU + 1 - !Write(6,*) "Hi" - CALL WRAPGRUP(GR,NTAU,PHASE) - !Write(6,*) "Hi1" - IF ( MOD(NTAU1,NWRAP ) .EQ. 0 ) THEN - NST = NINT( DBLE(NTAU1)/DBLE(NWRAP) ) - NT1 = NTAU1 - NWRAP - CALL WRAPUR(NT1, NTAU1,UR, DR, VR) - Z = cmplx(1.d0,0.d0) - Do nf = 1, N_FL - DO J = 1,Ndim - DO I = 1,Ndim - UL(I,J,nf) = UST(I,J,NST,nf) - VL(I,J,nf) = VST(I,J,NST,nf) - ENDDO - ENDDO - DO I = 1,Ndim - DL(I,nf) = DST(I,NST,nf) - ENDDO - ! Write in store Right prop from 1 to LTROT/NWRAP - !Write(6,*) 'Write UR, read UL ', NTAU1, NST - DO J = 1,Ndim - DO I = 1,Ndim - UST(I,J,NST,nf) = UR(I,J,nf) - VST(I,J,NST,nf) = VR(I,J,nf) - ENDDO - ENDDO - DO I = 1,Ndim - DST(I,NST,nf) = DR(I,nf) - ENDDO - NVAR = 1 - IF (NTAU1 .GT. LTROT/2) NVAR = 2 - !Write(6,*) ' Call Cgr' - do J = 1,Ndim - do I = 1,Ndim - TEST(I,J) = GR(I,J,nf) - enddo - enddo - CALL CGR(Z1, NVAR, GR(:,:,nf), UR(:,:,nf),DR(:,nf),VR(:,:,nf),UL(:,:,nf),DL(:,nf),VL(:,:,nf) ) - Z = Z*Z1 - !Write(6,*) 'Calling control ',NTAU1, Z1 - Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) - ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) - Call Control_PrecisionP(Z,Phase) - Phase = Z - ENDIF - - IF (NTAU1.GE. LOBS_ST .AND. NTAU1.LE. LOBS_EN ) THEN - !Write(6,*) 'Call obser ', Ntau1 - CALL Obser( GR, PHASE, Ntau1 ) - !Write(6,*) 'Return obser' - ENDIF - !Write(6,*) NTAU1 - ENDDO - - Do nf = 1,N_FL - CALL INITD(UL(:,:,nf),Z_ONE) - CALL INITD(VL(:,:,nf),Z_ONE) - Do n = 1,Ndim - DL(n,nf) = Z_ONE - Enddo - ENDDO - - DO NTAU = LTROT,1,-1 - NTAU1 = NTAU - 1 - CALL WRAPGRDO(GR,NTAU, PHASE) - IF (NTAU1.GE. LOBS_ST .AND. NTAU1.LE. LOBS_EN ) THEN - CALL Obser( GR, PHASE, Ntau1 ) - ENDIF - IF ( MOD(NTAU1,NWRAP).EQ.0 .AND. NTAU1.NE.0 ) THEN - ! WRITE(50,*) 'Recalc at :', NTAU1 - NST = NINT( DBLE(NTAU1)/DBLE(NWRAP) ) - NT1 = NTAU1 + NWRAP - !Write(6,*) 'Wrapul : ', NT1, NTAU1 - CALL WRAPUL(NT1,NTAU1, UL, DL, VL ) - !Write(6,*) 'Write UL, read UR ', NTAU1, NST - Z = cmplx(1.d0,0.d0) - do nf = 1,N_FL - DO J = 1,Ndim - DO I = 1,Ndim - UR(I,J,nf) = UST(I,J,NST,nf) - VR(I,J,nf) = VST(I,J,NST,nf) - ENDDO - ENDDO - DO I = 1,Ndim - DR(I,nf) = DST(I,NST,nf) - ENDDO - ! WRITE in store the left prop. from LTROT/NWRAP-1 to 1 - DO J = 1,Ndim - DO I = 1,Ndim - UST(I,J,NST,nf) = UL(I,J,nf) - VST(I,J,NST,nf) = VL(I,J,nf) - ENDDO - ENDDO - DO I = 1,Ndim - DST(I,NST,nf) = DL(I,nf) - ENDDO - NVAR = 1 - IF (NTAU1 .GT. LTROT/2) NVAR = 2 - !Write(6,*) ' Call Cgr' - do J = 1,Ndim - do I = 1,Ndim - TEST(I,J) = GR(I,J,nf) - enddo - enddo - CALL CGR(Z1, NVAR, GR(:,:,nf), UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf) ) - Z = Z*Z1 - !Write(6,*) 'Calling control: ', NTAU1, Z1 - Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) - ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) - Call Control_PrecisionP(Z,Phase) - Phase = Z - ENDIF - ENDDO - - !Calculate and compare green functions on time slice 0. - NT1 = 0 - CALL WRAPUL(NWRAP,NT1, UL, DL, VL ) - - do nf = 1,N_FL - CALL INITD(UR(:,:,nf),Z_ONE) - CALL INITD(VR(:,:,nf),Z_ONE) - DO I = 1,Ndim - DR(I,nf) = Z_ONE - ENDDO - ENDDO - Z = cmplx(1.d0,0.d0) - do nf = 1,N_FL - do J = 1,Ndim - do I = 1,Ndim - TEST(I,J) = GR(I,J,nf) - enddo - enddo - NVAR = 1 - CALL CGR(Z1, NVAR, GR(:,:,nf), UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf) ) - Z = Z*Z1 - !Write(6,*) 'Calling control 0', Z1 - Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) - ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) - Call Control_PrecisionP(Z,Phase) - Phase = Z - NST = NINT( DBLE(LTROT)/DBLE(NWRAP) ) - Do nf = 1,N_FL - DO I = 1,Ndim - DO J = 1,Ndim - UST(I,J,NST,nf) = CMPLX(0.D0,0.D0) - VST(I,J,NST,nf) = CMPLX(0.D0,0.D0) - ENDDO - ENDDO - DO I = 1,Ndim - DST(I ,NST,nf) = CMPLX(1.D0,0.D0) - UST(I,I,NST,nf) = CMPLX(1.D0,0.D0) - VST(I,I,NST,nf) = CMPLX(1.D0,0.D0) - ENDDO - enddo - IF ( LTAU == 1 ) then -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'Calling Tau_m', NWRAP, NSTM -!!$#else -!!$ Write(6,*) 'Calling Tau_m', NWRAP, NSTM -!!$#endif - - Call TAU_M( UST,DST,VST, GR, PHASE, NSTM, NWRAP ) -!!$#ifdef MPI -!!$ Write(6,*) Irank, 'Back Calling Tau_m' -!!$#else -!!$ Write(6,*) 'Back Calling Tau_m' -!!$#endif - endif - - ENDDO - Call Pr_obs(Ltau) - Call confout - Enddo - Call Control_Print - -#ifdef MPI - CALL MPI_FINALIZE(ierr) -#endif - -end Program Main diff --git a/Prog_8/nranf.f90 b/Prog_8/nranf.f90 deleted file mode 100644 index 4662b0f63..000000000 --- a/Prog_8/nranf.f90 +++ /dev/null @@ -1,12 +0,0 @@ - integer function nranf(N) - Use Random_wrap - implicit none - integer :: N - - nranf = nint(ranf()*dble(N) + 0.5) - - if (nranf .lt. 1 ) nranf = 1 - if (nranf .gt. N ) nranf = N - - end function nranf - diff --git a/Prog_8/outconfc.f90 b/Prog_8/outconfc.f90 deleted file mode 100644 index 1a8ec4e53..000000000 --- a/Prog_8/outconfc.f90 +++ /dev/null @@ -1,57 +0,0 @@ - SUBROUTINE confout - - Use Hamiltonian - - Implicit none - -#include "machine" - - -#ifdef MPI - INCLUDE 'mpif.h' - ! Local -#endif - - Integer :: I, IERR, ISIZE, IRANK, seed_in, K, iseed, Nt, nr - Integer, dimension(:), allocatable :: Seed_vec - Real (Kind=8) :: X - Logical :: lconf - character (len=64) :: file_sr, File_tg - -#ifdef MPI - INTEGER :: STATUS(MPI_STATUS_SIZE) - CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) - CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) - - Call Get_seed_Len(K) - Allocate(Seed_vec(K)) - Call Ranget(Seed_vec) - file_sr = "confout" - file_tg = File_i(file_sr,IRANK) - Open (Unit = 10, File=File_tg, status='unknown', ACTION='write') - Write(10,*) Seed_vec - do NT = 1,LTROT - do I = 1,Size(Nsigma,1) - write(10,*) NSIGMA(I,NT) - enddo - enddo - close(10) - Deallocate(Seed_vec) - -#else - Call Get_seed_Len(K) - Allocate(Seed_vec(K)) - Call Ranget(Seed_vec) - file_tg = "confout_0" - Open (Unit = 10, File=File_tg, status='unknown', ACTION='write') - Write(10,*) Seed_vec - do NT = 1,LTROT - do I = 1,Size(Nsigma,1) - write(10,*) Nsigma(I,NT) - enddo - enddo - close(10) - Deallocate(Seed_vec) -#endif - - END SUBROUTINE CONFOUT diff --git a/Prog_8/tau_m.f90 b/Prog_8/tau_m.f90 deleted file mode 100644 index b1056e18d..000000000 --- a/Prog_8/tau_m.f90 +++ /dev/null @@ -1,236 +0,0 @@ - Module Tau_m_mod - - Use Hamiltonian - Use Operator_mod - Use Precdef - Use Control - Use Hop_mod - - Contains - - SUBROUTINE TAU_M( UST,DST,VST, GR, PHASE, NSTM, NWRAP ) - - Implicit none - - Interface - SUBROUTINE WRAPUL(NTAU1, NTAU, UL ,DL, VL) - Use Hamiltonian - Implicit none - COMPLEX (KIND=8) :: UL(Ndim,Ndim,N_FL), VL(Ndim,Ndim,N_FL) - COMPLEX (KIND=8) :: DL(Ndim,N_FL) - Integer :: NTAU1, NTAU - END SUBROUTINE WRAPUL - SUBROUTINE CGR2_2(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ) - Use Precdef - Use MyMats - Use UDV_WRAP_mod - Implicit none - - ! Arguments - Integer, intent(in) :: LQ - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - end SUBROUTINE CGR2_2 - SUBROUTINE CGR2_1(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ, NVAR) - Use Precdef - Use MyMats - USe UDV_Wrap_mod - Implicit none - ! Arguments - Integer, intent(in) :: LQ, NVAR - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - end SUBROUTINE CGR2_1 - SUBROUTINE CGR2(GRT0, GR00, GRTT, GR0T, U2, D2, V2, U1, D1, V1, LQ) - - ! B2 = U2*D2*V2 - ! B1 = V1*D1*U1 - !Calc: ( 1 B1 )^-1 i.e. 2*LQ \times 2*LQ matrix - ! (-B2 1 ) - - - Use Precdef - Use UDV_WRAP_mod - Use MyMats - - Implicit none - - ! Arguments - Integer :: LQ - Complex (Kind=double), intent(in) :: U1(LQ,LQ), V1(LQ,LQ), U2(LQ,LQ), V2(LQ,LQ) - Complex (Kind=double), intent(in) :: D2(LQ), D1(LQ) - Complex (Kind=double), intent(inout) :: GRT0(LQ,LQ), GR0T(LQ,LQ), GR00(LQ,LQ), GRTT(LQ,LQ) - end SUBROUTINE CGR2 - end Interface - - Complex (Kind=double), Intent(in) :: UST(NDIM,NDIM,NSTM,N_FL), VST(NDIM,NDIM,NSTM,N_FL), DST(NDIM,NSTM,N_FL) - Complex (Kind=double), Intent(in) :: GR(NDIM,NDIM,N_FL), Phase - Integer, Intent(In) :: NSTM, NWRAP - - - ! Local - ! This could be placed as private for the module - Complex (Kind=double) :: GT0(NDIM,NDIM,N_FL), G00(NDIM,NDIM,N_FL), GTT(NDIM,NDIM,N_FL), G0T(NDIM,NDIM,N_FL) - Complex (Kind=double) :: UL(Ndim,Ndim,N_FL), DL(Ndim,N_FL), VL(Ndim,Ndim,N_FL) - Complex (Kind=double) :: UR(Ndim,Ndim,N_FL), DR(Ndim,N_FL), VR(Ndim,Ndim,N_FL) - Complex (Kind=double) :: HLP4(Ndim,Ndim), HLP5(Ndim,Ndim), HLP6(Ndim,Ndim) - - Complex (Kind=double) :: Z - Integer :: I, J, nf, NT, NT1, NTST, NST, NVAR - - !Tau = 0 - Do nf = 1, N_FL - DO J = 1,Ndim - DO I = 1,Ndim - Z = cmplx(0.d0,0.d0) - if (I == J ) Z = cone - G00(I,J,nf) = GR(I,J,nf) - GT0(I,J,nf) = GR(I,J,nf) - GTT(I,J,nf) = GR(I,J,nf) - G0T(I,J,nf) = -(Z - GR(I,J,nf)) - ENDDO - ENDDO - Enddo - NT = 0 - ! In Module Hamiltonian - CALL OBSERT(NT, GT0,G0T,G00,GTT, PHASE) - - Do nf = 1, N_FL - CALL INITD(UR(:,:,nf),cone) - CALL INITD(VR(:,:,nf),cone) - enddo - DR = cone - - - DO NT = 0,LTROT - 1 - ! Now wrapup: - NT1 = NT + 1 - CALL PROPR (GT0,NT1) - CALL PROPRM1 (G0T,NT1) - CALL PROPRM1 (GTT,NT1) - CALL PROPR (GTT,NT1) - ! In Module Hamiltonian - CALL OBSERT(NT1, GT0,G0T,G00,GTT,PHASE) - - IF ( MOD(NT1,NWRAP).EQ.0 .AND. NT1.NE.LTROT ) THEN - NTST = NT1 - NWRAP - NST = NT1/(NWRAP) - ! WRITE(6,*) 'NT1, NST: ', NT1,NST - CALL WRAPUR(NTST, NT1,UR, DR, VR) - DO nf = 1,N_FL - DO J = 1,NDIM - DO I = 1,NDIM - UL(I,J,nf) = UST(I,J,NST,nf) - VL(I,J,nf) = VST(I,J,NST,nf) - ENDDO - ENDDO - DO I = 1,NDIM - DL(I,nf) = DST(I,NST,nf) - ENDDO - Enddo - Do nf = 1,N_FL - Do J = 1,Ndim - DO I = 1,Ndim - HLP4(I,J) = GTT(I,J,nf) - HLP5(I,J) = GT0(I,J,nf) - HLP6(I,J) = G0T(I,J,nf) - Enddo - Enddo - NVAR = 1 - IF (NT1 > LTROT/2) NVAR = 2 - !DO I = 1,Ndim - ! Write(6,*) DL(I,nf)*DR(I,nf) - !enddo - !Write(6,*) 'Call CGR2' - Call CGR2_2(GT0(:,:,nf), G00(:,:,nf), GTT(:,:,nf), G0T(:,:,nf), & - & UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf),NDIM) - !Call CGR2 (GT0(:,:,nf), G00(:,:,nf), GTT(:,:,nf), G0T(:,:,nf), & - ! & UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf),NDIM) - - !Call CGR2_1(GT0(:,:,nf), G00(:,:,nf), GTT(:,:,nf), G0T(:,:,nf), & - ! & UR(:,:,nf),DR(:,nf),VR(:,:,nf), UL(:,:,nf),DL(:,nf),VL(:,:,nf),NDIM,NVAR) - - !Write(6,*) 'End Call CGR2' - !Write(6,*) ' Tau ', NT1 - !Write(6,*) ' G00 ' - Call Control_Precision_tau(GR(:,:,nf), G00(:,:,nf), Ndim) - !Write(6,*) ' GTT ' - Call Control_Precision_tau(HLP4 , GTT(:,:,nf), Ndim) - !Write(6,*) ' GT0 ' - Call Control_Precision_tau(HLP5 , GT0(:,:,nf), Ndim) - !Write(6,*) ' G0T ' - Call Control_Precision_tau(HLP6 , G0T(:,:,nf), Ndim) - Enddo - Endif - ENDDO - - END SUBROUTINE TAU_M - -!============================================================== - - SUBROUTINE PROPR(AIN,NT) - - ! Ain = B(NT-1, NT1) - ! Aout= Ain = B(NT , NT1) - - Implicit none - Complex (Kind=double), intent(INOUT) :: Ain(Ndim,Ndim,N_FL) - Integer, INTENT(IN) :: NT - - !Locals - Integer :: J,I,nf,n - Complex (Kind=double) :: HLP4(Ndim,Ndim) - Real (Kind=double) :: X - - Do nf = 1,N_FL - !CALL MMULT(HLP4,Exp_T(:,:,nf) ,Ain(:,:,nf)) - Call Hop_mod_mmthr(Ain(:,:,nf),HLP4,nf) - Do n = 1,Size(Op_V,1) - X = Phi(nsigma(n,nt),Op_V(n,nf)%type) - Call Op_mmultR(HLP4,Op_V(n,nf),X,Ndim) - ENDDO - Do J = 1,Ndim - do I = 1,Ndim - Ain(I,J,nf) = HLP4(I,J) - enddo - ENDDO - Enddo - - end SUBROUTINE PROPR -!============================================================== - SUBROUTINE PROPRM1(AIN,NT) - - !Ain = B^{-1}(NT-1, NT1) - !Aout= B^{-1}(NT , NT1) - - - Implicit none - - !Arguments - Complex (Kind=double), intent(Inout) :: AIN(Ndim, Ndim, N_FL) - Integer :: NT - - ! Locals - Integer :: J,I,nf,n - Complex (Kind=double) :: HLP4(Ndim,Ndim) - Real (Kind=double) :: X - - do nf = 1,N_FL - !Call MMULT(HLP4,Ain(:,:,nf),Exp_T_M1(:,:,nf) ) - Call Hop_mod_mmthl_m1(Ain(:,:,nf),HLP4,nf) - Do n =1,Size(Op_V,1) - X = -Phi(nsigma(n,nt),Op_V(n,nf)%type) - Call Op_mmultL(HLP4,Op_V(n,nf),X,Ndim) - Enddo - Do J = 1,Ndim - do I = 1,Ndim - Ain(I,J,nf) = HLP4(I,J) - enddo - Enddo - enddo - - END SUBROUTINE PROPRM1 -!============================================================== - end Module Tau_m_mod diff --git a/Prog_8/wrapul.f90 b/Prog_8/wrapul.f90 deleted file mode 100644 index 8d93cb17c..000000000 --- a/Prog_8/wrapul.f90 +++ /dev/null @@ -1,77 +0,0 @@ - SUBROUTINE WRAPUL(NTAU1, NTAU, UL ,DL, VL) - - !Given B(LTROT,NTAU1,Nf ) = VLUP, DLUP, ULUP - !Returns B(LTROT,NTAU, Nf ) = VLUP, DLUP, ULUP - - - !NOTE: NTAU1 > NTAU. - ! Does this for all replicas - Use Hamiltonian - Use Hop_mod - Use UDV_Wrap_mod - - Implicit none - - ! Arguments - COMPLEX (KIND=8) :: UL(Ndim,Ndim,N_FL), VL(Ndim,Ndim,N_FL) - COMPLEX (KIND=8) :: DL(Ndim,N_FL) - Integer :: NTAU1, NTAU - - - ! Working space. - COMPLEX (Kind=8) :: U(Ndim,Ndim), U1(Ndim,Ndim), V1(Ndim,Ndim), TMP(Ndim,Ndim), TMP1(Ndim,Ndim) - COMPLEX (Kind=8) :: D1(Ndim), Z_ONE - Integer :: I, J, NT, NCON, nr, n, nf - Real (Kind=8) :: X - - - - NCON = 0 ! Test for UDV :::: 0: Off, 1: On. - - Z_ONE = cmplx(1.d0,0.d0) - Do nf = 1, N_FL - CALL INITD(TMP,Z_ONE) - DO NT = NTAU1, NTAU+1 , -1 - Do n = Size(Op_V,1),1,-1 - X = Phi(nsigma(n,nt),Op_V(n,nf)%type) - Call Op_mmultL(Tmp,Op_V(n,nf),X,Ndim) - enddo - !CALL MMULT( TMP1,Tmp,Exp_T(:,:,nf) ) - Call Hop_mod_mmthl (Tmp, Tmp1,nf) - Tmp = Tmp1 - ENDDO - - !Carry out U,D,V decomposition. - DO J = 1,NDim - DO I = 1,NDim - TMP1(I,J) = CONJG( TMP(J,I) ) - U (I,J) = CONJG( UL (J,I,nf) ) - ENDDO - ENDDO - CALL MMULT(TMP,TMP1,U) - DO J = 1,NDim - DO I = 1,NDim - TMP(I,J) = TMP(I,J)*DL(J,nf) - ENDDO - ENDDO - CALL UDV_WRAP(TMP,U1,D1,V1,NCON) - !CALL UDV(TMP,U1,D1,V1,NCON) - DO J = 1,NDim - DO I = 1,NDim - UL (I,J,nf) = CONJG( U1(J,I) ) - TMP(I,J) = CONJG( V1(J,I) ) - ENDDO - ENDDO - CALL MMULT(TMP1,VL(:,:,nf),TMP) - DO J = 1,NDim - DO I = 1,NDim - VL(I,J,nf) = TMP1(I,J) - ENDDO - ENDDO - DO I = 1,NDim - DL(I,nf) = D1(I) - ENDDO - ENDDO - - END SUBROUTINE WRAPUL - diff --git a/Prog_8/wrapur.f90 b/Prog_8/wrapur.f90 deleted file mode 100644 index c53b28ce3..000000000 --- a/Prog_8/wrapur.f90 +++ /dev/null @@ -1,48 +0,0 @@ - SUBROUTINE WRAPUR(NTAU, NTAU1, UR, DR, VR) - - ! Given B(NTAU, 1 ) = UR, DR, VR - ! Returns B(NTAU1, 1 ) = UR, DR, VR - ! NOTE: NTAU1 > NTAU. - - Use Hamiltonian - Use UDV_Wrap_mod - Use Hop_mod - Implicit None - - ! Arguments - COMPLEX (KIND=8) :: UR(Ndim,Ndim,N_FL), VR(Ndim,Ndim,N_FL) - COMPLEX (KIND=8) :: DR(Ndim,N_FL) - Integer :: NTAU1, NTAU - - - ! Working space. - Complex (Kind=8) :: Z_ONE - COMPLEX (Kind=8) :: V1(Ndim,Ndim), TMP(Ndim,Ndim), TMP1(Ndim,Ndim) - Integer ::NCON, NT, I, J, n, nf - Real (Kind=8) :: X - - NCON = 0 ! Test for UDV :::: 0: Off, 1: On. - Z_ONE = cmplx(1.d0,0.d0) - - Do nf = 1,N_FL - CALL INITD(TMP,Z_ONE) - DO NT = NTAU + 1, NTAU1 - !CALL MMULT(TMP1,Exp_T(:,:,nf) ,TMP) - Call Hop_mod_mmthr(TMP,TMP1,nf) - TMP = TMP1 - Do n = 1,Size(Op_V,1) - X = Phi(nsigma(n,nt),Op_V(n,nf)%type) - Call Op_mmultR(Tmp,Op_V(n,nf),X,Ndim) - ENDDO - ENDDO - CALL MMULT(TMP1,TMP,UR(:,:,nf)) - DO J = 1,NDim - DO I = 1,NDim - TMP1(I,J) = TMP1(I,J)*DR(J,nf) - TMP(I,J) = VR(I,J,nf) - ENDDO - ENDDO - CALL UDV_WRAP(TMP1,UR(:,:,nf),DR(:,nf),V1,NCON) - CALL MMULT(VR(:,:,nf),V1,TMP) - ENDDO - END SUBROUTINE WRAPUR diff --git a/cmake/Modules/FindOpenMP_Fortran.cmake b/cmake/Modules/FindOpenMP_Fortran.cmake new file mode 100644 index 000000000..bca09d058 --- /dev/null +++ b/cmake/Modules/FindOpenMP_Fortran.cmake @@ -0,0 +1,106 @@ +# - Finds OpenMP support +# This module can be used to detect OpenMP support in a compiler. +# If the compiler supports OpenMP, the flags required to compile with +# openmp support are set. +# +# This module was modified from the standard FindOpenMP module to find Fortran +# flags. +# +# The following variables are set: +# OpenMP_Fortran_FLAGS - flags to add to the Fortran compiler for OpenMP +# support. In general, you must use these at both +# compile- and link-time. +# OMP_NUM_PROCS - the max number of processors available to OpenMP + +#============================================================================= +# Copyright 2009 Kitware, Inc. +# Copyright 2008-2009 André Rigland Brodtkorb +# +# Distributed under the OSI-approved BSD License (the "License"); +# see accompanying file Copyright.txt for details. +# +# This software is distributed WITHOUT ANY WARRANTY; without even the +# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +# See the License for more information. +#============================================================================= +# (To distribute this file outside of CMake, substitute the full +# License text for the above reference.) + +INCLUDE (${CMAKE_ROOT}/Modules/FindPackageHandleStandardArgs.cmake) + +SET (OpenMP_Fortran_FLAG_CANDIDATES + #Microsoft Visual Studio + "/openmp" + #Intel windows + "/Qopenmp" + #Intel + "-openmp" + # Intel 16+ + "-qopenmp" + #Gnu + "-fopenmp" + #Empty, if compiler automatically accepts openmp + " " + #Sun + "-xopenmp" + #HP + "+Oopenmp" + #IBM XL C/c++ + "-qsmp" + #Portland Group + "-mp" +) + +IF (DEFINED OpenMP_Fortran_FLAGS) + SET (OpenMP_Fortran_FLAG_CANDIDATES) +ENDIF (DEFINED OpenMP_Fortran_FLAGS) + +# check fortran compiler. also determine number of processors +FOREACH (FLAG ${OpenMP_Fortran_FLAG_CANDIDATES}) + SET (SAFE_CMAKE_REQUIRED_FLAGS "${CMAKE_REQUIRED_FLAGS}") + SET (CMAKE_REQUIRED_FLAGS "${FLAG}") + UNSET (OpenMP_FLAG_DETECTED CACHE) + MESSAGE (STATUS "Try OpenMP Fortran flag = [${FLAG}]") + FILE (WRITE "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/testFortranOpenMP.f90" +" +program TestOpenMP + use omp_lib + write(*,'(I2)',ADVANCE='NO') omp_get_num_procs() +end program TestOpenMP +") + SET (MACRO_CHECK_FUNCTION_DEFINITIONS + "-DOpenMP_FLAG_DETECTED ${CMAKE_REQUIRED_FLAGS}") + TRY_RUN (OpenMP_RUN_FAILED OpenMP_FLAG_DETECTED ${CMAKE_BINARY_DIR} + ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/testFortranOpenMP.f90 + COMPILE_DEFINITIONS ${CMAKE_REQUIRED_DEFINITIONS} + CMAKE_FLAGS -DCOMPILE_DEFINITIONS:STRING=${MACRO_CHECK_FUNCTION_DEFINITIONS} + COMPILE_OUTPUT_VARIABLE OUTPUT + RUN_OUTPUT_VARIABLE OMP_NUM_PROCS_INTERNAL) + IF (OpenMP_FLAG_DETECTED) + FILE (APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeOutput.log + "Determining if the Fortran compiler supports OpenMP passed with " + "the following output:\n${OUTPUT}\n\n") + SET (OpenMP_FLAG_DETECTED 1) + IF (OpenMP_RUN_FAILED) + MESSAGE (FATAL_ERROR "OpenMP found, but test code did not run") + ENDIF (OpenMP_RUN_FAILED) + SET (OMP_NUM_PROCS ${OMP_NUM_PROCS_INTERNAL} CACHE + STRING "Number of processors OpenMP may use" FORCE) + SET (OpenMP_Fortran_FLAGS_INTERNAL "${FLAG}") + BREAK () + ELSE () + FILE (APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log + "Determining if the Fortran compiler supports OpenMP failed with " + "the following output:\n${OUTPUT}\n\n") + SET (OpenMP_FLAG_DETECTED 0) + ENDIF (OpenMP_FLAG_DETECTED) +ENDFOREACH (FLAG ${OpenMP_Fortran_FLAG_CANDIDATES}) + +SET (OpenMP_Fortran_FLAGS "${OpenMP_Fortran_FLAGS_INTERNAL}" + CACHE STRING "Fortran compiler flags for OpenMP parallization") + +# handle the standard arguments for FIND_PACKAGE +FIND_PACKAGE_HANDLE_STANDARD_ARGS (OpenMP_Fortran DEFAULT_MSG + OpenMP_Fortran_FLAGS) + +MARK_AS_ADVANCED(OpenMP_Fortran_FLAGS) diff --git a/cmake/Modules/SetCompileFlag.cmake b/cmake/Modules/SetCompileFlag.cmake new file mode 100644 index 000000000..04ff3ffbd --- /dev/null +++ b/cmake/Modules/SetCompileFlag.cmake @@ -0,0 +1,112 @@ +############################################################################# +# Given a list of flags, this function will try each, one at a time, +# and choose the first flag that works. If no flags work, then nothing +# will be set, unless the REQUIRED key is given, in which case an error +# will be given. +# +# Call is: +# SET_COMPILE_FLAG(FLAGVAR FLAGVAL (Fortran|C|CXX) flag1 flag2...) +# +# For example, if you have the flag CMAKE_C_FLAGS and you want to add +# warnings and want to fail if this is not possible, you might call this +# function in this manner: +# SET_COMPILE_FLAGS(CMAKE_C_FLAGS "${CMAKE_C_FLAGS}" C REQUIRED +# "-Wall" # GNU +# "-warn all" # Intel +# ) +# The optin "-Wall" will be checked first, and if it works, will be +# appended to the CMAKE_C_FLAGS variable. If it doesn't work, then +# "-warn all" will be tried. If this doesn't work then checking will +# terminate because REQUIRED was given. +# +# The reasong that the variable must be given twice (first as the name then +# as the value in quotes) is because of the way CMAKE handles the passing +# of variables in functions; it is difficult to extract a variable's +# contents and assign new values to it from within a function. +############################################################################# + +INCLUDE(${CMAKE_ROOT}/Modules/CheckCCompilerFlag.cmake) +INCLUDE(${CMAKE_ROOT}/Modules/CheckCXXCompilerFlag.cmake) + +FUNCTION(SET_COMPILE_FLAG FLAGVAR FLAGVAL LANG) + + # Do some up front setup if Fortran + IF(LANG STREQUAL "Fortran") + # Create a list of error messages from compilers + SET(FAIL_REGEX + "ignoring unknown option" # Intel + "invalid argument" # Intel + "unrecognized .*option" # GNU + "[Uu]nknown switch" # Portland Group + "ignoring unknown option" # MSVC + "warning D9002" # MSVC, any lang + "[Uu]nknown option" # HP + "[Ww]arning: [Oo]ption" # SunPro + "command option .* is not recognized" # XL + ) + ENDIF(LANG STREQUAL "Fortran") + + # Make a variable holding the flags. Filter out REQUIRED if it is there + SET(FLAG_REQUIRED FALSE) + SET(FLAG_FOUND FALSE) + UNSET(FLAGLIST) + FOREACH (var ${ARGN}) + STRING(TOUPPER "${var}" UP) + IF(UP STREQUAL "REQUIRED") + SET(FLAG_REQUIRED TRUE) + ELSE() + SET(FLAGLIST ${FLAGLIST} "${var}") + ENDIF(UP STREQUAL "REQUIRED") + ENDFOREACH (var ${ARGN}) + + # Now, loop over each flag + FOREACH(flag ${FLAGLIST}) + + UNSET(FLAG_WORKS) + # Check the flag for the given language + IF(LANG STREQUAL "C") + CHECK_C_COMPILER_FLAG("${flag}" FLAG_WORKS) + ELSEIF(LANG STREQUAL "CXX") + CHECK_CXX_COMPILER_FLAG("${flag}" FLAG_WORKS) + ELSEIF(LANG STREQUAL "Fortran") + # There is no nice function to do this for FORTRAN, so we must manually + # create a test program and check if it compiles with a given flag. + SET(TESTFILE "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}") + SET(TESTFILE "${TESTFILE}/CMakeTmp/testFortranFlags.f90") + FILE(WRITE "${TESTFILE}" +" +program dummyprog + i = 5 +end program dummyprog +") + TRY_COMPILE(FLAG_WORKS ${CMAKE_BINARY_DIR} ${TESTFILE} + COMPILE_DEFINITIONS "${flag}" OUTPUT_VARIABLE OUTPUT) + + # Check that the output message doesn't match any errors + FOREACH(rx ${FAIL_REGEX}) + IF("${OUTPUT}" MATCHES "${rx}") + SET(FLAG_WORKS FALSE) + ENDIF("${OUTPUT}" MATCHES "${rx}") + ENDFOREACH(rx ${FAIL_REGEX}) + + ELSE() + MESSAGE(FATAL_ERROR "Unknown language in SET_COMPILE_FLAGS: ${LANG}") + ENDIF(LANG STREQUAL "C") + + # If this worked, use these flags, otherwise use other flags + IF(FLAG_WORKS) + # Append this flag to the end of the list that already exists + SET(${FLAGVAR} "${FLAGVAL} ${flag}" CACHE STRING + "Set the ${FLAGVAR} flags" FORCE) + SET(FLAG_FOUND TRUE) + BREAK() # We found something that works, so exit + ENDIF(FLAG_WORKS) + + ENDFOREACH(flag ${FLAGLIST}) + + # Raise an error if no flag was found + IF(FLAG_REQUIRED AND NOT FLAG_FOUND) + MESSAGE(FATAL_ERROR "No compile flags were found") + ENDIF(FLAG_REQUIRED AND NOT FLAG_FOUND) + +ENDFUNCTION() diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake new file mode 100644 index 000000000..02fadb289 --- /dev/null +++ b/cmake/Modules/SetFortranFlags.cmake @@ -0,0 +1,185 @@ +###################################################### +# Determine and set the Fortran compiler flags we want +###################################################### + +#################################################################### +# Make sure that the default build type is RELEASE if not specified. +#################################################################### +INCLUDE(${CMAKE_MODULE_PATH}/SetCompileFlag.cmake) + +# Make sure the build type is uppercase +STRING(TOUPPER "${CMAKE_BUILD_TYPE}" BT) + +IF(BT STREQUAL "RELEASE") + SET(CMAKE_BUILD_TYPE RELEASE CACHE STRING + "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + FORCE) +ELSEIF(BT STREQUAL "DEBUG") + SET (CMAKE_BUILD_TYPE DEBUG CACHE STRING + "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + FORCE) +ELSEIF(BT STREQUAL "TESTING") + SET (CMAKE_BUILD_TYPE TESTING CACHE STRING + "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + FORCE) +ELSEIF(NOT BT) + SET(CMAKE_BUILD_TYPE RELEASE CACHE STRING + "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + FORCE) + MESSAGE(STATUS "CMAKE_BUILD_TYPE not given, defaulting to RELEASE") +ELSE() + MESSAGE(FATAL_ERROR "CMAKE_BUILD_TYPE not valid, choices are DEBUG, RELEASE, or TESTING") +ENDIF(BT STREQUAL "RELEASE") + +######################################################### +# If the compiler flags have already been set, return now +######################################################### + +IF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG) + RETURN () +ENDIF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG) + +######################################################################## +# Determine the appropriate flags for this compiler for each build type. +# For each option type, a list of possible flags is given that work +# for various compilers. The first flag that works is chosen. +# If none of the flags work, nothing is added (unless the REQUIRED +# flag is given in the call). This way unknown compiles are supported. +####################################################################### + +##################### +### GENERAL FLAGS ### +##################### + +# Don't add underscores in symbols for C-compatability +#SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" +# Fortran "-fno-underscoring") + +# There is some bug where -march=native doesn't work on Mac +IF(APPLE) + SET(GNUNATIVE "-mtune=native") +ELSE() + SET(GNUNATIVE "-march=native") +ENDIF() + +#THe folowing does not seem to get added... +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" + Fortran "-fimplicit-none" # GNU +) + +# Optimize for the host's architecture +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" + Fortran "-xHost" # Intel + "/QxHost" # Intel Windows + ${GNUNATIVE} # GNU + "-ta=host" # Portland Group + ) + +#SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" +# Fortran "-ffree-form" # GNU +# "-free" #Intel Linux and MacOS +# "/free" #Intel Windows +#) + +#SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" +# Fortran "-ffree-line-length-none" # GNU +#) + +#SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" +# Fortran "-cpp" # GNU +# "-fpp" # Intel +# "/fpp" # Intel Windows +# "-Mpreprocess" # Portland Group +# ) + +################### +### DEBUG FLAGS ### +################### + +# NOTE: debugging symbols (-g or /debug:full) are already on by default + +# Disable optimizations +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" + Fortran REQUIRED "-O0" # All compilers not on Windows + "/Od" # Intel Windows + ) + +# Turn on all warnings +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" + Fortran "-warn all" # Intel + "/warn:all" # Intel Windows + "-Wall" # GNU + # Portland Group (on by default) + ) + +# Traceback +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" + Fortran "-traceback" # Intel/Portland Group + "/traceback" # Intel Windows + "-fbacktrace" # GNU (gfortran) + "-ftrace=full" # GNU (g95) + ) + +# Check array bounds +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" + Fortran "-check bounds" # Intel + "/check:bounds" # Intel Windows + "-fcheck=bounds" # GNU (New style) + "-fbounds-check" # GNU (Old style) + "-Mbounds" # Portland Group + ) + +##################### +### TESTING FLAGS ### +##################### + +# Optimizations +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_TESTING "${CMAKE_Fortran_FLAGS_TESTING}" + Fortran REQUIRED "-O2" # All compilers not on Windows + "/O2" # Intel Windows + ) + +##################### +### RELEASE FLAGS ### +##################### + +# NOTE: agressive optimizations (-O3) are already turned on by default + + +# Unroll loops +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-funroll-loops" # GNU + "-unroll" # Intel + "/unroll" # Intel Windows + "-Munroll" # Portland Group + ) + +# Inline functions +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-inline" # Intel + "/Qinline" # Intel Windows + "-finline-functions" # GNU + "-Minline" # Portland Group + ) + +# Interprocedural (link-time) optimizations +#SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" +# Fortran "-ipo" # Intel +# "/Qipo" # Intel Windows +# "-flto" # GNU +# "-Mipa" # Portland Group +# ) + +# Single-file optimizations +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-ip" # Intel + "/Qip" # Intel Windows + ) + +# Vectorize code +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-vec-report0" # Intel + "-qopt-report0" # Intel 16+ + "/Qvec-report0" # Intel Windows + "-Mvect" # Portland Group + ) diff --git a/cmake/Modules/SetParallelizationLibrary.cmake b/cmake/Modules/SetParallelizationLibrary.cmake new file mode 100644 index 000000000..603d4299c --- /dev/null +++ b/cmake/Modules/SetParallelizationLibrary.cmake @@ -0,0 +1,39 @@ +# Turns on either OpenMP or MPI +# If both are requested, the other is disabled +# When one is turned on, the other is turned off +# If both are off, we explicitly disable them just in case + +IF (USE_OPENMP AND USE_MPI) + MESSAGE (FATAL_ERROR "Cannot use both OpenMP and MPI") +ELSEIF (USE_OPENMP) + # Find OpenMP + IF (NOT OpenMP_Fortran_FLAGS) + FIND_PACKAGE (OpenMP_Fortran) + IF (NOT OpenMP_Fortran_FLAGS) + MESSAGE (FATAL_ERROR "Fortran compiler does not support OpenMP") + ENDIF (NOT OpenMP_Fortran_FLAGS) + ENDIF (NOT OpenMP_Fortran_FLAGS) + # Turn of MPI + UNSET (MPI_FOUND CACHE) + UNSET (MPI_COMPILER CACHE) + UNSET (MPI_LIBRARY CACHE) +ELSEIF (USE_MPI) + # Find MPI + IF (NOT MPI_Fortran_FOUND) + FIND_PACKAGE (MPI REQUIRED) + ENDIF (NOT MPI_Fortran_FOUND) + # Turn off OpenMP + SET (OMP_NUM_PROCS 0 CACHE + STRING "Number of processors OpenMP may use" FORCE) + UNSET (OpenMP_C_FLAGS CACHE) + UNSET (GOMP_Fortran_LINK_FLAGS CACHE) +ELSE () + # Turn off both OpenMP and MPI + SET (OMP_NUM_PROCS 0 CACHE + STRING "Number of processors OpenMP may use" FORCE) + UNSET (OpenMP_Fortran_FLAGS CACHE) + UNSET (GOMP_Fortran_LINK_FLAGS CACHE) + UNSET (MPI_FOUND CACHE) + UNSET (MPI_COMPILER CACHE) + UNSET (MPI_LIBRARY CACHE) +ENDIF (USE_OPENMP AND USE_MPI) diff --git a/cmake/Modules/SetUpLAPACK.cmake b/cmake/Modules/SetUpLAPACK.cmake new file mode 100644 index 000000000..ae5bdea52 --- /dev/null +++ b/cmake/Modules/SetUpLAPACK.cmake @@ -0,0 +1,11 @@ +# Find LAPACK (finds BLAS also) if not already found +IF(NOT LAPACK_FOUND) + ENABLE_LANGUAGE(C) # Some libraries need a C compiler to find + FIND_PACKAGE(LAPACK REQUIRED) + # Remember that LAPACK (and BLAS) was found. For some reason the + # FindLAPACK routine doesn't place these into the CACHE. + SET(BLAS_FOUND TRUE CACHE INTERNAL "BLAS was found" FORCE) + SET(LAPACK_FOUND TRUE CACHE INTERNAL "LAPACK was found" FORCE) + SET(BLAS_LIBRARIES ${BLAS_LIBRARIES} CACHE INTERNAL "BLAS LIBS" FORCE) + SET(LAPACK_LIBRARIES ${LAPACK_LIBRARIES} CACHE INTERNAL "LAPACK LIBS" FORCE) +ENDIF(NOT LAPACK_FOUND) diff --git a/distclean.cmake b/distclean.cmake new file mode 100644 index 000000000..8e24f9e49 --- /dev/null +++ b/distclean.cmake @@ -0,0 +1,68 @@ +# This CMake script will delete build directories and files to bring the +# package back to it's distribution state + +# We want to start from the top of the source dir, so if we are in build +# we want to start one directory up +GET_FILENAME_COMPONENT(BASEDIR ${CMAKE_SOURCE_DIR} NAME) +IF(${BASEDIR} STREQUAL "build") + SET(TOPDIR "${CMAKE_SOURCE_DIR}/..") +ELSE() + SET(TOPDIR "${CMAKE_SOURCE_DIR}") +ENDIF() + +MACRO(GET_PARENT_DIRECTORIES search_string return_list grandparents) + FILE(GLOB_RECURSE new_list ${search_string}) + SET(dir_list "") + FOREACH(file_path ${new_list}) + GET_FILENAME_COMPONENT(dir_path ${file_path} PATH) + # Remove an extra directory component to return grandparent + IF(${grandparents}) + # Tack on a fake extension to trick CMake into removing a second + # path component + SET(dir_path "${dir_path}.tmp") + GET_FILENAME_COMPONENT(dir_path ${dir_path} PATH) + ENDIF(${grandparents}) + SET(dir_list ${dir_list} ${dir_path}) + ENDFOREACH() + LIST(REMOVE_DUPLICATES dir_list) + SET(${return_list} ${dir_list}) +ENDMACRO() + +# Find directories and files that we will want to remove +FILE(GLOB_RECURSE CMAKECACHE "${TOPDIR}/*CMakeCache.txt") +FILE(GLOB_RECURSE CMAKEINSTALL "${TOPDIR}/*cmake_install.cmake" + "${TOPDIR}/*install_manifest.txt") +FILE(GLOB_RECURSE MAKEFILE "${TOPDIR}/*Makefile") +FILE(GLOB_RECURSE CMAKETESTFILES "${TOPDIR}/*CTestTestfile.cmake") +SET(TOPDIRECTORIES "${TOPDIR}/lib" + "${TOPDIR}/test" + "${TOPDIR}/bin" +) + +# CMake has trouble finding directories recursively, so locate these +# files and then save the parent directory of the files +GET_PARENT_DIRECTORIES(Makefile.cmake CMAKEFILES 0) +GET_PARENT_DIRECTORIES(LastTest.log CMAKETESTING 1) + +# Place these files and directories into a list +SET(DEL ${TOPDIRECTORIES} + ${CMAKECACHE} + ${CMAKEINSTALL} + ${MAKEFILE} + ${CMAKEFILES} + ${CMAKETESTING} + ${CMAKETESTFILES} +) + +# If we are not in the build dir, delete that as well +IF(NOT (${BASEDIR} STREQUAL "build")) + FILE(GLOB BUILD "${TOPDIR}/build") + SET(DEL ${DEL} ${BUILD}) +ENDIF() + +# Loop over the directories and delete each one +FOREACH(D ${DEL}) + IF(EXISTS ${D}) + FILE(REMOVE_RECURSE ${D}) + ENDIF() +ENDFOREACH() diff --git a/src/Analysis/CMakeLists.txt b/src/Analysis/CMakeLists.txt new file mode 100644 index 000000000..76d395f25 --- /dev/null +++ b/src/Analysis/CMakeLists.txt @@ -0,0 +1,8 @@ +add_executable(jackv5 jackv5.f90) +add_executable(cov_eq cov_eq.F90) +add_executable(cov_tau cov_tau.f90) + +foreach(analysisbinary jackv5 cov_eq cov_tau) +target_link_libraries(${analysisbinary} ${MODULES} ${MYEIS} ${MYNAG} ${MYLIN}) +target_link_libraries(${analysisbinary} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}) +endforeach() diff --git a/Analysis_7/cov_eq.f90 b/src/Analysis/cov_eq.F90 similarity index 100% rename from Analysis_7/cov_eq.f90 rename to src/Analysis/cov_eq.F90 diff --git a/Analysis_7/cov_tau.f90 b/src/Analysis/cov_tau.f90 similarity index 100% rename from Analysis_7/cov_tau.f90 rename to src/Analysis/cov_tau.f90 diff --git a/Analysis_7/jackv5.f90 b/src/Analysis/jackv5.f90 similarity index 100% rename from Analysis_7/jackv5.f90 rename to src/Analysis/jackv5.f90 diff --git a/Libraries/Modules/BIDON b/src/Modules/BIDON similarity index 100% rename from Libraries/Modules/BIDON rename to src/Modules/BIDON diff --git a/src/Modules/CMakeLists.txt b/src/Modules/CMakeLists.txt new file mode 100644 index 000000000..190350489 --- /dev/null +++ b/src/Modules/CMakeLists.txt @@ -0,0 +1,4 @@ +#Modules library +SET(Modules_src ${SRCMODULES}/Files_mod.f90 ${SRCMODULES}/Histogram.f90 ${SRCMODULES}/Histogram_v2.f90 ${SRCMODULES}/Natural_constants.f90 ${SRCMODULES}/Random_Wrap.f90 ${SRCMODULES}/errors.f90 ${SRCMODULES}/fourier.f90 ${SRCMODULES}/lattices_v3.f90 +${SRCMODULES}/log_mesh.f90 ${SRCMODULES}/mat_mod.f90 ${SRCMODULES}/matrix.f90 ${SRCMODULES}/maxent.f90 ${SRCMODULES}/maxent_stoch.f90 ${SRCMODULES}/maxent_stoch.f90 ${SRCMODULES}/precdef.mod.f90 ${SRCMODULES}/smooth_stoch.f90 ${SRCMODULES}/tmp.f90) +ADD_LIBRARY(${MODULES} STATIC ${Modules_src}) diff --git a/Libraries/Modules/Files_mod.f90 b/src/Modules/Files_mod.f90 similarity index 100% rename from Libraries/Modules/Files_mod.f90 rename to src/Modules/Files_mod.f90 diff --git a/Libraries/Modules/Histogram.f90 b/src/Modules/Histogram.f90 similarity index 100% rename from Libraries/Modules/Histogram.f90 rename to src/Modules/Histogram.f90 diff --git a/Libraries/Modules/Histogram_v2.f90 b/src/Modules/Histogram_v2.f90 similarity index 100% rename from Libraries/Modules/Histogram_v2.f90 rename to src/Modules/Histogram_v2.f90 diff --git a/Libraries/Modules/Natural_constants.f90 b/src/Modules/Natural_constants.f90 similarity index 100% rename from Libraries/Modules/Natural_constants.f90 rename to src/Modules/Natural_constants.f90 diff --git a/Libraries/Modules/Random_Wrap.f90 b/src/Modules/Random_Wrap.f90 similarity index 100% rename from Libraries/Modules/Random_Wrap.f90 rename to src/Modules/Random_Wrap.f90 diff --git a/Libraries/Modules/errors.f90 b/src/Modules/errors.f90 similarity index 100% rename from Libraries/Modules/errors.f90 rename to src/Modules/errors.f90 diff --git a/Libraries/Modules/fourier.f90 b/src/Modules/fourier.f90 similarity index 100% rename from Libraries/Modules/fourier.f90 rename to src/Modules/fourier.f90 diff --git a/Libraries/Modules/lattices_v3.f90 b/src/Modules/lattices_v3.f90 similarity index 100% rename from Libraries/Modules/lattices_v3.f90 rename to src/Modules/lattices_v3.f90 diff --git a/Libraries/Modules/log_mesh.f90 b/src/Modules/log_mesh.f90 similarity index 100% rename from Libraries/Modules/log_mesh.f90 rename to src/Modules/log_mesh.f90 diff --git a/Libraries/Modules/machine b/src/Modules/machine similarity index 100% rename from Libraries/Modules/machine rename to src/Modules/machine diff --git a/Libraries/Modules/mat_mod.f90 b/src/Modules/mat_mod.f90 similarity index 100% rename from Libraries/Modules/mat_mod.f90 rename to src/Modules/mat_mod.f90 diff --git a/Libraries/Modules/matrix.f90 b/src/Modules/matrix.f90 similarity index 100% rename from Libraries/Modules/matrix.f90 rename to src/Modules/matrix.f90 diff --git a/Libraries/Modules/maxent.f90 b/src/Modules/maxent.f90 similarity index 100% rename from Libraries/Modules/maxent.f90 rename to src/Modules/maxent.f90 diff --git a/Libraries/Modules/maxent_stoch.G90 b/src/Modules/maxent_stoch.G90 similarity index 100% rename from Libraries/Modules/maxent_stoch.G90 rename to src/Modules/maxent_stoch.G90 diff --git a/Libraries/Modules/maxent_stoch.f90 b/src/Modules/maxent_stoch.f90 similarity index 100% rename from Libraries/Modules/maxent_stoch.f90 rename to src/Modules/maxent_stoch.f90 diff --git a/Libraries/Modules/maxent_stoch_w.f90 b/src/Modules/maxent_stoch_w.f90 similarity index 100% rename from Libraries/Modules/maxent_stoch_w.f90 rename to src/Modules/maxent_stoch_w.f90 diff --git a/Libraries/Modules/precdef.mod.f90 b/src/Modules/precdef.mod.f90 similarity index 100% rename from Libraries/Modules/precdef.mod.f90 rename to src/Modules/precdef.mod.f90 diff --git a/Libraries/Modules/smooth_stoch.f90 b/src/Modules/smooth_stoch.f90 similarity index 100% rename from Libraries/Modules/smooth_stoch.f90 rename to src/Modules/smooth_stoch.f90 diff --git a/Libraries/Modules/tmp.f90 b/src/Modules/tmp.f90 similarity index 100% rename from Libraries/Modules/tmp.f90 rename to src/Modules/tmp.f90 diff --git a/src/MyEis/CMakeLists.txt b/src/MyEis/CMakeLists.txt new file mode 100644 index 000000000..1c25eeaee --- /dev/null +++ b/src/MyEis/CMakeLists.txt @@ -0,0 +1,8 @@ +# Eispack library +SET(EIS_src ${SRCEIS}/balanc.f ${SRCEIS}/balbak.f ${SRCEIS}/cbabk2.f ${SRCEIS}/cbal.f ${SRCEIS}/cdiv.f ${SRCEIS}/cg.f ${SRCEIS}/ch.f +${SRCEIS}/comqr.f ${SRCEIS}/comqr2.f ${SRCEIS}/corth.f ${SRCEIS}/csroot.f ${SRCEIS}/elmhes.f ${SRCEIS}/eltran.f ${SRCEIS}/epslon.f ${SRCEIS}/hqr.f + ${SRCEIS}/hqr2.f ${SRCEIS}/htribk.f ${SRCEIS}/htridi.f ${SRCEIS}/pythag.f ${SRCEIS}/rg.f ${SRCEIS}/rs.f ${SRCEIS}/tql2.f + ${SRCEIS}/tqlrat.f ${SRCEIS}/tred1.f ${SRCEIS}/tred2.f) + + +ADD_LIBRARY(${MYEIS} STATIC ${EIS_src}) diff --git a/Libraries/MyEis/Makefile b/src/MyEis/Makefile similarity index 100% rename from Libraries/MyEis/Makefile rename to src/MyEis/Makefile diff --git a/Libraries/MyEis/balanc.f b/src/MyEis/balanc.f similarity index 100% rename from Libraries/MyEis/balanc.f rename to src/MyEis/balanc.f diff --git a/Libraries/MyEis/balbak.f b/src/MyEis/balbak.f similarity index 100% rename from Libraries/MyEis/balbak.f rename to src/MyEis/balbak.f diff --git a/Libraries/MyEis/cbabk2.f b/src/MyEis/cbabk2.f similarity index 100% rename from Libraries/MyEis/cbabk2.f rename to src/MyEis/cbabk2.f diff --git a/Libraries/MyEis/cbal.f b/src/MyEis/cbal.f similarity index 100% rename from Libraries/MyEis/cbal.f rename to src/MyEis/cbal.f diff --git a/Libraries/MyEis/cdiv.f b/src/MyEis/cdiv.f similarity index 100% rename from Libraries/MyEis/cdiv.f rename to src/MyEis/cdiv.f diff --git a/Libraries/MyEis/cg.f b/src/MyEis/cg.f similarity index 100% rename from Libraries/MyEis/cg.f rename to src/MyEis/cg.f diff --git a/Libraries/MyEis/ch.f b/src/MyEis/ch.f similarity index 100% rename from Libraries/MyEis/ch.f rename to src/MyEis/ch.f diff --git a/Libraries/MyEis/comp b/src/MyEis/comp similarity index 100% rename from Libraries/MyEis/comp rename to src/MyEis/comp diff --git a/Libraries/MyEis/comqr.f b/src/MyEis/comqr.f similarity index 100% rename from Libraries/MyEis/comqr.f rename to src/MyEis/comqr.f diff --git a/Libraries/MyEis/comqr2.f b/src/MyEis/comqr2.f similarity index 100% rename from Libraries/MyEis/comqr2.f rename to src/MyEis/comqr2.f diff --git a/Libraries/MyEis/corth.f b/src/MyEis/corth.f similarity index 100% rename from Libraries/MyEis/corth.f rename to src/MyEis/corth.f diff --git a/Libraries/MyEis/csroot.f b/src/MyEis/csroot.f similarity index 100% rename from Libraries/MyEis/csroot.f rename to src/MyEis/csroot.f diff --git a/Libraries/MyEis/elmhes.f b/src/MyEis/elmhes.f similarity index 100% rename from Libraries/MyEis/elmhes.f rename to src/MyEis/elmhes.f diff --git a/Libraries/MyEis/eltran.f b/src/MyEis/eltran.f similarity index 100% rename from Libraries/MyEis/eltran.f rename to src/MyEis/eltran.f diff --git a/Libraries/MyEis/epslon.f b/src/MyEis/epslon.f similarity index 100% rename from Libraries/MyEis/epslon.f rename to src/MyEis/epslon.f diff --git a/Libraries/MyEis/hqr.f b/src/MyEis/hqr.f similarity index 100% rename from Libraries/MyEis/hqr.f rename to src/MyEis/hqr.f diff --git a/Libraries/MyEis/hqr2.f b/src/MyEis/hqr2.f similarity index 100% rename from Libraries/MyEis/hqr2.f rename to src/MyEis/hqr2.f diff --git a/Libraries/MyEis/htribk.f b/src/MyEis/htribk.f similarity index 100% rename from Libraries/MyEis/htribk.f rename to src/MyEis/htribk.f diff --git a/Libraries/MyEis/htridi.f b/src/MyEis/htridi.f similarity index 100% rename from Libraries/MyEis/htridi.f rename to src/MyEis/htridi.f diff --git a/Libraries/MyEis/pythag.f b/src/MyEis/pythag.f similarity index 100% rename from Libraries/MyEis/pythag.f rename to src/MyEis/pythag.f diff --git a/Libraries/MyEis/rg.f b/src/MyEis/rg.f similarity index 100% rename from Libraries/MyEis/rg.f rename to src/MyEis/rg.f diff --git a/Libraries/MyEis/rs.f b/src/MyEis/rs.f similarity index 100% rename from Libraries/MyEis/rs.f rename to src/MyEis/rs.f diff --git a/Libraries/MyEis/tql2.f b/src/MyEis/tql2.f similarity index 100% rename from Libraries/MyEis/tql2.f rename to src/MyEis/tql2.f diff --git a/Libraries/MyEis/tqlrat.f b/src/MyEis/tqlrat.f similarity index 100% rename from Libraries/MyEis/tqlrat.f rename to src/MyEis/tqlrat.f diff --git a/Libraries/MyEis/tred1.f b/src/MyEis/tred1.f similarity index 100% rename from Libraries/MyEis/tred1.f rename to src/MyEis/tred1.f diff --git a/Libraries/MyEis/tred2.f b/src/MyEis/tred2.f similarity index 100% rename from Libraries/MyEis/tred2.f rename to src/MyEis/tred2.f diff --git a/src/MyLin/CMakeLists.txt b/src/MyLin/CMakeLists.txt new file mode 100644 index 000000000..07ef34f71 --- /dev/null +++ b/src/MyLin/CMakeLists.txt @@ -0,0 +1,12 @@ +# Eispack library +SET(LIN_src ${SRCLIN}/cgedi.f +${SRCLIN}/cgefa.f +${SRCLIN}/dgedi.f +${SRCLIN}/dgefa.f +${SRCLIN}/zgedi.f +${SRCLIN}/zgefa.f +${SRCLIN}/zqrdc.f +${SRCLIN}/zqrsl.f +) + +ADD_LIBRARY(${MYLIN} STATIC ${LIN_src}) diff --git a/Libraries/MyLin/Makefile b/src/MyLin/Makefile similarity index 100% rename from Libraries/MyLin/Makefile rename to src/MyLin/Makefile diff --git a/Libraries/MyLin/bidon b/src/MyLin/bidon similarity index 100% rename from Libraries/MyLin/bidon rename to src/MyLin/bidon diff --git a/Libraries/MyLin/cgedi.f b/src/MyLin/cgedi.f similarity index 100% rename from Libraries/MyLin/cgedi.f rename to src/MyLin/cgedi.f diff --git a/Libraries/MyLin/cgefa.f b/src/MyLin/cgefa.f similarity index 100% rename from Libraries/MyLin/cgefa.f rename to src/MyLin/cgefa.f diff --git a/Libraries/MyLin/dgedi.f b/src/MyLin/dgedi.f similarity index 100% rename from Libraries/MyLin/dgedi.f rename to src/MyLin/dgedi.f diff --git a/Libraries/MyLin/dgefa.f b/src/MyLin/dgefa.f similarity index 100% rename from Libraries/MyLin/dgefa.f rename to src/MyLin/dgefa.f diff --git a/Libraries/MyLin/work.pc b/src/MyLin/work.pc similarity index 100% rename from Libraries/MyLin/work.pc rename to src/MyLin/work.pc diff --git a/Libraries/MyLin/work.pcl b/src/MyLin/work.pcl similarity index 100% rename from Libraries/MyLin/work.pcl rename to src/MyLin/work.pcl diff --git a/Libraries/MyLin/zgedi.f b/src/MyLin/zgedi.f similarity index 100% rename from Libraries/MyLin/zgedi.f rename to src/MyLin/zgedi.f diff --git a/Libraries/MyLin/zgefa.f b/src/MyLin/zgefa.f similarity index 100% rename from Libraries/MyLin/zgefa.f rename to src/MyLin/zgefa.f diff --git a/Libraries/MyLin/zqrdc.f b/src/MyLin/zqrdc.f similarity index 100% rename from Libraries/MyLin/zqrdc.f rename to src/MyLin/zqrdc.f diff --git a/Libraries/MyLin/zqrsl.f b/src/MyLin/zqrsl.f similarity index 100% rename from Libraries/MyLin/zqrsl.f rename to src/MyLin/zqrsl.f diff --git a/src/MyNag/CMakeLists.txt b/src/MyNag/CMakeLists.txt new file mode 100644 index 000000000..7eb8c6eaa --- /dev/null +++ b/src/MyNag/CMakeLists.txt @@ -0,0 +1,26 @@ +# Eispack library +SET(NAG_src ${SRCNAG}/F01QCF.f +${SRCNAG}/F01QDF.f +${SRCNAG}/F01QEF.f +${SRCNAG}/F01RCF.f +${SRCNAG}/F01REF.f +${SRCNAG}/F06AAZ.f +${SRCNAG}/F06FBF.f +${SRCNAG}/F06FJF.f +${SRCNAG}/F06FRF.f +${SRCNAG}/F06HBF.f +${SRCNAG}/F06HRF.f +${SRCNAG}/F06QHF.f +${SRCNAG}/F06KJF.f +${SRCNAG}/F06THF.f +${SRCNAG}/P01ABF.f +${SRCNAG}/P01ABW.f +${SRCNAG}/P01ABZ.f +${SRCNAG}/P01ABY.f +${SRCNAG}/P01ACF.f +${SRCNAG}/X02AJF.f +${SRCNAG}/X04AAF.f +${SRCNAG}/X04BAF.f +) + +ADD_LIBRARY(${MYNAG} STATIC ${NAG_src}) diff --git a/Libraries/MyNag/F01QCF.f b/src/MyNag/F01QCF.f similarity index 100% rename from Libraries/MyNag/F01QCF.f rename to src/MyNag/F01QCF.f diff --git a/Libraries/MyNag/F01QDF.f b/src/MyNag/F01QDF.f similarity index 100% rename from Libraries/MyNag/F01QDF.f rename to src/MyNag/F01QDF.f diff --git a/Libraries/MyNag/F01QEF.f b/src/MyNag/F01QEF.f similarity index 100% rename from Libraries/MyNag/F01QEF.f rename to src/MyNag/F01QEF.f diff --git a/Libraries/MyNag/F01RCF.f b/src/MyNag/F01RCF.f similarity index 100% rename from Libraries/MyNag/F01RCF.f rename to src/MyNag/F01RCF.f diff --git a/Libraries/MyNag/F01REF.f b/src/MyNag/F01REF.f similarity index 100% rename from Libraries/MyNag/F01REF.f rename to src/MyNag/F01REF.f diff --git a/Libraries/MyNag/F06AAZ.f b/src/MyNag/F06AAZ.f similarity index 100% rename from Libraries/MyNag/F06AAZ.f rename to src/MyNag/F06AAZ.f diff --git a/Libraries/MyNag/F06FBF.f b/src/MyNag/F06FBF.f similarity index 100% rename from Libraries/MyNag/F06FBF.f rename to src/MyNag/F06FBF.f diff --git a/Libraries/MyNag/F06FJF.f b/src/MyNag/F06FJF.f similarity index 100% rename from Libraries/MyNag/F06FJF.f rename to src/MyNag/F06FJF.f diff --git a/Libraries/MyNag/F06FRF.f b/src/MyNag/F06FRF.f similarity index 100% rename from Libraries/MyNag/F06FRF.f rename to src/MyNag/F06FRF.f diff --git a/Libraries/MyNag/F06FRF.f~ b/src/MyNag/F06FRF.f~ similarity index 100% rename from Libraries/MyNag/F06FRF.f~ rename to src/MyNag/F06FRF.f~ diff --git a/Libraries/MyNag/F06HBF.f b/src/MyNag/F06HBF.f similarity index 100% rename from Libraries/MyNag/F06HBF.f rename to src/MyNag/F06HBF.f diff --git a/Libraries/MyNag/F06HRF.f b/src/MyNag/F06HRF.f similarity index 100% rename from Libraries/MyNag/F06HRF.f rename to src/MyNag/F06HRF.f diff --git a/Libraries/MyNag/F06KJF.f b/src/MyNag/F06KJF.f similarity index 100% rename from Libraries/MyNag/F06KJF.f rename to src/MyNag/F06KJF.f diff --git a/Libraries/MyNag/F06QHF.f b/src/MyNag/F06QHF.f similarity index 100% rename from Libraries/MyNag/F06QHF.f rename to src/MyNag/F06QHF.f diff --git a/Libraries/MyNag/F06THF.f b/src/MyNag/F06THF.f similarity index 100% rename from Libraries/MyNag/F06THF.f rename to src/MyNag/F06THF.f diff --git a/Libraries/MyNag/Makefile b/src/MyNag/Makefile similarity index 100% rename from Libraries/MyNag/Makefile rename to src/MyNag/Makefile diff --git a/Libraries/MyNag/P01ABF.f b/src/MyNag/P01ABF.f similarity index 100% rename from Libraries/MyNag/P01ABF.f rename to src/MyNag/P01ABF.f diff --git a/Libraries/MyNag/P01ABW.f b/src/MyNag/P01ABW.f similarity index 100% rename from Libraries/MyNag/P01ABW.f rename to src/MyNag/P01ABW.f diff --git a/Libraries/MyNag/P01ABY.f b/src/MyNag/P01ABY.f similarity index 100% rename from Libraries/MyNag/P01ABY.f rename to src/MyNag/P01ABY.f diff --git a/Libraries/MyNag/P01ABZ.f b/src/MyNag/P01ABZ.f similarity index 100% rename from Libraries/MyNag/P01ABZ.f rename to src/MyNag/P01ABZ.f diff --git a/Libraries/MyNag/P01ACF.f b/src/MyNag/P01ACF.f similarity index 100% rename from Libraries/MyNag/P01ACF.f rename to src/MyNag/P01ACF.f diff --git a/Libraries/MyNag/X02AJF.f b/src/MyNag/X02AJF.f similarity index 100% rename from Libraries/MyNag/X02AJF.f rename to src/MyNag/X02AJF.f diff --git a/Libraries/MyNag/X04AAF.f b/src/MyNag/X04AAF.f similarity index 100% rename from Libraries/MyNag/X04AAF.f rename to src/MyNag/X04AAF.f diff --git a/Libraries/MyNag/X04BAF.f b/src/MyNag/X04BAF.f similarity index 100% rename from Libraries/MyNag/X04BAF.f rename to src/MyNag/X04BAF.f diff --git a/Libraries/MyNag/comp b/src/MyNag/comp similarity index 100% rename from Libraries/MyNag/comp rename to src/MyNag/comp diff --git a/Libraries/MyNag/work.pc b/src/MyNag/work.pc similarity index 100% rename from Libraries/MyNag/work.pc rename to src/MyNag/work.pc diff --git a/Libraries/MyNag/work.pcl b/src/MyNag/work.pcl similarity index 100% rename from Libraries/MyNag/work.pcl rename to src/MyNag/work.pcl diff --git a/src/Prog/CMakeLists.txt b/src/Prog/CMakeLists.txt new file mode 100644 index 000000000..5cde6af0b --- /dev/null +++ b/src/Prog/CMakeLists.txt @@ -0,0 +1,68 @@ +######################################## +# Set up how to compile the source files +######################################## + +# Add the source files +SET(Common_src ${SRCPROG}/Hop_mod.f90 + ${SRCPROG}/Operator.f90 + ${SRCPROG}/UDV_WRAP.F90 + ${SRCPROG}/cgr1.f90 + ${SRCPROG}/cgr2.f90 + ${SRCPROG}/cgr2_1.f90 + ${SRCPROG}/cgr2_2.f90 + ${SRCPROG}/control_mod.F90 + ${SRCPROG}/gperp.F90 + ${SRCPROG}/inconfc.F90 + ${SRCPROG}/main.F90 + ${SRCPROG}/nranf.f90 + ${SRCPROG}/outconfc.F90 + ${SRCPROG}/print_bin_mod.F90 + ${SRCPROG}/tau_m.f90 + ${SRCPROG}/upgrade.f90 + ${SRCPROG}/wrapgrdo.f90 + ${SRCPROG}/wrapgrup.f90 + ${SRCPROG}/wrapul.f90 + ${SRCPROG}/wrapur.f90 +) + +##################################### +# Tell how to install this executable +##################################### + +IF(WIN32) + SET(CMAKE_INSTALL_PREFIX "C:\\Program Files") +ELSE() + SET(CMAKE_INSTALL_PREFIX /usr/local) +ENDIF(WIN32) +#INSTALL(TARGETS ${PROGEXE} RUNTIME DESTINATION bin) + +# Define the executables in terms of the source files +#Add in this loop further models +FOREACH(MODELTYPE Hub SPT Ising) +SET(${MODELTYPE}_src ${SRCPROG}/Hamiltonian_${MODELTYPE}.F90 ${Common_src}) +ADD_EXECUTABLE(${PROGEXE}_${MODELTYPE} ${${MODELTYPE}_src}) + +##################################################### +# Add the needed libraries and special compiler flags +##################################################### +TARGET_LINK_LIBRARIES(${PROGEXE}_${MODELTYPE} ${MODULES} ${MYEIS} ${MYNAG} ${MYLIN}) + +# Uncomment if you need to link to BLAS and LAPACK +TARGET_LINK_LIBRARIES(${PROGEXE}_${MODELTYPE} ${BLAS_LIBRARIES} + ${LAPACK_LIBRARIES} + ${CMAKE_THREAD_LIBS_INIT}) + +# Uncomment if you have parallization +IF(USE_OPENMP) + SET_TARGET_PROPERTIES(${PROGEXE}_${MODELTYPE} PROPERTIES + COMPILE_FLAGS "${OpenMP_Fortran_FLAGS}" + LINK_FLAGS "${OpenMP_Fortran_FLAGS}") +ELSEIF(USE_MPI) + SET_TARGET_PROPERTIES(${PROGEXE}_${MODELTYPE} PROPERTIES + COMPILE_FLAGS "${MPI_Fortran_COMPILE_FLAGS}" + LINK_FLAGS "${MPI_Fortran_LINK_FLAGS}") + INCLUDE_DIRECTORIES(${MPI_Fortran_INCLUDE_PATH}) + TARGET_LINK_LIBRARIES(${PROGEXE}_${MODELTYPE} ${MPI_Fortran_LIBRARIES}) +ENDIF(USE_OPENMP) +INSTALL(TARGETS ${PROGEXE}_${MODELTYPE} RUNTIME DESTINATION bin) +ENDFOREACH() diff --git a/Prog_7/Hamiltonian_Hub.f90 b/src/Prog/Hamiltonian_Hub.F90 similarity index 100% rename from Prog_7/Hamiltonian_Hub.f90 rename to src/Prog/Hamiltonian_Hub.F90 diff --git a/Prog_8/Hamiltonian_Ising.f90 b/src/Prog/Hamiltonian_Ising.F90 similarity index 100% rename from Prog_8/Hamiltonian_Ising.f90 rename to src/Prog/Hamiltonian_Ising.F90 diff --git a/Prog_8/Hamiltonian_SPT.f90 b/src/Prog/Hamiltonian_SPT.F90 similarity index 100% rename from Prog_8/Hamiltonian_SPT.f90 rename to src/Prog/Hamiltonian_SPT.F90 diff --git a/Prog_7/Hop_mod.f90 b/src/Prog/Hop_mod.f90 similarity index 100% rename from Prog_7/Hop_mod.f90 rename to src/Prog/Hop_mod.f90 diff --git a/Prog_8/Makefile b/src/Prog/Makefile similarity index 100% rename from Prog_8/Makefile rename to src/Prog/Makefile diff --git a/Prog_8/Operator.f90 b/src/Prog/Operator.f90 similarity index 98% rename from Prog_8/Operator.f90 rename to src/Prog/Operator.f90 index e2b1e9fa0..82dc0012e 100644 --- a/Prog_8/Operator.f90 +++ b/src/Prog/Operator.f90 @@ -287,7 +287,7 @@ Subroutine Op_Wrapup(Mat,Op,spin,Ndim,N_Type) ! Local Complex (Kind=8) :: VH(Ndim,Op%N), Z, Z1 - Integer :: n, i, m, m1 + Integer :: n, i, m, m1, oppm, oppn @@ -329,17 +329,19 @@ Subroutine Op_Wrapup(Mat,Op,spin,Ndim,N_Type) enddo enddo Do n = 1,Op%N + oppn = Op%P(n) Do I = 1,Ndim - Mat(Op%P(n),I) = VH(I,n) + Mat(oppn,I) = VH(I,n) Enddo Enddo elseif (N_Type == 2) then VH = cmplx(0.d0,0.d0) do n = 1,Op%N Do m = 1,Op%N - Z1 = conjg(Op%U(n,m)) + Z1 = conjg(Op%U(n,m)) + oppm = Op%P(m) DO I = 1,Ndim - VH(I,n) = VH(I,n) + Mat(I,Op%P(m)) * Z1 + VH(I,n) = VH(I,n) + Mat(I,oppm) * Z1 Enddo enddo Enddo diff --git a/Prog_7/UDV_WRAP.f90 b/src/Prog/UDV_WRAP.F90 similarity index 100% rename from Prog_7/UDV_WRAP.f90 rename to src/Prog/UDV_WRAP.F90 diff --git a/Prog_7/cgr1.f90 b/src/Prog/cgr1.f90 similarity index 100% rename from Prog_7/cgr1.f90 rename to src/Prog/cgr1.f90 diff --git a/Prog_7/cgr2.f90 b/src/Prog/cgr2.f90 similarity index 100% rename from Prog_7/cgr2.f90 rename to src/Prog/cgr2.f90 diff --git a/Prog_7/cgr2_1.f90 b/src/Prog/cgr2_1.f90 similarity index 100% rename from Prog_7/cgr2_1.f90 rename to src/Prog/cgr2_1.f90 diff --git a/Prog_7/cgr2_2.f90 b/src/Prog/cgr2_2.f90 similarity index 100% rename from Prog_7/cgr2_2.f90 rename to src/Prog/cgr2_2.f90 diff --git a/Prog_7/control_mod.f90 b/src/Prog/control_mod.F90 similarity index 100% rename from Prog_7/control_mod.f90 rename to src/Prog/control_mod.F90 diff --git a/Prog_7/gperp.f90 b/src/Prog/gperp.F90 similarity index 100% rename from Prog_7/gperp.f90 rename to src/Prog/gperp.F90 diff --git a/Prog_7/inconfc.f90 b/src/Prog/inconfc.F90 similarity index 100% rename from Prog_7/inconfc.f90 rename to src/Prog/inconfc.F90 diff --git a/Prog_7/machine b/src/Prog/machine similarity index 100% rename from Prog_7/machine rename to src/Prog/machine diff --git a/Prog_7/main.f90 b/src/Prog/main.F90 similarity index 100% rename from Prog_7/main.f90 rename to src/Prog/main.F90 diff --git a/Prog_7/nranf.f90 b/src/Prog/nranf.f90 similarity index 100% rename from Prog_7/nranf.f90 rename to src/Prog/nranf.f90 diff --git a/Prog_7/outconfc.f90 b/src/Prog/outconfc.F90 similarity index 100% rename from Prog_7/outconfc.f90 rename to src/Prog/outconfc.F90 diff --git a/Prog_8/print_bin_mod.f90 b/src/Prog/print_bin_mod.F90 similarity index 100% rename from Prog_8/print_bin_mod.f90 rename to src/Prog/print_bin_mod.F90 diff --git a/Prog_7/tau_m.f90 b/src/Prog/tau_m.f90 similarity index 100% rename from Prog_7/tau_m.f90 rename to src/Prog/tau_m.f90 diff --git a/Prog_8/upgrade.f90 b/src/Prog/upgrade.f90 similarity index 100% rename from Prog_8/upgrade.f90 rename to src/Prog/upgrade.f90 diff --git a/Prog_8/wrapgrdo.f90 b/src/Prog/wrapgrdo.f90 similarity index 100% rename from Prog_8/wrapgrdo.f90 rename to src/Prog/wrapgrdo.f90 diff --git a/Prog_8/wrapgrup.f90 b/src/Prog/wrapgrup.f90 similarity index 100% rename from Prog_8/wrapgrup.f90 rename to src/Prog/wrapgrup.f90 diff --git a/Prog_7/wrapul.f90 b/src/Prog/wrapul.f90 similarity index 100% rename from Prog_7/wrapul.f90 rename to src/Prog/wrapul.f90 diff --git a/Prog_7/wrapur.f90 b/src/Prog/wrapur.f90 similarity index 100% rename from Prog_7/wrapur.f90 rename to src/Prog/wrapur.f90