qform Subroutine

    subroutine qform(m, n, q, Ldq, Wa)
        implicit none

        integer, intent(in) :: m !! a positive integer input variable set to the number
                                !! of rows of a and the order of q.
        integer, intent(in) :: n !! a positive integer input variable set to the number
                                !! of columns of a.
        integer, intent(in) :: Ldq !! a positive integer input variable not less than m
                                  !! which specifies the leading dimension of the array q.
        real(wp), intent(inout) :: q(Ldq, m) !! an m by m array. on input the full lower trapezoid in
                                            !! the first min(m,n) columns of q contains the factored form.
                                            !! on output q has been accumulated into a square matrix.
        real(wp), intent(inout) :: Wa(m) !! a work array of length m.

        integer :: i, j, jm1, k, l, minmn, np1
        real(wp) :: sum, temp

        ! zero out upper triangle of q in the first min(m,n) columns.

        minmn = min(m, n)
        if (minmn >= 2) then
            do j = 2, minmn
                jm1 = j - 1
                do i = 1, jm1
                    q(i, j) = zero
                end do
            end do
        end if

        ! initialize remaining columns to those of the identity matrix.

        np1 = n + 1
        if (m >= np1) then
            do j = np1, m
                do i = 1, m
                    q(i, j) = zero
                end do
                q(j, j) = one
            end do
        end if

        ! accumulate q from its factored form.

        do l = 1, minmn
            k = minmn - l + 1
            do i = k, m
                Wa(i) = q(i, k)
                q(i, k) = zero
            end do
            q(k, k) = one
            if (Wa(k) /= zero) then
                do j = k, m
                    sum = zero
                    do i = k, m
                        sum = sum + q(i, j)*Wa(i)
                    end do
                    temp = sum/Wa(k)
                    do i = k, m
                        q(i, j) = q(i, j) - temp*Wa(i)
                    end do
                end do
            end if
        end do

    end subroutine qform