-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathget_gauge_order.f90
232 lines (185 loc) · 8.06 KB
/
get_gauge_order.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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
MODULE get_gauge_order
IMPLICIT NONE
PRIVATE
PUBLIC :: getgaugeorder
INTERFACE getgaugeorder
MODULE PROCEDURE getgaugeorder4
END INTERFACE
CONTAINS
SUBROUTINE getgaugeorder4( Ta, Tb, Dbond, Gauge, TrunError )
!
! Purpose:
! Contract tensor Ta and Tb and perform svd decomposition to get
! the unitary matrix U1 and U2.
! Ta and Tb are input tensors associated with sites.
! Dbond is the (expected) bond dimension of the tensor.
! Gauge is the output truncated unitary matrix.
! TrunError is the truncation error.
! Details definitions can be found in PRB 86 045139(2012)
! http://arxiv.org/abs/1201.1144
! By Yuzhi Liu on 04/26/2013; v1.
USE set_precisions
USE generic_contract
IMPLICIT NONE
! Data dictionary: declare calling parameter types & definitions
REAL(kind=DBL), INTENT(IN), DIMENSION(:,:,:,:) :: Ta ! Input array A
REAL(kind=DBL), INTENT(IN), DIMENSION(:,:,:,:) :: Tb ! Input array B
INTEGER, INTENT(IN) :: Dbond ! Bond dimension
REAL(kind=DBL), INTENT(OUT), DIMENSION(:,:,:) :: Gauge ! Output U matrix
REAL(kind=DBL), INTENT(OUT) :: TrunError ! Truncation error
! Data dictionary: declare local variable types & definitions
INTEGER :: dim1, dim2 ! Dim of Ta and Tb
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:,:,:) :: T1 ! Temp arrays
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:,:,:) :: T2 ! Temp arrays
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:) :: M1 ! Temp arrays
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:) :: U1 ! UL arrays
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:) :: U2 ! UR arrays
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:) :: Spectra1 ! A array
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:) :: Spectra2 ! A array
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:,:) :: V ! V arrays
REAL(kind=DBL) :: Tr1, Tr2 ! Traces of A
INTEGER :: err, istat ! error flag
dim1 = SIZE( Ta, 1 )
dim2 = SIZE( Tb, 1 )
!
! Calculate U2 ( right side)
!
! T1(2,4,2',4')
ALLOCATE( T1(SIZE(Ta,2), SIZE(Ta,4), SIZE(Ta,2), SIZE(Ta,4)), STAT=istat)
ifAllocatedT1: IF( ALLOCATED(T1) ) THEN
CALL TensorContract( Ta, Ta, 4, (/1,3/), 4, (/1,3/), error=err,outT=T1)
ELSE ifAllocatedT1
WRITE(*,*) 'Warning-Array T1 not allocated in getgaugeorder4.'
END IF ifAllocatedT1
! T2(6,4,6',4')
ALLOCATE( T2(SIZE(Tb,2), SIZE(Tb,3), SIZE(Tb,2), SIZE(Tb,3)), STAT=istat)
ifAllocatedT2: IF( ALLOCATED(T2) ) THEN
CALL TensorContract( Tb, Tb, 4, (/1,4/), 4, (/1,4/), error=err,outT=T2 )
ELSE ifAllocatedT2
WRITE(*,*) 'Warning-Array T2 not allocated in getgaugeorder4.'
END IF ifAllocatedT2
! M1(26,2'6')
ALLOCATE( M1(SIZE(T1,1)*SIZE(T2,1), SIZE(T1,3)*SIZE(T2,3)), STAT=istat )
ifAllocatedM1: IF( ALLOCATED(M1) ) THEN
CALL TensorContract( T1, T2, 4, (/2,4/), 4, (/2,4/), (/1,3,2,4/), &
(/2,2/), error=err, outT=M1 )
ELSE ifAllocatedM1
WRITE(*,*) 'Warning-Array M1 not allocated in getgaugeorder4.'
END IF ifAllocatedM1
DEALLOCATE(T1, STAT=istat)
DEALLOCATE(T2, STAT=istat)
! Reduce the potential round off error
M1 = ( M1 + TRANSPOSE(M1) ) / 2.0_DBL
!WRITE(*,*) 'The matrix M1 is: ', M1
! SVD decompose M1 to get U2
ALLOCATE( U2( SIZE(M1,1), SIZE(M1,1) ), STAT = istat )
ALLOCATE( V( SIZE(M1,1), SIZE(M1,1) ), STAT = istat )
ALLOCATE( Spectra2( SIZE(M1,1) ), STAT = istat )
CALL SVD( M1, U2, Spectra2, V, SIZE(M1,1), SIZE(M1,1) )
! Notice that the contents of M1 are destroyed.
!WRITE(*,*) 'The matrix M1 after SVD is: ', M1
DEALLOCATE( V, STAT = istat )
DEALLOCATE( M1, STAT = istat )
!WRITE(*,*) 'The matrix U2 is: ', U2
!WRITE(*,*) 'The singular values Spectra2: ', Spectra2
!
! Calculate U1 ( left side)
!
! T1(1,4,1',4')
ALLOCATE( T1(SIZE(Ta,2), SIZE(Ta,4), SIZE(Ta,2), SIZE(Ta,4)), STAT=istat)
IF( ALLOCATED(T1) ) THEN
CALL TensorContract( Ta, Ta, 4, (/2,3/), 4, (/2,3/), error=err,outT=T1)
ELSE
WRITE(*,*) 'Warning-Array T1 not allocated in getgaugeorder4.'
END IF
! T2(5,4,5',4')
ALLOCATE( T2(SIZE(Tb,2), SIZE(Tb,3), SIZE(Tb,2), SIZE(Tb,3)), STAT=istat)
IF( ALLOCATED(T2) ) THEN
CALL TensorContract( Tb, Tb, 4, (/2,4/), 4, (/2,4/), error=err,outT=T2 )
ELSE
WRITE(*,*) 'Warning-Array T2 not allocated in getgaugeorder4.'
END IF
! M1(15,1'5')
ALLOCATE( M1(SIZE(T1,1)*SIZE(T2,1), SIZE(T1,3)*SIZE(T2,3)), STAT=istat )
IF( ALLOCATED(M1) ) THEN
CALL TensorContract( T1, T2, 4, (/2,4/), 4, (/2,4/), (/1,3,2,4/), &
(/2,2/), error=err, outT=M1 )
ELSE
WRITE(*,*) 'Warning-Array M1 not allocated in getgaugeorder4.'
END IF
DEALLOCATE(T1, STAT=istat)
DEALLOCATE(T2, STAT=istat)
! Reduce the potential round off error
M1 = ( M1 + TRANSPOSE(M1) ) / 2.0_DBL
!WRITE(*,*) 'The matrix M1 is: ', M1
! SVD decompose M1 to get U2
ALLOCATE( U1( SIZE(M1,1), SIZE(M1,1) ), STAT = istat )
ALLOCATE( V( SIZE(M1,1), SIZE(M1,1) ), STAT = istat )
ALLOCATE( Spectra1( SIZE(M1,1) ), STAT = istat )
CALL SVD( M1, U1, Spectra1, V, SIZE(M1,1), SIZE(M1,1) )
! Notice that the contents of M1 are destroyed.
!WRITE(*,*) 'The matrix M1 after SVD is: ', M1
DEALLOCATE( V, STAT = istat )
DEALLOCATE( M1, STAT = istat )
!WRITE(*,*) 'The matrix U1 is: ', U1
!WRITE(*,*) 'The singular values Spectra1: ', Spectra1
! Calculate truncation error
Tr1 = SUM( Spectra1( Dbond+1: ) ) / SUM( Spectra1 )
Tr2 = SUM( Spectra2( Dbond+1: ) ) / SUM( Spectra2 )
TrunError = MIN( Tr1, Tr2 )
DEALLOCATE( Spectra1, STAT = istat )
DEALLOCATE( Spectra2, STAT = istat )
! Determine the Gauge
IF ( Tr1 <= Tr2 ) THEN
GAUGE = RESHAPE( U1(:, 1:Dbond), (/dim1, dim2, Dbond/) )
ELSE
GAUGE = RESHAPE( U2(:, 1:Dbond), (/dim1, dim2, Dbond/) )
END IF
DEALLOCATE( U1, STAT = istat )
DEALLOCATE( U2, STAT = istat )
END SUBROUTINE getgaugeorder4
SUBROUTINE SVD(A,U,S,V,M,N)
!
! Program computes the matrix singular value decomposition.
! Using Lapack library.
!
! Programmed by sukhbinder Singh
! 14th January 2011
! Modified by Yuzhi Liu on 04/26/2013.
!
USE set_precisions
IMPLICIT NONE
! Data dictionary: declare calling parameter types & definitions
INTEGER :: M, N
REAL(kind=DBL), INTENT(INOUT), DIMENSION(M,N) :: A
REAL(kind=DBL), INTENT(OUT), DIMENSION(M,M) :: U
REAL(kind=DBL), INTENT(OUT), DIMENSION(N,N) :: V
REAL(kind=DBL), INTENT(OUT), DIMENSION(MIN(M,N)) :: S
! Data dictionary: declare local variable types & definitions
REAL(kind=DBL), DIMENSION(N,N) :: VT
REAL(kind=DBL), ALLOCATABLE, DIMENSION(:) :: WORK
INTEGER :: LDA, LDU, LWORK, LDVT, INFO
CHARACTER :: JOBU, JOBVT
INTEGER :: istat
JOBU = 'A'
JOBVT = 'A'
LDA = M
LDU = M
LDVT = N
LWORK = MAX( 1, 3*MIN(M,N) + MAX(M,N), 5*MIN(M,N) )
ALLOCATE(WORK(LWORK), STAT=istat)
ifAllocatedWORK: IF ( ALLOCATED(WORK) ) THEN
CALL DGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, &
WORK, LWORK, INFO )
ifInfo: IF ( INFO == 0 ) THEN
!WRITE(*,*) 'SVD performed successfully!'
ELSE ifInfo
WRITE(*,*) 'SVD failed. Error INFO is : ', INFO
END IF ifInfo
ELSE ifAllocatedWORK
WRITE(*,*) 'Warning-Array WORK not allocated in SVD.'
END IF ifAllocatedWORK
DEALLOCATE(WORK, STAT = istat)
V = TRANSPOSE(VT)
END SUBROUTINE SVD
END MODULE get_gauge_order