construct_gs_transform_matrix Subroutine

    subroutine construct_gs_transform_matrix(overlap, S, S_tilde, k, N, matrix_size)

        ! Construct the matrix, S, which transforms the Krylov vectors to a set of orthogonal vectors.

        ! We use the notation from the Appendix of Phys. Rev. B. 85, 205119 for the matrices and the
        ! indices, except we use 'm' instead of 'n' and 'l' instead of 'k' for the indices.

        ! Usage: overlap on input should contain the overlap matrix of the Krylov vectors that
        ! you want to produce a transformation matrix for. All other matrices input should be
        ! allocated on input, and should be the same size as overlap.

        real(dp), intent(inout) :: overlap(:, :), S(:, :), S_tilde(:, :), k(:, :), N(:)
        integer, intent(in) :: matrix_size
        integer :: m, i, l, p, q

        S = 0.0_dp
        S_tilde = 0.0_dp
        k = 0.0_dp
        N = 0.0_dp

        do m = 1, matrix_size
            S(m, m) = 1.0_dp
            S_tilde(m, m) = 1.0_dp
        end do

        N(1) = overlap(1, 1)
        S(1, 1) = 1 / sqrt(N(1))

        do m = 2, matrix_size
            do i = 1, m - 1
                do l = 1, i
                    k(i, m) = k(i, m) + S(l, i) * overlap(m, l)
                end do
            end do
            do p = 1, m - 1
                do q = p, m - 1
                    S_tilde(p, m) = S_tilde(p, m) - k(q, m) * S(p, q)
                end do
            end do
            do q = 1, m
                do p = 1, m
                    N(m) = N(m) + S_tilde(p, m) * S_tilde(q, m) * overlap(p, q)
                end do
            end do
            if (N(m) < 0.0_dp) return
            S(:, m) = S_tilde(:, m) / sqrt(N(m))
        end do

    end subroutine construct_gs_transform_matrix