-
Notifications
You must be signed in to change notification settings - Fork 0
/
distance_norm_ip.f90
156 lines (130 loc) · 7.83 KB
/
distance_norm_ip.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
!###################################################################################################################################
! VECTOR NORMS & INNER PRODUCTS
!###################################################################################################################################
complex(8) function inner_prod(d, psi, phi) ! Returns the Euclidean inner product of two complex VECTORS psi and phi
implicit none
integer :: d ! Dimension of the vectors
complex(8) :: psi(1:d), phi(1:d) ! Complex vectors whose inner product is to be computed
integer :: j ! Auxiliary variable
inner_prod = (0.d0,0.d0) ; do j = 1, d ; inner_prod = inner_prod + conjg(psi(j))*phi(j) ; enddo
! complexity: ~ O(16d**2), for large d
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function inner_prod_r(d, psi, phi) ! Returns the Euclidean inner product of two real VECTORS psi and phi
implicit none
integer :: d ! Dimension of the vectors
real(8) :: psi(1:d), phi(1:d) ! Real vectors whose inner product is to be computed
integer :: j ! Auxiliary variable
inner_prod_r = 0.d0 ; do j = 1, d ; inner_prod_r = inner_prod_r + psi(j)*phi(j) ; enddo
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function norm(d, psi) ! Returns the Euclidean norm of the complex VECTOR psi
implicit none
integer :: d ! Dimension of the vector
complex(8) :: psi(1:d) ! Complex vector whose norm is to be computed
complex(8) :: inner_prod ! Inner product function to be used
norm = sqrt(dble(inner_prod(d, psi, psi)))
! complexity: ~ O(16d**2), for large d
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function norm_r(d, psi) ! Returns the Euclidean norm of the real VECTOR psi
implicit none
integer :: d ! Dimension of the vector
real(8) :: psi(1:d) ! Real vector whose norm is to be computed
real(8) :: inner_prod_r ! Inner product function to be used
norm_r = sqrt( inner_prod_r(d, psi, psi) )
end
!###################################################################################################################################
! FIDELITIES
!###################################################################################################################################
real(8) function fidelity_pp(d, psi, phi) ! Returns the fidelity between two pure states
implicit none
integer :: d ! Dimension of the vectors
complex(8) :: psi(1:d), phi(1:d) ! Complex vectors whose fidelity is to be computed
complex(8) :: inner_prod, ip ! Inner product function
ip = inner_prod(d, psi, phi) ; fidelity_pp = (dreal(ip))**2.d0 + (dimag(ip))**2.d0
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function fidelity_pm(d, psi, rho) ! Returns the fidelity between a pure and a mixed state
implicit none
integer :: d ! Dimension of the vector and matrix
complex(8) :: psi(1:d) ! Pure state vector
complex(8) :: rho(1:d,1:d) ! Density matrix
complex(8) :: fid ! Auxiliary variable
integer :: j, k ! Auxiliary variables for counters
fid = 0.d0
do j = 1, d ; do k = 1, d ; fid = fid + conjg(psi(j))*rho(j,k)*psi(k) ; enddo ; enddo
fidelity_pm = dble(fid)
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function fidelity_mm(d, rho, zeta) ! Returns the fidelity between two mixed states
implicit none
integer :: d ! Dimension of the density matrices
complex(8) :: rho(1:d,1:d), zeta(1:d,1:d) ! Density matrices
real(8) :: r(1:d), z(1:d) ! Vectors for the eigenvalues of rho and zeta, respectively.
complex(8) :: A(1:d,1:d) ! Auxiliary matrix. We use r also for the eigenvalues if A=sqrt(rho)*zeta*sqrt(rho)
integer :: j, k, l ! Auxiliary variables for counters
complex(8) :: op(1:d,1:d) ! For the outer product
complex(8) :: inner_prod ! For the inner product function
call lapack_zheevd('V', d, rho, r) ; call lapack_zheevd('V', d, zeta, z)
A = 0.d0
do j = 1, d ; do k = 1, d ; do l = 1, d
call outer_product(d, rho(:,j), rho(:,l), op)
A = A + sqrt(r(j)*r(l))*z(k)*inner_prod(d, rho(:,j), zeta(:,k))*inner_prod(d, zeta(:,k), rho(:,l))*op
enddo ; enddo ; enddo
call lapack_zheevd('N', d, A, r) ; fidelity_mm = 0.d0 ; do j = 1, d ; fidelity_mm = fidelity_mm + sqrt(r(j)) ; enddo
fidelity_mm = (fidelity_mm)**2.d0
end
!###################################################################################################################################
! NORMS
!###################################################################################################################################
real(8) function norm_hs(m, n, A) ! Returns the Hilbert-Schmidt norm (or 2-norm) of a complex matrix A
implicit none
integer :: m, n ! Dimensions of the matrix A
complex(8) :: A(1:m,1:n) ! Matrix of which want to compute the 2-norm
integer :: j, k ! Auxiliary variables for counters
norm_hs = 0.d0
do j = 1, m ; do k = 1, n
norm_hs = norm_hs + (dble(A(j,k)))**2.d0 + (aimag(A(j,k)))**2.d0
enddo ; enddo
norm_hs = sqrt(norm_hs)
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function norm_tr(d, A) ! Returns the trace norm (or 1-norm) of an hermitian matrix A
! ||A||_1 = \sum_j |a_j|, where a_j are the eigenvalues of A
implicit none
integer :: d ! Dimension of the matrix A
complex(8) :: A(1:d,1:d) ! Matrix of which want to compute the 1-norm
real(8) :: egv(1:d) ! Eigenvalues of A
integer :: j ! Auxiliary variable for counters
call lapack_zheevd('N', d, A, egv) ; norm_tr = 0.d0 ; do j = 1, d ; norm_tr = norm_tr + abs(egv(j)) ; enddo
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function norm_tr_ge(d, A) ! Returns the trace norm (or 1-norm) of a general normal matrix A
! ||A||_1 = \sum_j |a_j|, where a_j are the eigenvalues of A
implicit none
integer :: d ! Dimension of the matrix A
complex(8) :: A(1:d,1:d) ! Matrix of which want to compute the 1-norm
complex(8) :: egv(1:d) ! Eigenvalues of A
integer :: j ! Auxiliary variable for counters
call lapack_zgeev('N', d, A, egv) ; norm_tr_ge = 0.d0 ; do j = 1, d ; norm_tr_ge = norm_tr_ge + abs(egv(j)) ; enddo
end
!------------------------------------------------------------------------------------------------------------------------------------
real(8) function norm_l1(m, n, A) ! Returns the l1-norm of a complex mxn complex matrix A
implicit none
integer :: m, n ! Dimensions of the matrix A
complex(8) :: A(1:m,1:n) ! The matrix
integer :: j, k
norm_l1 = 0.d0 ; do j=1,m ; do k=1,n ; norm_l1 = norm_l1 + abs(A(j,k)) ; enddo ; enddo
end
!###################################################################################################################################
! DISTANCES
!###################################################################################################################################
real(8) function d_hs(d, rho, zeta) ! Returns the Hilbert-Schmidt (2-norm) distance between two density matrices
implicit none
integer :: d ! Dimension of the density matrices
complex(8) :: rho(1:d,1:d), zeta(1:d,1:d) ! Density matrices
real(8) :: norm_hs ! Function for the Hilbert-Schmidt norm
d_hs = norm_hs(d, d, rho - zeta)
end
!###################################################################################################################################