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

!*****************************************************************************80
!
!! BVLS solves a least squares problem with bounds on the variables.
!
!  Discussion:
!
!    Given an M by N matrix, A, and an M-vector, B,  compute an
!    N-vector, X, that solves the least-squares problem 
!      A * X = B
!    subject to the lower and upper bounds:
!      BND(1,1:N) <= X(J) <= BND(2,1:N)
!
!    In cases where no bound is intended, simply set the corresponding entry
!    of BND to a very negative or very positive value.
!
!    This algorithm is a generalization of NNLS, which solves
!    the least-squares problem,
!      A * X = B
!    subject to the constraints
!      0 <= X(1:N).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    22 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns in A.
!
!    Input/output, real ( kind = 8 ) A(M,N).
!    On input, A() contains the M by N matrix, A.
!    On output, 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).
!
!    Input/output, real ( kind = 8 ) B(M).
!    On input, B() contains the M-vector, B.
!    On output, B() contains Q*B.  The same Q multiplies A.
!
!    Input, real ( kind = 8 ) BND(2,N), lower and upper bounds for X.
!    BND(1,J) <= X(J) <= BND(2,J).
!
!    Output, real ( kind = 8 ) X(N), the solution vector.
!
!    Output, real ( kind = 8 ) RNORM, the Euclidean norm of the residual 
!    vector, b - A*X.
!
!    Output, integer ( kind = 4 ) NSETP, the number of components of the 
!    solution vector, X(), that are not at their constraint values.
!
!    Output, real ( kind = 8 ) W(N), 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.
!
!    Output, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F:
!    * 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.
!
!    Output, integer ( kind = 4 ) IERR, status flag.
!    * 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.
!
!  Local Parameters:
!
!   ITMAX  [integer]
!       Set to 3*N.  Maximum number of iterations permitted.
!       This is usually larger than required.  Library software will
!       likely make this an optional argument.
!
!   ITER   [integer]
!       Iteration counter.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  logical ( kind = 4 ) find
  logical ( kind = 4 ) free
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) iter
  integer ( kind = 4 ) itmax
  integer ( kind = 4 ) iz
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) iz2
  integer ( kind = 4 ) j
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ), parameter :: one = 1.0D+00
  real ( kind = 8 ) rnorm
  real ( kind = 8 ) up
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) z(m)

  call initialize ( m, n, a, b, bnd, ierr, index, iter, itmax, &
    iz, iz1, iz2, j, npp1, nsetp, w, x )

  do
!
!  Quit (failure) on error flag.
!
    if ( ierr /= 0 ) then
      exit
    end if
!
!  Quit (success) if all coefficients are already in the solution
!
    if ( iz2 < iz1 ) then
      exit
    end if
!
!  Quit (success) if M columns of A have been triangularized.
!
    if ( m <= nsetp ) then
      exit
    end if

    call select_another_coeff ( m, n, iz, iz1, iz2, j, npp1, nsetp, free, &
      up, b, bnd, index, x, a, w, find, z )
!
!  Quit if no index was found to move from set Z to set P.
!
    if ( .not. find ) then
      exit
    end if

    call move_index_j ( m, n, iz, iz2, j, up, iz1, npp1, a, index, w, &
      z, nsetp, b )

    call test_set_p_against_constraints ( m, n, a, b, bnd, index, iter, &
      itmax, iz1, npp1, nsetp, x, z, ierr )

  end do

  call termination ( ierr, m, n, npp1, b, w, rnorm )

  return
end
subroutine bvls_report ( m, n, a, b, bnd, x, rnorm, nsetp, w, index, ierr )

!*****************************************************************************80
!
!! BVLS_REPORT reports on the results of a successful call to BVLS.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    25 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns in A.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input, real ( kind = 8 ) B(M), the right-hand-side vector, B.
!
!    Input, real ( kind = 8 ) BND(2,N), lower and upper bounds on X.
!
!    Input, real ( kind = 8 ) X(N), the solution vector.
!
!    Input, real ( kind = 8 ) RNORM, the Euclidean norm of the residual 
!    vector, b - A*X.
!
!    Input, integer ( kind = 4 ) NSETP, the number of components of the 
!    solution vector, X(), that are not at their constraint values.
!
!    Input, real ( kind = 8 ) W(N), the dual solution vector. 
!
!    Input, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F:
!    * 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.
!
!    Input, integer ( kind = 4 ) IERR, status flag.
!    * 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.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) r(m)
  real ( kind = 8 ) rnorm
  real ( kind = 8 ) rnorm2
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) w2(n)
  real ( kind = 8 ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'BVLS_REPORT:'  

  if ( ierr /= 0 ) then
    write ( *, '(a,i8)' ) '  Abnormal error flag, IERR = ', ierr
    return
  end if

  write ( *, '(a,i8)' ) '  Number of components not at constraints =', nsetp
!
!  X
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Solution vector, X:'
  write ( *, '(a)' ) ''
  write ( *, '(2x,5g14.6)' ) x(1:n)
!
!  INDEX
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Variable index INDEX:'
  write ( *, '(a)' ) ''
  write ( *, '(2x,5i14)' ) index(1:n)
!
!  R = B - A * X
!
  r(1:m) = b(1:m) - matmul ( a(1:m,1:n), x(1:n) )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Residual R = B - A*X:'
  write ( *, '(a)' ) ' '
  write ( *, '(2x,5g14.6)' ) r(1:m)
!
!  ||R||
!
  rnorm2 = sqrt ( dot_product ( r(1:m), r(1:m) ) )
  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.5)') '  Residual norm           = ', rnorm2
  write ( *, '(a,g17.5)') '  Residual norm from BVLS = ', rnorm
!
!  W2 = A'*R
!
  w2(1:n) = matmul ( r(1:m), a(1:m,1:n) )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Dual vector: W = (A'')*R:'
  write ( *, '(a)' ) ' '
  write ( *, '(2x,5g14.6)' ) w2(1:n)
!
!  W: Dual vector.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Dual vector from BVLS: W'
  write ( *, '(a)' ) ' '
  write ( *, '(2x,5g14.6)' ) w(1:n)

  return
end
subroutine constrained_feasible ( nsetp, index, bnd, x, z, hitbnd, alpha, &
  jj, ibound )

!*****************************************************************************80
!
!! CONSTRAINED_FEASIBLE determines if each coefficient is interior.
!
!  Discussion:
!
!    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.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) NSETP.
!
!    Input, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Input, real ( kind = 8 ) BND(2,N).
!
!    Input, real ( kind = 8 ) X(N).
!
!    Input, real ( kind = 8 ) Z(M).
!
!    Output, logical HITBND.
!
!    Output, real ( kind = 8 ) ALPHA.
!
!    Output, integer ( kind = 4 ) JJ.
!
!    Output, integer ( kind = 4 ) IBOUND.
!
  implicit none

  integer ( kind = 4 ) nsetp

  real ( kind = 8 ) alpha
  real ( kind = 8 ) bnd(2,*)
  logical ( kind = 4 ) hitbnd
  integer ( kind = 4 ) ibound
  integer ( kind = 4 ) index(*)
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lbound
  real ( kind = 8 ) t
  real ( kind = 8 ) x(*)
  real ( kind = 8 ) z(*)

  alpha = 2.0D+00

  do ip = 1, nsetp

    l = index(ip)
!
!  Z(IP) hits lower bound
!
    if ( z(ip) <= bnd(1,l) ) then

      lbound = 1
!
!  Z(IP) hits upper bound
!
    else if ( bnd(2,l) <= z(ip) ) then

      lbound = 2
!
!  Z(IP) is in the interior:
!
    else
      lbound = 0
    end if

    if ( lbound /= 0 ) then
      t = ( bnd(lbound,l) - x(l) ) / ( z(ip) - x(l) )
      if ( t < alpha ) then
        alpha = t
        jj = ip
        ibound = lbound
      end if
    end if

  end do

  hitbnd = ( 0.0D+00 < abs ( alpha - 2.0D+00 ) )

  return
end
subroutine htc ( p, m, u, up )

!*****************************************************************************80
!
!! HTC constructs a Householder transformation.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) P, an index into the U vector.
!
!    Input, integer ( kind = 4 ) M, the length of U.
!
!    Input/output, real ( kind = 8 ) U(M).
!
!    Output, real ( kind = 8 ) UP.
!
  implicit none

  integer ( kind = 4 ) m

  real ( kind = 8 ) nrm2
  integer ( kind = 4 ) p
  real ( kind = 8 ) u(m)
  real ( kind = 8 ) up
  real ( kind = 8 ) vnorm

  vnorm = nrm2 ( m + 1 - p, u(p:m) )

  if ( 0.0D+00 < u(p) ) then
    vnorm = - vnorm
  end if

  up = u(p) - vnorm
  u(p) = vnorm

  return
end
subroutine initialize ( m, n, a, b, bnd, ierr, index, iter, itmax, &
  iz, iz1, iz2, j, npp1, nsetp, w, x )

!*****************************************************************************80
!
!! INITIALIZE initializes internal data.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    22 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the dimensions of the matrix.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input, real ( kind = 8 ) B(M).
!
!    Input, real ( kind = 8 ) BND(2,N).
!
!    Output, integer ( kind = 4 ) IERR.
!
!    Output, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Output, integer ( kind = 4 ) ITER.
!
!    Output, integer ( kind = 4 ) ITMAX.
!
!    Output, integer ( kind = 4 ) IZ.
!
!    Output, integer ( kind = 4 ) IZ1.
!
!    Output, integer ( kind = 4 ) IZ2.
!
!    Output, integer ( kind = 4 ) J.
!
!    Output, integer ( kind = 4 ) NPP1.
!
!    Output, integer ( kind = 4 ) NSETP.
!
!    Output, real ( kind = 8 ) W(N).
!
!    Output, real ( kind = 8 ) X(N).
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) iter
  integer ( kind = 4 ) itmax
  integer ( kind = 4 ) iz
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) iz2
  integer ( kind = 4 ) j
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ), parameter :: one = 1.0D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) x(n)

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

  iz2 = n
  iz1 = 1
  nsetp = 0
  npp1 = 1
!
!  Loop on IZ to initialize X().
!
  iz = iz1

  do

    if ( iz2 < iz ) then
      exit
    end if

    j = index(iz)

    if ( bnd(1,j) <= - huge ( one ) ) then

      if ( huge ( one ) <= bnd(2,j) ) then
        x(j) = 0.0D+00
      else
        x(j) = min ( 0.0D+00, bnd(2,j) )
      end if

    else if ( huge ( one ) <= bnd(2,j) ) then

      x(j) = max ( 0.0D+00, bnd(1,j) )

    else

      r = bnd(2,j) - bnd(1,j)
!
!  Here X(J) is constrained to a single value.
!
      if ( r <= 0.0D+00 ) then

        index(iz) = index(iz2)
        index(iz2) = j
        iz = iz - 1
        iz2 = iz2 - 1
        x(j) = bnd(1,j)
        w(j) = 0.0D+00
!
!  Sets X(J) to 0 if the constraint interval includes 0;
!  otherwise set X(J) to the endpoint of the constraint interval closest to 0.
!
      else if ( 0.0D+00 < r ) then

        x(j) = max ( bnd(1,j), min ( bnd(2,j), 0.0D+00 ) )

      else

        ierr = 3
        return

      end if

    end if
!
!  Change B() to reflect a nonzero starting value for X(J).
!
    if ( 0.0D+00 < abs ( x(j) ) ) then
      b(1:m) = b(1:m) - a(1:m,j) * x(j)
    end if

    iz = iz + 1

  end do

  return
end
subroutine move_coef_i ( m, n, i, ibound, bnd, jj, b, a, x, nsetp, index, &
  npp1, iz1 )

!*****************************************************************************80
!
!! MOVE_COEF_I moves coefficient I from set P to set Z.
!
!  Discussion:
!
!    Modify A(*,*), B(*) and the index arrays to move coefficient I
!    from set P to set Z.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of rows of A.
!
!    Input, integer ( kind = 4 ) N, the number of columns of A.
!
!    Input, integer ( kind = 4 ) I, the index of the coefficient to be moved.
!
!    Input, integer ( kind = 4 ) IBOUND.
!
!    Input, real ( kind = 8 ) BND(2,N).
!
!    Input, integer ( kind = 4 ) JJ.
!
!    Input/output, real ( kind = 8 ) B(M).
!
!    Input/output, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input/output, real ( kind = 8 ) X(N).
!
!    Input/output, integer ( kind = 4 ) NSETP.
!
!    Input/output, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Input/output, integer ( kind = 4 ) NPP1.
!
!    Input/output, integer ( kind = 4 ) IZ1.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  real ( kind = 8 ) cc
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ibound
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) s(n)
  real ( kind = 8 ) sm
  real ( kind = 8 ) ss
  real ( kind = 8 ) x(n)

  x(i) = bnd(ibound,i)

  if ( 0.0D+00 < abs ( x(i) ) .and. 0 < jj ) then
    b(1:jj) = b(1:jj) - a(1:jj,i) * x(i)
  end if

  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(1:n) = a(j-1,1:n)
    a(j-1,1:n) = cc * s(1:n) + ss * a(j,1:n)
    a(j,1:n) = cc * a(j,1:n) - ss * s(1:n)

    a(j-1,ii) = sm
    a(j,ii) = 0.0D+00

    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

  return
end
subroutine move_index_j ( m, n, iz, iz2, j, up, iz1, npp1, a, index, w, &
  z, nsetp, b )

!*****************************************************************************80
!
!! MOVE_INDEX_J moves the index J from set Z to set P.
!
!  Discussion:
!
!    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.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M.
!
!    Input, integer ( kind = 4 ) N.
!
!    Input, integer ( kind = 4 ) IZ.
!
!    Input, integer ( kind = 4 ) IZ2.
!
!    Input, integer ( kind = 4 ) J, the index to be moved.
!
!    Input, real ( kind = 8 ) UP.
!
!    Input/output, integer ( kind = 4 ) IZ1.
!
!    Input/output, integer ( kind = 4 ) NPP1.
!
!    Input/output, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input/output, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Input/output, real ( kind = 8 ) W(N).
!
!    Input/output, real ( kind = 8 ) Z(M).
!
!    Output, integer ( kind = 4 ) NSETP.
!
!    Output, real ( kind = 8 ) B(M).
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) iz
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) iz2
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) jz
  integer ( kind = 4 ) l
  real ( kind = 8 ) norm
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) sm
  real ( kind = 8 ) up
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) z(m)

  b(1:m) = z(1:m)

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

  norm = a(nsetp,j)
  a(nsetp,j) = up

  if ( 0.0D+00 < abs ( norm ) ) 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

  a(npp1:m,j) = 0.0D+00

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

  return
end
function nrm2 ( n, x )

!*****************************************************************************80
!
!! NRM2 returns the Euclidean norm of a vector.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vector.
!
!    Input, real ( kind = 8 ) X(N), the vector.
!
!    Output, real ( kind = 8 ) NRM2, the norm of the vector.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) absxi
  integer ( kind = 4 ) ix
  real ( kind = 8 ) norm
  real ( kind = 8 ) nrm2
  real ( kind = 8 ) scale
  real ( kind = 8 ) ssq
  real ( kind = 8 ) x(n)

  if ( n < 1 ) then
    norm = 0.0D+00
  else if ( n == 1 ) then
    norm = abs ( x(1) )
  else
    scale = 0.0D+00
    ssq = 1.0D+00

    do ix = 1, n
      absxi = abs ( x(ix) )
      if ( 0.0D+00 < absxi ) then
        if ( scale < absxi ) then
          ssq = 1.0D+00 + ssq * ( scale / absxi )**2
          scale = absxi
        else
          ssq = ssq + ( absxi / scale )**2
        end if
      end if
    end do
    norm  = scale * sqrt ( ssq )
  end if

  nrm2 = norm

  return
end
subroutine r8mat_uniform_01 ( m, n, seed, r )

!*****************************************************************************80
!
!! R8MAT_UNIFORM_01 fills an R8MAT with unit pseudorandom numbers.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns in
!    the array.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(M,N), the array of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_huge = 2147483647
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real ( kind = 8 ) r(m,n)

  do j = 1, n

    do i = 1, m

      k = seed / 127773

      seed = 16807 * ( seed - k * 127773 ) - k * 2836

      if ( seed < 0 ) then
        seed = seed + i4_huge
      end if

      r(i,j) = real ( seed, kind = 8 ) * 4.656612875D-10

    end do
  end do

  return
end
subroutine r8vec_uniform_01 ( n, seed, r )

!*****************************************************************************80
!
!! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    05 July 2006
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(N), the vector of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real ( kind = 8 ) r(n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop 1
  end if

  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + 2147483647
    end if

    r(i) = real ( seed, kind = 8 ) * 4.656612875D-10

  end do

  return
end
subroutine rotg ( sa, sb, c, s )

!*****************************************************************************80
!
!! ROTG defines a Givens rotation.
!
!  Discussion:
!
!    Given values A and B, this routine computes
!
!    SIGMA = sign ( A ) if abs ( A ) >  abs ( B )
!          = sign ( B ) if abs ( A ) <= abs ( B );
!
!    R     = SIGMA * ( A * A + B * B );
!
!    C = A / R if R is not 0
!      = 1     if R is 0;
!
!    S = B / R if R is not 0,
!        0     if R is 0.
!
!    The computed numbers then satisfy the equation
!
!    (  C  S ) ( A ) = ( R )
!    ( -S  C ) ( B ) = ( 0 )
!
!    The routine also computes
!
!    Z = S     if abs ( A ) > abs ( B ),
!      = 1 / C if abs ( A ) <= abs ( B ) and C is not 0,
!      = 1     if C is 0.
!
!    The single value Z encodes C and S, and hence the rotation:
!
!    If Z = 1, set C = 0 and S = 1;
!    If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z;
!    if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C );
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    15 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) SA, the value of A.
!    On output, SA is overwritten with R.
!
!    Input, real ( kind = 8 ) SB, the value of B.
!
!    Output, real ( kind = 8 ) C, S, the cosine and sine of the
!    Givens rotation.
!
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) r
  real ( kind = 8 ) roe
  real ( kind = 8 ) s
  real ( kind = 8 ) sa
  real ( kind = 8 ) sb
  real ( kind = 8 ) scale

  roe = sb

  if ( abs ( sb ) < abs ( sa ) ) then
    roe = sa
  end if

  scale = abs ( sa ) + abs ( sb )

  if ( scale <= 0.0D+00 ) then
    c = 1.0D+00
    s = 0.0D+00
    return
  end if

  r = scale * sqrt ( ( sa / scale ) ** 2 + ( sb / scale ) ** 2 )

  if ( roe < 0.0D+00 ) then
    r = - r
  end if

  c = sa / r
  s = sb / r
  sa = r

  return
end
subroutine select_another_coeff ( m, n, iz, iz1, iz2, j, npp1, nsetp, free, &
  up, b, bnd, index, x, a, w, find, z )

!*****************************************************************************80
!
!! SELECT_ANOTHER_COEFF searches set Z for a new coefficient to solve for.
!
!  Discussion:
!
!    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.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    18 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of rows of A.
!
!    Input, integer ( kind = 4 ) N, the number of columns of A.
!
!    Output, integer ( kind = 4 ) IZ.
!
!    Input, integer ( kind = 4 ) IZ1, IZ2.
!
!    Input, integer ( kind = 4 ) J.
!
!    Input, integer ( kind = 4 ) NPP1.
!
!    Input, integer ( kind = 4 ) NSETP.
!
!    Input, logical ( kind = 4 ) FREE.
!
!    Input, real ( kind = 8 ) UP.
!
!    Input, real ( kind = 8 ) B(M).
!
!    Input, real ( kind = 8 ) BND(2,N).
!
!    Input, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Input, real ( kind = 8 ) X(N).
!
!    Input/output, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input/output, real ( kind = 8 ) W(N).
!
!    Output, logical ( kind = 4 ) FIND.
!
!    Output, real ( kind = 8 ) Z(M).
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n 

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  logical ( kind = 4 ) find
  logical ( kind = 4 ) free
  logical ( kind = 4 ) free1
  logical ( kind = 4 ) free2
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) iz
  integer ( kind = 4 ) iz_index
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) iz2
  integer ( kind = 4 ) j
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) up
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) z(m)

  find = .false.

  do iz_index = iz1, iz2

    iz = iz_index

    j = index(iz)
!
!  FREE1 = true if X(J) is not at the left end-point of its constraint region.
!  FREE2 = true if X(J) is not at the right end-point of its constraint region.
!  FREE = true if X(J) is not at either end-point of its constraint region.
!
    free1 = ( bnd(1,j) < x(j) )
    free2 = ( x(j) < bnd(2,j) )
    free = ( free1 .and. free2 )

    if ( free ) then

      call test_coef_j ( m, n, j, npp1, nsetp, free, up, b, x, a, w, &
        find, z )
!
!  Compute dual coefficient W(J).
!
    else

      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) < 0.0D+00 ) then
        if ( free1 ) then
          call test_coef_j ( m, n, j, npp1, nsetp, free, up, b, x, a, &
            w, find, z )
        end if
      else if ( 0.0D+00 < w(j) ) then
        if ( free2 ) then
          call test_coef_j ( m, n, j, npp1, nsetp, free, up, b, x, a, &
            w, find, z )
        end if
      end if

    end if

    if ( find ) then
      return
    end if

  end do

  return
end
subroutine termination ( ierr, m, n, npp1, b, w, rnorm )

!*****************************************************************************80
!
!! TERMINATION handles the termination of the BVLS algorithm.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) IERR, error flag.
!
!    Input, integer ( kind = 4 ) M.
!
!    Input, integer ( kind = 4 ) N.
!
!    Input, integer ( kind = 4 ) NPP1.
!
!    Input, real ( kind = 8 ) B(M).
!
!    Output, real ( kind = 8 ) W(N).
!
!    Output, real ( kind = 8 ) RNORM.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) b(m)
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) npp1
  real ( kind = 8 ) nrm2
  real ( kind = 8 ) rnorm
  real ( kind = 8 ) w(n)
!
!  Compute the norm of the residual vector.
!
  rnorm = 0.0D+00

  if ( ierr <= 0 ) then

    if ( npp1 <= m ) then
      rnorm = nrm2 ( m + 1 - npp1, b(npp1:m) )
    else
      w(1:n) = 0.0D+00
    end if

  end if

  return
end
subroutine test_coef_j ( m, n, j, npp1, nsetp, free, up, b, x, a, w, &
  find, z )

!*****************************************************************************80
!
!! TEST_COEF_J
!
!  Discussion:
!
!    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.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 June 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of rows of A.
!
!    Input, integer ( kind = 4 ) N, the number of columns of A.
!
!    Input, integer ( kind = 4 ) J.
!
!    Input, integer ( kind = 4 ) NPP1.
!
!    Input, integer ( kind = 4 ) NSETP.
!
!    Input, logical ( kind = 4 ) FREE.
!
!    Input, real ( kind = 8 ) UP.
!
!    Input, real ( kind = 8 ) B(M).
!
!    Input, real ( kind = 8 ) X(N).
!
!    Input/output, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input/output, real ( kind = 8 ) W(N).
!
!    Output, logical ( kind = 4 ) FIND.
!
!    Output, real ( kind = 8 ) Z(M).
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n 

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) asave
  real ( kind = 8 ) b(m)
  real ( kind = 8 ), parameter :: eps = 2.220446049250313D-016
  logical ( kind = 4 ) find
  logical ( kind = 4 ) free
  integer ( kind = 4 ) j
  real ( kind = 8 ) norm
  integer ( kind = 4 ) npp1
  real ( kind = 8 ) nrm2
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) sm
  real ( kind = 8 ) unorm
  real ( kind = 8 ) up
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) z(m)
  real ( kind = 8 ) ztest

  asave = a(npp1,j)

  call htc ( npp1, m, a(1:m,j), up )

  unorm = nrm2 ( nsetp, a(1:nsetp,j) )

  if ( eps * unorm < abs ( a(npp1,j) ) ) then
!
!  Column J is sufficiently independent.  Copy B into Z, update Z.
!
    z(1:m) = b(1:m)
!
!  Compute product of transformation and updated right-hand side.
!
    norm = a(npp1,j)
    a(npp1,j) = up

    if ( 0.0D+00 < abs ( norm ) ) 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 ( 0.0D+00 < abs ( x(j) ) ) then
      z(1:npp1) = z(1:npp1) + a(1:npp1,j) * x(j)
    end if
!
!  Adjust Z as though X(J) had been reset to zero.
!
    if ( free ) then

      find = .true.

    else
!
!  Solve for ZTEST, the proposed new value for X(J).
!  FIND indicates 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) < 0.0D+00 .and. ztest < x(j) ) .or. &
        ( 0.0D+00 < w(j) .and. x(j) < ztest )

    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) = 0.0D+00
  end if

  return
end
subroutine test_set_p_against_constraints ( m, n, a, b, bnd, index, iter, &
  itmax, iz1, npp1, nsetp, x, z, ierr )

!*****************************************************************************80
!
!! TEST_SET_P_AGAINST_CONSTRAINTS
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    13 August 2014
!
!  Author:
!
!    Original FORTRAN90 version by Charles Lawson, Richard Hanson.
!    This FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Charles Lawson, Richard Hanson,
!    Solving Least Squares Problems,
!    SIAM, 1995,
!    ISBN: 0898713560,
!    LC: QA275.L38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input, real ( kind = 8 ) B(M).
!
!    Input, real ( kind = 8 ) BND(2,N).
!
!    Input/output, integer ( kind = 4 ) INDEX(N), defines the sets P, Z, and F.
!
!    Input/output, integer ( kind = 4 ) ITER.
!
!    Input, integer ( kind = 4 ) ITMAX.
!
!    Input/output, integer ( kind = 4 ) IZ1.
!
!    Input/output, integer ( kind = 4 ) NPP1.
!
!    Input/output, integer ( kind = 4 ) NSETP.
!
!    Input/output, real ( kind = 8 ) X(N).
!
!    Input/output, real ( kind = 8 ) Z(M).
!
!    Output, integer IERR.
!
!  Local Parameters:
!
!    ALPHA
!    HITBND
!    I
!    IBOUND
!    II
!    IP
!    JJ
!    L
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  real ( kind = 8 ) alpha
  real ( kind = 8 ) b(m)
  real ( kind = 8 ) bnd(2,n)
  logical ( kind = 4 ) hitbnd
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ibound
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) index(n)
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) iter
  integer ( kind = 4 ) itmax
  integer ( kind = 4 ) iz1
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) jjj
  integer ( kind = 4 ) l
  integer ( kind = 4 ) npp1
  integer ( kind = 4 ) nsetp
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) z(m)

  loopb: do
!
!  The solution obtained by solving the current set P is in the array Z().
!
    iter = iter + 1
    if ( itmax < iter ) then
      ierr = 4
      exit loopb
    end if
!
!  Set HITBND.
!  If HITBND is true, also set ALPHA, JJ, and IBOUND.
!
    call constrained_feasible ( nsetp, index, bnd, x, z, hitbnd, alpha, &
      jj, ibound )

    if ( .not. hitbnd ) then
      exit loopb
    end if
!
!  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)
!
!  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 ( m, n, i, ibound, bnd, jj, b, a, x, nsetp, index, &
        npp1, iz1 )

      if ( nsetp <= 0 ) then
        exit loopb
      end if
!
!  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
      jj = nsetp

      do jjj = 1, nsetp

        i = index(jjj)

        if ( x(i) <= bnd(1,i) ) then
          ibound = 1
          jj = jjj
          exit
        else if ( bnd(2,i) <= x(i) ) then
          ibound = 2
          jj = jjj
          exit
        end if

      end do

      if ( ibound <= 0 ) then
        exit
      end if

    end do
!
!  Copy B into Z.
!
    z(1:m) = b(1:m)
!
!  Solve and loop back again.
!
    do i = nsetp, 1, -1
      if ( i /= nsetp ) then
        z(1:i) = z(1:i) - a(1:i,ii) * z(i+1)
      end if
      ii = index(i)
      z(i) = z(i) / a(i,ii)
    end do

  end do loopb

  do ip = 1, nsetp
    i = index(ip)
    x(i) = z(ip)
  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end