subroutine Check_XY_orthogonality(X, Y)
implicit none
real(dp), intent(in) :: X(ov_space, ov_space), Y(ov_space, ov_space)
integer :: i, i_p, m, m_p, v, mi_ind, mp_ip_ind, mu, mu_p, j
real(dp) :: temp
character(len=*), parameter :: t_r = "Check_XY_orthogonality"
!Eigenvectors orthonormal
do mu = 1, ov_space
do mu_p = 1, ov_space
temp = 0.0_dp
do i = 1, ov_space
temp = temp + X(i, mu) * X(i, mu_p) - Y(i, mu) * Y(i, mu_p)
end do
if ((mu == mu_p) .and. ((abs(temp) - 1.0_dp) > 1.0e-6_dp)) then
write(stdout, *) mu, mu_p, temp
call writevector(X(:, mu), 'X(:,mu)')
call writevector(Y(:, mu), 'Y(:,mu)')
call stop_all(t_r, 'X/Y not normalized')
else if ((mu /= mu_p) .and. (abs(temp) > 1.0e-6_dp)) then
write(stdout, *) mu, mu_p
call stop_all(t_r, 'X/Y not orthogonal')
end if
end do
end do
!Rows orthonormal too
do i = 1, ov_space
do j = 1, ov_space
temp = 0.0_dp
do mu = 1, ov_space
temp = temp + X(i, mu) * X(j, mu) - Y(i, mu) * Y(j, mu)
end do
if ((i /= j) .and. (abs(temp) > 1.0e-7_dp)) then
write(stdout, *) i, j, temp
call stop_all(t_r, 'X/Y rows not orthogonal')
else if ((i == j) .and. (abs(temp) - 1.0_dp) > 1.0e-7_dp) then
write(stdout, *) i, j, temp
call stop_all(t_r, 'X/Y rows not normalized')
end if
end do
end do
end subroutine Check_XY_orthogonality