-
Notifications
You must be signed in to change notification settings - Fork 3
/
ED_HAMILTONIAN_DIRECT_HxV.f90
176 lines (160 loc) · 6.07 KB
/
ED_HAMILTONIAN_DIRECT_HxV.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
! > SPARSE MAT-VEC DIRECT ON-THE-FLY PRODUCT
MODULE ED_HAMILTONIAN_DIRECT_HxV
USE ED_HAMILTONIAN_COMMON
USE ED_AUX_FUNX
implicit none
private
integer :: iiup,iidw
integer :: iud,jj
integer :: i,iup,idw
integer :: j,jup,jdw
integer :: m,mup,mdw
integer :: ishift
integer :: isector,jsector
integer :: ms
integer :: impi
integer :: ilat,jlat,iorb,jorb,ispin,jspin,is,js,ibath
integer :: k1,k2,k3,k4
integer :: ialfa,ibeta
real(8) :: sg1,sg2,sg3,sg4
complex(8) :: htmp,htmpup,htmpdw
logical :: Jcondition
!>Sparse Mat-Vec direct on-the-fly product
public :: directMatVec_main
#ifdef _MPI
public :: directMatVec_MPI_main
#endif
contains
subroutine directMatVec_main(Nloc,vin,Hv)
integer :: Nloc
complex(8),dimension(Nloc) :: vin
complex(8),dimension(Nloc) :: Hv
complex(8),dimension(:),allocatable :: vt,Hvt
integer,dimension(Ns) :: ibup,ibdw
integer,dimension(2*Ns_Ud) :: Indices,Jndices ![2-2*Norb]
integer,dimension(Ns_Ud,Ns_Orb) :: Nups,Ndws ![1,Ns]-[Norb,1+Nbath]
integer,dimension(Nlat,Norb) :: Nup,Ndw
real(8),dimension(Nlat,Nspin,Norb,Nbath) :: diag_hybr
real(8),dimension(Nlat,Nspin,Norb,Nbath) :: bath_diag
complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Nbath) :: Hbath_reconstructed
!
if(.not.Hstatus)stop "directMatVec_cc ERROR: Hsector NOT set"
isector=Hsector
!
if(Nloc/=getdim(isector))stop "directMatVec_cc ERROR: Nloc != dim(isector)"
!
!Get diagonal hybridization, bath energy
diag_hybr=zero
bath_diag=zero
do ibath=1,Nbath
Hbath_reconstructed(:,:,:,:,:,:,ibath)=Hbath_build(dmft_bath%item(ibath)%lambda)
do ilat=1,Nlat
do ispin=1,Nspin
do iorb=1,Norb
diag_hybr(ilat,ispin,iorb,ibath)=dmft_bath%item(ibath)%v(index_stride_lso(ilat,ispin,iorb))
bath_diag(ilat,ispin,iorb,ibath)=DREAL(Hbath_Reconstructed(ilat,ilat,ispin,ispin,iorb,iorb,ibath))
enddo
enddo
enddo
enddo
!
Hv=zero
!
!-----------------------------------------------!
!LOCAL HAMILTONIAN PART: H_loc*vin = vout
include "ED_HAMILTONIAN/direct/HxV_local.f90"
!
!UP HAMILTONIAN TERMS
include "ED_HAMILTONIAN/direct/HxV_up.f90"
!
!DW HAMILTONIAN TERMS
include "ED_HAMILTONIAN/direct/HxV_dw.f90"
!
!NON-LOCAL HAMILTONIAN PART: H_non_loc*vin = vout
if(Jhflag)then
include "ED_HAMILTONIAN/direct/HxV_non_local.f90"
endif
!-----------------------------------------------!
!
return
end subroutine directMatVec_main
#ifdef _MPI
subroutine directMatVec_MPI_main(Nloc,vin,Hv)
integer :: Nloc,N
complex(8),dimension(Nloc) :: Vin
complex(8),dimension(Nloc) :: Hv
complex(8),dimension(:),allocatable :: vt,Hvt
integer,dimension(Ns) :: ibup,ibdw
integer,dimension(2*Ns_Ud) :: Indices,Jndices ![2-2*Norb]
integer,dimension(Ns_Ud,Ns_Orb) :: Nups,Ndws ![1,Ns]-[Norb,1+Nbath]
integer,dimension(Nlat,Norb) :: Nup,Ndw
real(8),dimension(Nlat,Nspin,Norb,Nbath) :: diag_hybr
real(8),dimension(Nlat,Nspin,Norb,Nbath) :: bath_diag
complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Nbath) :: Hbath_reconstructed
!
integer :: MpiIerr
integer,allocatable,dimension(:) :: Counts
integer,allocatable,dimension(:) :: Offset
!
if(.not.Hstatus)stop "directMatVec_cc ERROR: Hsector NOT set"
isector=Hsector
!
!
if(MpiComm==MPI_UNDEFINED.OR.MpiComm==Mpi_Comm_Null)&
stop "directMatVec_MPI_cc ERRROR: MpiComm = MPI_UNDEFINED"
! if(.not.MpiStatus)stop "directMatVec_MPI_cc ERROR: MpiStatus = F"
!
!Get diagonal hybridization, bath energy
diag_hybr=zero
bath_diag=zero
do ibath=1,Nbath
Hbath_reconstructed(:,:,:,:,:,:,ibath)=Hbath_build(dmft_bath%item(ibath)%lambda)
do ilat=1,Nlat
do ispin=1,Nspin
do iorb=1,Norb
diag_hybr(ilat,ispin,iorb,ibath)=dmft_bath%item(ibath)%v(index_stride_lso(ilat,ispin,iorb))
bath_diag(ilat,ispin,iorb,ibath)=Hbath_Reconstructed(ilat,ilat,ispin,ispin,iorb,iorb,ibath)
enddo
enddo
enddo
enddo
!
Hv=zero
!
!-----------------------------------------------!
!LOCAL HAMILTONIAN PART: H_loc*vin = vout
include "ED_HAMILTONIAN/direct_mpi/HxV_local.f90"
!
!UP HAMILTONIAN TERMS: MEMORY CONTIGUOUS
include "ED_HAMILTONIAN/direct_mpi/HxV_up.f90"
!
!DW HAMILTONIAN TERMS: MEMORY NON-CONTIGUOUS
mpiQup=DimUp/MpiSize
if(MpiRank<mod(DimUp,MpiSize))MpiQup=MpiQup+1
allocate(vt(mpiQup*DimDw)) ;vt=zero
allocate(Hvt(mpiQup*DimDw));Hvt=zero
call vector_transpose_MPI(DimUp,MpiQdw,Vin,DimDw,MpiQup,vt) !Vin^T --> Vt
include "ED_HAMILTONIAN/direct_mpi/HxV_dw.f90"
deallocate(vt) ; allocate(vt(DimUp*mpiQdw)) ;vt=zero !reallocate Vt
call vector_transpose_MPI(DimDw,mpiQup,Hvt,DimUp,mpiQdw,vt) !Hvt^T --> Vt
Hv = Hv + Vt
deallocate(vt)
!
!NON-LOCAL HAMILTONIAN PART: H_non_loc*vin = vout
if(Jhflag)then
N = 0
call AllReduce_MPI(MpiComm,Nloc,N)
!
allocate(vt(N)) ; vt = zero
call allgather_vector_MPI(MpiComm,vin,vt)
!
include "ED_HAMILTONIAN/direct_mpi/HxV_non_local.f90"
!
deallocate(Vt)
endif
!-----------------------------------------------!
!
return
end subroutine directMatVec_MPI_main
#endif
END MODULE ED_HAMILTONIAN_DIRECT_HXV