bvls Subroutine

    subroutine bvls (a, b, bnd, x, rnorm, nsetp, w, index, ierr, max_iter)

    implicit none

    real(wp),dimension(:,:),intent(inout) :: a  !! On entry A() contains the M by N matrix, A.
                                                !! On return A() contains the product matrix, Q*A, where
                                                !! Q is an M by M orthogonal matrix generated by this
                                                !! subroutine.  The dimensions are M=size(A,1) and N=size(A,2).
    real(wp),dimension(:),intent(inout)   :: b  !! On entry B() contains the M-vector, B.
                                                !! On return, B() contains Q*B.  The same Q multiplies A.
    real(wp),dimension(:,:),intent(in) :: bnd   !! BND(1,J) is the lower bound for X(J).
                                                !! BND(2,J) is the upper bound for X(J).
                                                !!
                                                !! Require:  BND(1,J)  <=  BND(2,J).
                                                !!
                                                !! The values BND(1,J) = -huge(ONE) and BND(2,J) = huge(ONE) are
                                                !! suggested choices to designate that there is no constraint in that
                                                !! direction.  The parameter ONE is 1.0 in the working precision.
    real(wp),dimension(:),intent(out)   :: x    !! On entry X() need not be initialized.  On return,
                                                !! X() will contain the solution N-vector.
    real(wp),intent(out)                :: rnorm !! Euclidean norm of the residual vector, b - A*X.
    integer,intent(out)                 :: nsetp !! Indicates the number of components of the solution
                                                 !! vector, X(), that are not at their constraint values.
    real(wp),dimension(:),intent(out)   :: w    !! An N-array.  On return, W() will contain the dual solution
                                                !! vector.   Using Set definitions below:
                                                !!
                                                !! * W(J) = 0 for all j in Set P,
                                                !! * W(J)  <=  0 for all j in Set Z, such that X(J) is at its
                                                !!   lower bound, and
                                                !! * W(J)  >=  0 for all j in Set Z, such that X(J) is at its
                                                !!   upper bound.
                                                !!
                                                !! If BND(1,J) = BND(2,J), so the variable X(J) is fixed,
                                                !! then W(J) will have an arbitrary value.
    integer,dimension(:),intent(out)    :: index    !! An INTEGER working array of size N.  On exit the contents
                                                    !! of this array define the sets P, Z, and F as follows:
                                                    !!
                                                    !! * INDEX(1)   through INDEX(NSETP)    =  Set P.
                                                    !! * INDEX(IZ1) through INDEX(IZ2)      = Set Z.
                                                    !! * INDEX(IZ2+1) through INDEX(N)      = Set F.
                                                    !!
                                                    !! IZ1 = NSETP + 1 = NPP1
                                                    !!
                                                    !! Any of these sets may be empty.  Set F is those components
                                                    !! that are constrained to a unique value by the given
                                                    !! constraints.   Sets P and Z are those that are allowed a non-
                                                    !! zero range of values.  Of these, set Z are those whose final
                                                    !! value is a constraint value, while set P are those whose
                                                    !! final value is not a constraint.  The value of IZ2 is not returned.
                                                    !!
                                                    !! It is computable as the number of bounds constraining a component
                                                    !! of X uniquely.
    integer,intent(out)                 :: ierr !! Indicates status on return:
                                                !!
                                                !! * 0 -- Solution completed.
                                                !! * 1 -- M  <=  0 or N  <=  0
                                                !! * 2 -- B(:), X(:), BND(:,:), W(:), or INDEX(:) size or shape violation.
                                                !! * 3 -- Input bounds are inconsistent.
                                                !! * 4 -- Exceed maximum number of iterations.
    integer,intent(in) :: max_iter !! maximum number of iterations (if <=0, then `3*n` is used)

    logical :: find, hitbnd, free1, free2, free
    integer :: iter     !! Iteration counter.
    integer :: itmax    !! Maximum number of iterations permitted.
                        !! Defaults to 3*n if `max_iter<=0`.
                        !! This is usually larger than required.
    integer :: m, n, i, ibound, ii, ip, iz, iz1, iz2, &
               j, jj, jz, l, lbound, npp1
    real(wp),dimension(size(a,2)) :: s
    real(wp),dimension(size(a,1)) :: z
    real(wp) :: alpha, asave, cc, range, &
                norm, sm, ss, t, unorm, up, ztest

    real(wp),parameter :: eps = epsilon(one)  !! Determines the relative linear dependence of a column vector
                                              !! for a variable moved from its initial value.  This is used in
                                              !! one place with the default value EPS=EPSILON(ONE).  Other
                                              !! values, larger or smaller may be needed for some problems.
                                              !! Library software will likely make this an optional argument.

    call initialize()
    ! The above call will set IERR.

    loopa: do

        ! Quit on error flag, or if all coefficients are already in the
        ! solution, .or. if M columns of A have been triangularized.
        if (ierr /=  0  .or.  iz1  > iz2 .or. nsetp >= m) exit loopa

        call select_another_coeff_to_solve_for()

        ! See if no index was found to be moved from set Z to set P.
        ! Then go to termination.
        if ( .not. find ) exit loopa

        call move_j_from_set_z_to_set_p()

        call test_set_p_against_constraints()

        ! The above call may set IERR.
        ! All coefficients in set P are strictly feasible.  Loop back.

    end do loopa

    call termination()

contains  ! These are internal subroutines.

subroutine initialize()

   m=size(a,1)
   n=size(a,2)
   if (m  <=  0 .or. n  <=  0) then
      ierr = 1
      return
   end if

   ! Check array sizes for consistency and with M and N.
   if (size(x) < n) then
      ierr=2
      return
   end if
   if (size(b) < m) then
      ierr=2
      return
   end if
   if (size(bnd,1) /= 2) then
      ierr=2
      return
   end if
   if (size(bnd,2) < n) then
      ierr=2
      return
   end if
   if (size(w) < n) then
      ierr=2
      return
   end if
   if (size(index) < n) then
      ierr=2
      return
   end if

   ierr = 0
   if (max_iter<=0) then
      itmax = 3*n
   else
      itmax = max_iter
   end if
   iter=0
   ! Initialize the array index().
   do i=1,n
      index(i)=i
   end do

   iz2=n
   iz1=1
   nsetp=0
   npp1=1

   ! Begin:  Loop on IZ to initialize  X().
   iz=iz1
   do
      if (iz  >  iz2 ) exit
      j=index(iz)
      if ( bnd(1,j)   <= -huge(one)) then
         if (bnd(2,j) >=  huge(one)) then
            x(j) = zero
         else
            x(j) = min(zero,bnd(2,j))
         end if
     else  if ( bnd(2,j) >= huge(one)) then
        x(j) = max(zero,bnd(1,j))
     else
        range = bnd(2,j) - bnd(1,j)
        if ( range <= zero ) then
           ! Here X(J) is constrained to a single value.
           index(iz)=index(iz2)
           index(iz2)=j
           iz=iz-1
           iz2=iz2-1
           x(j)=bnd(1,j)
           w(j)=zero
         else  if ( range  >  zero) then
            !! The following statement sets X(J) to 0 if the constraint interval
            !! includes 0, and otherwise sets X(J) to the endpoint of the
            !! constraint interval that is closest to 0.
            x(j) = max(bnd(1,j), min(bnd(2,j),zero))
         else
            ierr = 3
            return
         end if
   end if
      if ( abs(x(j)) > zero ) then
         ! Change B() to reflect a nonzero starting value for X(J).
         b(1:m)=b(1:m)-a(1:m,j)*x(j)
      end if
      iz=iz+1
   end do

end subroutine initialize

subroutine select_another_coeff_to_solve_for()

   !! 1. Search through set z for a new coefficient to solve for.
   !!    First select a candidate that is either an unconstrained
   !!    coefficient or else a constrained coefficient that has room
   !!    to move in the direction consistent with the sign of its dual
   !!    vector component.  Components of the dual (negative gradient)
   !!    vector will be computed as needed.
   !! 2. For each candidate start the transformation to bring this
   !!    candidate into the triangle, and then do two tests:  Test size
   !!    of new diagonal value to avoid extreme ill-conditioning, and
   !!    the value of this new coefficient to be sure it moved in the
   !!    expected direction.
   !! 3. If some coefficient passes all these conditions, set FIND = true,
   !!    The index of the selected coefficient is J = INDEX(IZ).
   !! 4. If no coefficient is selected, set FIND = false.

   find = .false.
   do iz=iz1,iz2
      j=index(iz)

      ! Set FREE1 = true if X(J) is not at the left end-point of its
      ! constraint region.
      ! Set FREE2 = true if X(J) is not at the right end-point of its
      ! constraint region.
      ! Set FREE = true if X(J) is not at either end-point of its
      ! constraint region.

     free1 = x(j) > bnd(1,j)
     free2 = x(j) < bnd(2,j)
     free = free1 .and. free2

     if ( free ) then
        call test_coef_j_for_diag_elt_and_direction_of_change()
     else
        ! Compute dual coefficient W(J).
        w(j)=dot_product(a(npp1:m,j),b(npp1:m))

        ! Can X(J) move in the direction indicated by the sign of W(J)?
        if ( w(j) < zero ) then
           if ( free1 ) call test_coef_j_for_diag_elt_and_direction_of_change()
        else  if ( w(j) > zero ) then
           if ( free2 ) call test_coef_j_for_diag_elt_and_direction_of_change()
        end if
     end if
     if ( find ) return
   end do

end subroutine select_another_coeff_to_solve_for

subroutine test_coef_j_for_diag_elt_and_direction_of_change()

   !! 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 htc (npp1, a(1:m,j), up)
   unorm = nrm2(a(1:nsetp,j))
   if ( abs(a(npp1,j)) > eps * unorm) then

      ! Column J is sufficiently independent.  Copy b into Z, update Z.
      z(1:m)=b(1:m)
      ! Compute product of transormation and updated right-hand side.
      norm=a(npp1,j)
      a(npp1,j)=up
      if (abs(norm) > zero) then
         sm=dot_product(a(npp1:m,j)/norm, z(npp1:m))/up
         z(npp1:m)=z(npp1:m)+sm*a(npp1:m,j)
         a(npp1,j)=norm
      end if

      if (abs(x(j)) >  zero) z(1:npp1)=z(1:npp1)+a(1:npp1,j)*x(j)
      ! Adjust Z() as though X(J) had been reset to zero.
      if ( free ) then
         find = .true.
      else
         !! Solve for ZTEST ( proposed new value for X(J) ).
         !! Then set FIND to indicate if ZTEST has moved away from X(J) in
         !! the expected direction indicated by the sign of W(J).
         ztest=z(npp1)/a(npp1,j)
         find = ( w(j) < zero  .and.  ztest < x(j) ) .or. &
                ( w(j) > zero  .and.  ztest > x(j) )
      end if
   end if

   ! If J was not accepted to be moved from set Z to set P,
   ! restore A(NNP1,J).  Failing these tests may mean the computed
   ! sign of W(J) is suspect, so here we set W(J) = 0.  This will
   ! not affect subsequent computation, but cleans up the W() array.
   if ( .not. find ) then
      a(npp1,j)=asave
      w(j)=zero
   end if

end subroutine test_coef_j_for_diag_elt_and_direction_of_change

subroutine move_j_from_set_z_to_set_p()

   !! The index  J=index(IZ)  has been selected to be moved from
   !! set Z to set P.  Z() contains the old B() adjusted as though X(J) = 0.
   !! A(*,J) contains the new Householder transformation vector.

   b(1:m)=z(1:m)
   index(iz)=index(iz1)
   index(iz1)=j
   iz1=iz1+1
   nsetp=npp1
   npp1=npp1+1
   ! The following loop can be null or not required.
   norm=a(nsetp,j)
   a(nsetp,j)=up
   if (abs(norm) > zero) then
      do jz=iz1,iz2
         jj=index(jz)
         sm=dot_product(a(nsetp:m,j)/norm, a(nsetp:m,jj))/up
         a(nsetp:m,jj)=a(nsetp:m,jj)+sm*a(nsetp:m,j)
      end do
   a(nsetp,j)=norm
   end if
   ! The following loop can be null.
   do l=npp1,m
      a(l,j)=zero
   end do!  L

   w(j)=zero

   ! Solve the triangular system.  Store this solution temporarily in Z().
   do i = nsetp, 1, -1
      if (i /= nsetp) z(1:i)=z(1:i)-a(1:i,ii)*z(i+1)
      ii=index(i)
      z(i)=z(i)/a(i,ii)
   end do

end subroutine move_j_from_set_z_to_set_p

subroutine test_set_p_against_constraints()

   loopb: do
      ! The solution obtained by solving the current set P is in the array Z().

      iter=iter+1
      if (iter > itmax) then
         ierr = 4
         exit loopb
      end if

      call see_if_all_constrained_coeffs_are_feasible()

      ! The above call sets HITBND.  If HITBND = true then it also sets
      ! ALPHA, JJ, and IBOUND.
      if ( .not. hitbnd ) exit loopb

      ! Here ALPHA will be between 0 and 1 for interpolation
      ! between the old X() and the new Z().
      do ip=1,nsetp
         l=index(ip)
         x(l)=x(l)+alpha*(z(ip)-x(l))
      end do

      i=index(jj)
      ! Note:  The exit test is done at the end of the loop, so the loop
      ! will always be executed at least once.
      do

         ! Modify A(*,*), B(*) and the index arrays to move coefficient I
         ! from set P to set Z.

         call move_coef_i_from_set_p_to_set_z

         if (nsetp <= 0) exit loopb

         ! See if the remaining coefficients 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 infeasible or on a boundary
         ! will be set to the boundary value and moved from set P to set Z.

          ibound = 0
          do jj=1,nsetp
             i=index(jj)
             if ( x(i) <= bnd(1,i)) then
                 ibound=1
                 exit
             else if ( x(i) >= bnd(2,i)) then
                 ibound=2
                 exit
             end if
         end do
         if (ibound <= 0) exit
      end do

      ! Copy B( ) into Z( ).  Then solve again and loop back.
      z(1:m)=b(1:m)

      do i = nsetp, 1, -1
         if (i /= nsetp) z(1:i)=z(1:i)-a(1:i,ii)*z(i+1)
         ii=index(i)
         z(i)=z(i)/a(i,ii)
      end do
   end do loopb

   ! The following loop can be null.
   do ip=1,nsetp
      i=index(ip)
      x(i)=z(ip)
   end do

end subroutine test_set_p_against_constraints

subroutine see_if_all_constrained_coeffs_are_feasible()

   !! See if each coefficient in set P is strictly interior to its constraint region.
   !! If so, set HITBND = false.
   !! If not, set HITBND = true, and also set ALPHA, JJ, and IBOUND.
   !! Then ALPHA will satisfy  0.  < ALPHA  <=  1.

   alpha=two
   do ip=1,nsetp
      l=index(ip)
      if (z(ip) <= bnd(1,l)) then
         ! Z(IP) HITS LOWER BOUND
         lbound=1
      else  if (z(ip) >= bnd(2,l)) then
         ! Z(IP) HITS UPPER BOUND
         lbound=2
      else
         lbound = 0
      end if

      if ( lbound /= 0 ) then
         t=(bnd(lbound,l)-x(l))/(z(ip)-x(l))
         if (alpha > t) then
           alpha=t
           jj=ip
           ibound=lbound
         end if
      end if
   end do
   hitbnd = abs(alpha - two) > zero

end subroutine see_if_all_constrained_coeffs_are_feasible

subroutine move_coef_i_from_set_p_to_set_z()

   x(i)=bnd(ibound,i)
   if (abs(x(i)) > zero .and. jj > 0) b(1:jj)=b(1:jj)-a(1:jj,i)*x(i)

   ! The following loop can be null.
   do j = jj+1, nsetp
      ii=index(j)
      index(j-1)=ii
      call rotg (a(j-1,ii),a(j,ii),cc,ss)
      sm=a(j-1,ii)
      ! The plane rotation is applied to two rows of A and the right-hand
      ! side.  One row is moved to the scratch array S and then the updates
      ! are computed.  The intent is for array operations to be performed
      ! and minimal extra data movement.  One extra rotation is applied
      ! to column II in this approach.
      s=a(j-1,1:n)
      a(j-1,1:n)=cc*s+ss*a(j,1:n)
      a(j,1:n)=cc*a(j,1:n)-ss*s
      a(j-1,ii)=sm
      a(j,ii)=zero
      sm=b(j-1)
      b(j-1)=cc*sm+ss*b(j)
      b(j)=cc*b(j)-ss*sm
   end do

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

end subroutine move_coef_i_from_set_p_to_set_z

subroutine termination()

   if (ierr <= 0) then
      ! Compute the norm of the residual vector.
      sm=zero
      if (npp1 <= m) then
         sm=nrm2(b(npp1:m))
      else
         w(1:n)=zero
      end if
      rnorm=sm
   end if

end subroutine termination

pure subroutine rotg(sa,sb,c,s)

   real(wp),intent(inout) :: sa
   real(wp),intent(in)    :: sb
   real(wp),intent(out)   :: c
   real(wp),intent(out)   :: s

   real(wp) :: roe,scale,r

   roe = sb
   if ( abs(sa) > abs(sb) ) roe = sa
   scale = abs(sa) + abs(sb)
   if ( scale <= zero ) then
      c = one
      s = zero
   else
      r = scale*sqrt((sa/scale)**2 + (sb/scale)**2)
      if (roe < zero) r=-r
      c = sa/r
      s = sb/r
      sa = r
   end if

end subroutine rotg

pure function nrm2 (x) result(norm)

   !!  NRM2 returns the Euclidean norm of a vector via the function
   !!  name, so that
   !!
   !!  `NRM2 := sqrt( x'*x )`
   !!
   !!### See also
   !!  * [[dnrm2]]

   real(wp),dimension(:),intent(in) :: x
   real(wp) :: norm

   real(wp) :: absxi, scale, ssq
   integer :: n, ix

   n=size(x)
   if ( n < 1) then
      norm  = zero
   else if ( n == 1 ) then
      norm  = abs( x( 1 ) )
   else
      scale = zero
      ssq   = one
      do ix = 1, n
         absxi = abs( x( ix ) )
          if (absxi > zero ) then
             if ( scale < absxi ) then
                ssq   = one + ssq*( scale/absxi )**2
                scale = absxi
             else
                ssq   = ssq + ( absxi/scale )**2
             end if
         end if
     end do
     norm  = scale * sqrt( ssq )
   end if

end function nrm2

pure subroutine htc (p, u, up)

   !! Construct a Householder transformation.

   integer,intent(in)  :: p
   real(wp),dimension(:),intent(inout) :: u
   real(wp),intent(out) :: up

   real(wp) :: vnorm

   vnorm=nrm2(u(p:size(u)))
   if (u(p) > zero) vnorm=-vnorm
   up=u(p)-vnorm
   u(p)=vnorm

end subroutine htc

end subroutine bvls