nnls Subroutine

    subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,mode,max_iter)

    implicit none

    integer,intent(in)                      :: mda     !! first dimensioning parameter for the array `a`.
    integer,intent(in)                      :: n
    real(wp),dimension(mda,n),intent(inout) :: a       !! on entry, contains the `m` by `n`
                                                       !! matrix, `a`. on exit, contains
                                                       !! the product matrix, `q*a`, where `q` is an
                                                       !! `m` by `m` orthogonal matrix generated implicitly by
                                                       !! this subroutine.
    integer,intent(in)                      :: m
    real(wp),dimension(m),intent(inout)     :: b       !! on entry, contains the m-vector `b`. on exit, contains `q*b`.
    real(wp),dimension(n),intent(out)       :: x       !! the solution vector.
    real(wp),intent(out)                    :: rnorm   !! euclidean norm of the residual vector.
    real(wp),dimension(n),intent(inout)     :: w       !! array of working space.  on exit `w` will contain
                                                       !! the dual solution vector. `w` will satisfy `w(i) = 0`
                                                       !! for all `i` in set `p` and `w(i) <= 0` for all `i` in set `z`.
    real(wp),dimension(m),intent(inout)     :: zz      !! an m-array of working space.
    integer,intent(out)                     :: mode    !! this is a success-failure flag with the following meanings:
                                                       !!
                                                       !! * ***1*** the solution has been computed successfully.
                                                       !! * ***2*** the dimensions of the problem are bad. either `m<=0` or `n<=0`.
                                                       !! * ***3*** iteration count exceeded. more than `3*n` iterations.
    integer,intent(in)                      :: max_iter !! maximum number of iterations (if <=0, then `3*n` is used)

    integer :: i,ii,ip,iter,itmax,iz,iz1,iz2,izmax,j,jj,jz,l,npp1,nsetp,rtnkey
    real(wp) :: alpha,asave,cc,sm,ss,t,temp,unorm,up,wmax,ztest
    real(wp),dimension(1) :: dummy
    integer,dimension(n) :: index   !! an integer working array.
                                    !! the contents of this array define the sets
                                    !! `p` and `z` as follows:
                                    !!
                                    !! * `index(1:nsetp) = set p`.
                                    !! * `index(iz1:iz2) = set z`.
                                    !!
                                    !! where: `iz1 = nsetp + 1 = npp1`, `iz2 = n`

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

    mode = 1
    if ( m<=0 .or. n<=0 ) then
        mode = 2
        return
    end if
    iter = 0

    if (max_iter<=0) then
        itmax = 3*n
    else
        itmax = max_iter
    end if

    ! initialize the arrays index(1:n) and x(1:n).
    x = zero
    index = [(i, i=1,n)]
    iz2 = n
    iz1 = 1
    nsetp = 0
    npp1 = 1

    main : do

        ! ******  main loop begins here  ******
        ! quit if all coefficients are already in the solution.
        ! or if m cols of a have been triangularized.

        if ( iz1<=iz2 .and. nsetp<m ) then

            ! compute components of the dual (negative gradient) vector w().

            do iz = iz1 , iz2
                j = index(iz)
                sm = zero
                do l = npp1 , m
                    sm = sm + a(l,j)*b(l)
                end do
                w(j) = sm
            end do

            do
                ! find largest positive w(j).
                wmax = zero
                do iz = iz1 , iz2
                    j = index(iz)
                    if ( w(j)>wmax ) then
                        wmax = w(j)
                        izmax = iz
                    end if
                end do

                ! if wmax <= 0. go to termination.
                ! this indicates satisfaction of the kuhn-tucker conditions.
                if ( wmax<=zero ) then
                    call termination()
                    return
                end if

                iz = izmax
                j = index(iz)

                ! the sign of w(j) is ok for j to be moved to set p.
                ! begin the transformation and check new diagonal element to avoid
                ! near linear dependence.

                asave = a(npp1,j)
                call h12(1,npp1,npp1+1,m,a(1,j),1,up,dummy,1,1,0)
                unorm = zero
                if ( nsetp/=0 ) then
                    do l = 1 , nsetp
                        unorm = unorm + a(l,j)**2
                    end do
                end if
                unorm = sqrt(unorm)
                if (abs(a(npp1,j))*factor >= unorm*eps) then

                    ! col j is sufficiently independent.  copy b into zz, update zz
                    ! and solve for ztest ( = proposed new value for x(j) ).

                    do l = 1 , m
                        zz(l) = b(l)
                    end do
                    call h12(2,npp1,npp1+1,m,a(1,j),1,up,zz,1,1,1)
                    ztest = zz(npp1)/a(npp1,j)

                    ! see if ztest is positive
                    if ( ztest>zero ) then

                        ! the index j=index(iz) has been selected to be moved from
                        ! set z to set p. update b, update indices, apply householder
                        ! transformations to cols in new set z, zero subdiagonal elts in
                        ! col j, set w(j)=0.

                        do l = 1 , m
                            b(l) = zz(l)
                        end do

                        index(iz) = index(iz1)
                        index(iz1) = j
                        iz1 = iz1 + 1
                        nsetp = npp1
                        npp1 = npp1 + 1

                        if ( iz1<=iz2 ) then
                            do jz = iz1 , iz2
                                jj = index(jz)
                                call h12(2,nsetp,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
                            end do
                        end if

                        if ( nsetp/=m ) then
                            do l = npp1 , m
                                a(l,j) = zero
                            end do
                        end if

                        w(j) = zero
                        ! solve the triangular system.
                        ! store the solution temporarily in zz().
                        rtnkey = 1
                        exit
                    end if
                end if

                ! reject j as a candidate to be moved from set z to set p.
                ! restore a(npp1,j), set w(j)=0., and loop back to test dual
                ! coeffs again.

                a(npp1,j) = asave
                w(j) = zero

            end do

        else
            call termination()
            return
        end if

        ! ******  end of main loop  ******

        secondary : do

            ! the following block of code is used as an internal subroutine
            ! to solve the triangular system, putting the solution in zz().
            do l = 1 , nsetp
                ip = nsetp + 1 - l
                if ( l/=1 ) then
                    do ii = 1 , ip
                        zz(ii) = zz(ii) - a(ii,jj)*zz(ip+1)
                    end do
                end if
                jj = index(ip)
                zz(ip) = zz(ip)/a(ip,jj)
            end do
            if (rtnkey/=1 .and. rtnkey/=2) return

            ! ******  secondary loop begins here ******

            ! iteration counter.

            iter = iter + 1
            if ( iter>itmax ) then
                mode = 3
                !write (*,'(/a)') ' nnls quitting on iteration count.'
                call termination()
                return
            end if

            ! see if all new constrained coeffs are feasible.
            ! if not compute alpha.

            alpha = two
            do ip = 1 , nsetp
                l = index(ip)
                if ( zz(ip)<=zero ) then
                    t = -x(l)/(zz(ip)-x(l))
                    if ( alpha>t ) then
                        alpha = t
                        jj = ip
                    end if
                end if
            end do

            ! if all new constrained coeffs are feasible then alpha will
            ! still = 2.    if so exit from secondary loop to main loop.

            if ( abs(alpha-two)<=zero ) then
                ! ******  end of secondary loop  ******

                do ip = 1 , nsetp
                    i = index(ip)
                    x(i) = zz(ip)
                end do
                ! all new coeffs are positive.  loop back to beginning.
                cycle main
            end if

            ! otherwise use alpha which will be between 0. and 1. to
            ! interpolate between the old x and the new zz.

            do ip = 1 , nsetp
                l = index(ip)
                x(l) = x(l) + alpha*(zz(ip)-x(l))
            end do

            ! modify a and b and the index arrays to move coefficient i
            ! from set p to set z.

            i = index(jj)

            move_p : do
                x(i) = zero

                if ( jj/=nsetp ) then
                    jj = jj + 1
                    do j = jj , nsetp
                        ii = index(j)
                        index(j-1) = ii
                        call g1(a(j-1,ii),a(j,ii),cc,ss,a(j-1,ii))
                        a(j,ii) = zero
                        do l = 1 , n
                            if ( l/=ii ) then
                                ! apply procedure g2 (cc,ss,a(j-1,l),a(j,l))
                                temp = a(j-1,l)
                                a(j-1,l) = cc*temp + ss*a(j,l)
                                a(j,l) = -ss*temp + cc*a(j,l)
                            end if
                        end do
                        ! apply procedure g2 (cc,ss,b(j-1),b(j))
                        temp = b(j-1)
                        b(j-1) = cc*temp + ss*b(j)
                        b(j) = -ss*temp + cc*b(j)
                    end do
                end if

                npp1 = nsetp
                nsetp = nsetp - 1
                iz1 = iz1 - 1
                index(iz1) = i

                ! see if the remaining coeffs in set p are feasible.  they should
                ! be because of the way alpha was determined.
                ! if any are infeasible it is due to round-off error.  any
                ! that are nonpositive will be set to zero
                ! and moved from set p to set z.

                do jj = 1 , nsetp
                    i = index(jj)
                    if ( x(i)<=zero ) cycle move_p
                end do
                exit move_p

            end do move_p

            ! copy b( ) into zz( ).  then solve again and loop back.

            do i = 1 , m
                zz(i) = b(i)
            end do
            rtnkey = 2

        end do secondary

    end do main

    contains

    subroutine termination()
        !! come to here for termination.
        !! compute the norm of the final residual vector.

        sm = zero
        if ( npp1<=m ) then
            do i = npp1 , m
                sm = sm + b(i)**2
            end do
        else
            do j = 1 , n
                w(j) = zero
            end do
        end if
        rnorm = sqrt(sm)
    end subroutine termination

    end subroutine nnls