SUBROUTINE nonnegls(nin, b)
! Non-negative least squares fitting after INCLUD has been used to form
! an orthogonal reduction of the data.

! Author: Alan.Miller @ vic.cmis.csiro.au
! Latest version - 19 October 1998

USE lsq
IMPLICIT NONE

INTEGER, INTENT(OUT)         :: nin
REAL (lsq_kind), INTENT(OUT) :: b(:)

! Local variables

REAL (lsq_kind)            :: alpha, sxy(ncol), sxx(ncol), temp, wmax, z(ncol)
INTEGER                    :: i, ier, t
REAL (lsq_kind), PARAMETER :: zero = 0.0_lsq_kind

! The step numbers below are those for Lawson & Hanson's NNLS algorithm
! from page 161 of their book:
! Lawson, C.L. & Hanson, R.J. `Solving least squares problems', Prentice-Hall,
! 1974, republished in SIAM Classics in Applied Mathematics 1995.

! Step 1

nin = 0
b(:ncol) = zero

! Step 2

2 CALL sums_of_sq_and_prod(nin+1, ncol, sxy, sxx)

! Step 4
! We find the max of (sxy)^2 / sxx conditional upon sxy > 0.
! This is better than L & H's (sxy - b'.sxx) which is dependent upon the size
! of the X-variables and so tends to pick X-variables which have large values
! ahead of X-variables with small values.

t = 0
wmax = zero
DO i = nin+1, ncol
  IF (sxy(i) > zero) THEN
    temp = sxy(i)**2 / sxx(i)
    IF (temp > wmax) THEN
      wmax = temp
      t = i
    END IF
  END IF
END DO

! Step 3

IF (t == 0) RETURN

! Step 5, add the variable from position t to the model.

nin = nin + 1
IF (t > nin) CALL vmove(t, nin, ier)

! Step 6, calculate the regression coefficients.

6 CALL regcf(z, nin, ier)

! Step 7, if all coefficients > 0, set b = z, then go back to step 2, unless
! all the variables have been included.

IF (ALL(z(:nin) > zero)) THEN
  b(:nin) = z(:nin)
  IF (nin == ncol) RETURN
  GO TO 2
END IF

! Steps 8-11, one or more negative regression coefficients.
! Steps 8-9, calculate shrinkage factor alpha.

t = 0
alpha = 1.0
nin = nin - 1
DO i = 1, nin
  IF (z(i) <= zero) THEN
    temp = b(i) / (b(i) - z(i))
    IF (temp < alpha) THEN
      alpha = temp
      t = i
    END IF
  END IF
END DO

! Step 10, shrink regression coefficients

IF (t > 0) THEN
  b(1:nin) = b(1:nin) + alpha * (z(1:nin) - b(1:nin))
ELSE
  GO TO 2                    ! This should never happen
END IF

! Step 11, remove variables with coefficients <= 0.

CALL vmove(t, nin, ier)
nin = nin - 1
GO TO 6


CONTAINS


SUBROUTINE sums_of_sq_and_prod(first, last, sxy, sxx)
! Calculate the sums of squares, and cross-products, of those parts of the
! X and Y variables which are orthogonal to the X-variables already in the
! model (if any) in positions 1 .. first-1.  Ignore any variables after last.

INTEGER, INTENT(IN)          :: first, last
REAL (lsq_kind), INTENT(OUT) :: sxy(:), sxx(:)

!     Local variables

INTEGER         :: inc, pos, row, col
REAL (lsq_kind) :: zero = 0.0, diag, dy

!     Accumulate sums of squares & products from row FIRST.
!     The orthogonal reduction gave:
!     X = Q.sqrt(D).R   where Q is orthonormal and D is diagonal and such
!     that the diagonal elements of R are 1's.
!     The sums of squares and products are given by  R'DR
!     where only rows & columns FIRST .. LAST of R & D are used.

sxx(first:last) = zero
sxy(first:last) = zero
inc = ncol - last
pos = row_ptr(first)
DO row = first, last
  diag = d(row)
  dy = diag * rhs(row)
  sxx(row) = sxx(row) + diag
  sxy(row) = sxy(row) + dy
  DO col = row+1, last
    sxx(col) = sxx(col) + diag * r(pos)**2
    sxy(col) = sxy(col) + dy * r(pos)
    pos = pos + 1
  END DO
  pos = pos + inc
END DO

RETURN
END SUBROUTINE sums_of_sq_and_prod

END SUBROUTINE nonnegls