QR_Householder_decomposition Subroutine

    subroutine QR_Householder_decomposition(A, Q, R)
        real(dp), dimension(:, :), intent(in) :: A
        real(dp), dimension(size(A, 1), size(A, 2)), intent(out) :: Q, R
        real(dp), dimension(size(A, 1), size(A, 2)) :: Id, H, v_mat_tmp
        real(dp), dimension(size(A, 1)) :: v, u, x
        integer :: N, i, j, k
        real(dp) :: alpha, w, signe, norm_u

        N = size(A, 1)

        R = A

        Id = Identity_n(N)
        Q = Identity_n(N)

        do k = 1, N

            x = 0.d0
            u = 0.d0
            v = 0.d0
            v_mat_tmp = 0.d0
            x(k:N) = R(K:N, K)

            alpha = norm2(R(k:N, k))

            signe = -sign(alpha, x(k))
            u(k:N) = x(k:N) - signe * Id(k:N, k)

            norm_u = norm2(u)
            if (norm_u < TOL_CONVERGENCE_dp) cycle
            v(k:N) = u(k:N) / norm_u

            w = 1.d0
            do i = k, N
                do j = k, N
                    v_mat_tmp(i, j) = v(i) * v(j)
                end do
            end do

            H = Id
            H(k:N, k:N) = Id(k:N, k:N) - (1.d0 + w) * v_mat_tmp(k:N, k:N)

            Q = matmul(Q, H)

            R(k:N, k:N) = matmul(H(k:N, k:N), R(k:N, k:N))

        end do
    end subroutine QR_Householder_decomposition