hfti Subroutine

    subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g)

    implicit none

    integer,intent(in)                       :: mda   !! the first dimensioning parameter of matrix `a` (mda >= m).
    integer,intent(in)                       :: m
    integer,intent(in)                       :: n
    integer,intent(in)                       :: mdb   !! first dimensioning parameter of matrix `b` (mdb>=max(m,n))
    integer,intent(in)                       :: nb
    real(wp),dimension(mda,n),intent(inout)  :: a     !! the array `a` initially contains the \( m \times n \) matrix \(\mathbf{A}\)
                                                      !! of the least squares problem \( \mathbf{A} \mathbf{x} = \mathbf{b} \).
                                                      !! either `m >= n` or `m < n` is permitted.
                                                      !! there is no restriction on the rank of `a`.
                                                      !! the matrix `a` will be modified by the subroutine.
    real(wp),intent(in)                      :: tau   !! absolute tolerance parameter for pseudorank
                                                      !! determination, provided by the user.
    integer,intent(out)                      :: krank !! pseudorank of `a`, set by the subroutine.
    real(wp),dimension(nb),intent(out)       :: rnorm !! on exit, `rnorm(j)` will contain the euclidian
                                                      !! norm of the residual vector for the problem
                                                      !! defined by the `j-th` column vector of the array `b`.
    real(wp),dimension(n),intent(inout)      :: h     !! array of working space
    real(wp),dimension(n),intent(inout)      :: g     !! array of working space
    real(wp),dimension(mdb,nb),intent(inout) :: b     !! if `nb = 0` the subroutine will make no reference
                                                      !! to the array `b`. if `nb > 0` the array `b` must
                                                      !! initially contain the `m x nb` matrix `b` of the
                                                      !! the least squares problem `ax = b` and on return
                                                      !! the array `b` will contain the `n x nb` solution `x`.

    integer  :: i, ii, ip1, j, jb, jj, k, kp1, l, ldiag, lmax
    real(wp) :: hmax, sm, tmp
    logical  :: need_lmax
    integer,dimension(n) :: ip    !! integer array of working space
                                  !! recording permutation indices of column vectors

    real(wp),parameter :: factor = 0.001_wp

    k = 0
    ldiag = min(m,n)
    if ( ldiag<=0 ) then

        ! the solution vectors, x, are now
        ! in the first  n  rows of the array b(,).
        krank = k
        return

    else

        do j = 1 , ldiag

            need_lmax = .true.
            if ( j/=1 ) then
                ! update squared column lengths and find lmax
                lmax = j
                do l = j , n
                    h(l) = h(l) - a(j-1,l)**2
                    if ( h(l)>h(lmax) ) lmax = l
                end do
                if (factor*h(lmax) >= hmax*eps) need_lmax = .false.

            end if

            if (need_lmax) then
                ! compute squared column lengths and find lmax
                lmax = j
                do l = j , n
                    h(l) = zero
                    do i = j , m
                        h(l) = h(l) + a(i,l)**2
                    end do
                    if ( h(l)>h(lmax) ) lmax = l
                end do
                hmax = h(lmax)
            end if

            ! lmax has been determined

            ! do column interchanges if needed.

            ip(j) = lmax
            if ( ip(j)/=j ) then
                do i = 1 , m
                    tmp = a(i,j)
                    a(i,j) = a(i,lmax)
                    a(i,lmax) = tmp
                end do
                h(lmax) = h(j)
            end if

            ! compute the j-th transformation and apply it to a and b.

            call h12(1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
            call h12(2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb)

        end do

        ! determine the pseudorank, k, using the tolerance, tau.
        do j=1,ldiag
            if (abs(a(j,j))<=tau) exit
        end do
        k=j-1
        kp1=j

    end if

    ! compute the norms of the residual vectors.

    if ( nb>0 ) then
        do jb = 1 , nb
            tmp = zero
            if ( kp1<=m ) then
                do i = kp1 , m
                    tmp = tmp + b(i,jb)**2
                end do
            end if
            rnorm(jb) = sqrt(tmp)
        end do
    end if
    ! special for pseudorank = 0
    if ( k>0 ) then

        ! if the pseudorank is less than n compute householder
        ! decomposition of first k rows.

        if ( k/=n ) then
            do ii = 1 , k
                i = kp1 - ii
                call h12(1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1)
            end do
        end if

        if ( nb>0 ) then
            do jb = 1 , nb

                ! solve the k by k triangular system.

                do l = 1 , k
                    sm = zero
                    i = kp1 - l
                    if ( i/=k ) then
                        ip1 = i + 1
                        do j = ip1 , k
                            sm = sm + a(i,j)*b(j,jb)
                        end do
                    end if
                    b(i,jb) = (b(i,jb)-sm)/a(i,i)
                end do

                ! complete computation of solution vector.

                if ( k/=n ) then
                    do j = kp1 , n
                        b(j,jb) = zero
                    end do
                    do i = 1 , k
                        call h12(2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1)
                    end do
                end if

                ! re-order the solution vector to compensate for the
                ! column interchanges.

                do jj = 1 , ldiag
                    j = ldiag + 1 - jj
                    if ( ip(j)/=j ) then
                        l = ip(j)
                        tmp = b(l,jb)
                        b(l,jb) = b(j,jb)
                        b(j,jb) = tmp
                    end if
                end do
            end do
        end if

    else if ( nb>0 ) then

        do jb = 1 , nb
            do i = 1 , n
                b(i,jb) = zero
            end do
        end do

    end if
    krank = k

    end subroutine hfti