SUBROUTINE lfnorm(n, m, x, y, beta, z, ky, ifault)

!     ALGORITHM AS 135  APPL. STATIST. (1979) VOL.28, NO. 1

!     Min-Max estimates for a linear multiple regression problem

! Input arguments
! ---------------
! N      Number of cases
! M      Number of predictors
! X(,)   Values of the predictor variables, one case per row.
! Y()    Values of the dependent variable

! Output arguments
! ----------------
! BETA() Final estimates of the regression coefficients
! Z      The least maximum absolute deviation
! KY     The number of iterations
! IFAULT Error indicator = 0 for normal termination
!                        = 1 if the observation matrix has less than full rank
!
! N.B. Arguments NDIM & MDIM have been removed.

! Fortran 90 version by Alan.Miller @ vic.cmis.csiro.au
! Latest revision - 5 January 1999

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN)     :: m
REAL (dp), INTENT(IN)   :: x(:, :)   ! x(ndim, mdim)
REAL (dp), INTENT(IN)   :: y(:)
REAL (dp), INTENT(OUT)  :: beta(:)
REAL (dp), INTENT(OUT)  :: z
INTEGER, INTENT(OUT)    :: ky
INTEGER, INTENT(OUT)    :: ifault

REAL (dp) :: hilo(m), xrxf(m), xsxf(m), lu(m, m), dev1, sigr, sumxr,  &
             deviat, yest, ydev, sumxs, ratio, test, delta, div, swing,    &
             SAVE, sigs, top, acu, big
INTEGER   :: i, ibase(m), ii, ii1, irow, k, kk, kkk, k1, k2, m1, sss, rrr,  &
             indx(m)
LOGICAL   :: intl

REAL (dp), PARAMETER :: zero = 0._dp, one = 1._dp, two = 2._dp

acu = SQRT( EPSILON(1.0_dp) )
big = SQRT( HUGE(1.0_dp) )
ifault = 0
ky = 0
z = zero
m1 = m - 1

!     Set up initial LU decomposition

DO i = 1, m
  indx(i) = i
END DO
intl = .true.
kkk = 1
CALL update(kkk, x, lu, ibase, indx, intl, n, m, ifault)
IF (ifault /= 0) RETURN
intl = .false.
irow = kkk

!     Calculate beta value

k = indx(1)
k1 = ibase(1)
beta(k) = y(k1) / lu(k, 1)
DO ii = 2, m
  k = indx(ii)
  k1 = ibase(ii)
  beta(k) = y(k1)
  ii1 = ii - 1
  DO i = 1, ii1
    kk = indx(i)
    beta(k) = beta(k) - lu(kk, ii) * beta(kk)
  END DO
  beta(k) = beta(k) / lu(k, ii)
END DO
DO ii = 1, m1
  k1 = m - ii
  k = indx(k1)
  DO i = 1, ii
    kk = m - i + 1
    k2 = indx(kk)
    beta(k) = beta(k) - lu(k2, k1) * beta(k2)
  END DO
END DO

!     Search for and set first violated constraint as R-th constraint

50 irow = irow + 1
IF (irow > n) RETURN
dev1 = zero
DO i = 1, m
  dev1 = dev1 + x(irow, i) * beta(i)
END DO
dev1 = dev1 - y(irow)
IF (ABS(dev1) < acu) GO TO 50
sigr = SIGN(one, dev1)
rrr = irow

!     Adjust for the R-th constraint

k = indx(1)
xrxf(1) = x(rrr, k)
DO ii = 2, m
  k = indx(ii)
  xrxf(ii) = x(rrr, k)
  ii1 = ii - 1
  DO i = 1, ii1
    xrxf(ii) = xrxf(ii) - lu(k, i) * xrxf(i)
  END DO
END DO
k = indx(m)
xrxf(m) = xrxf(m) / lu(k, m)
hilo(m) = SIGN(one, -sigr * xrxf(m))
sumxr = sigr - hilo(m) * xrxf(m)
DO ii = 1, m1
  k1 = m - ii
  k = indx(k1)
  DO i = 1, ii
    k2 = m - i + 1
    xrxf(k1) = xrxf(k1) - lu(k, k2) * xrxf(k2)
  END DO
  xrxf(k1) = xrxf(k1) / lu(k, k1)
  hilo(k1) = SIGN(one, -sigr * xrxf(k1))
  sumxr = sumxr - hilo(k1) * xrxf(k1)
END DO
z = ABS(dev1 / sumxr)

!     Start of main iterative loop.
!     Search for the most violated S-th constraint

110 sss = 0
deviat = acu

!     Calculate beta value

k = indx(1)
k1 = ibase(1)
beta(k) = (y(k1) + z * hilo(1)) / lu(k, 1)
DO ii = 2, m
  k = indx(ii)
  k1 = ibase(ii)
  beta(k) = y(k1) + z * hilo(ii)
  ii1 = ii - 1
  DO i = 1, ii1
    kk = indx(i)
    beta(k) = beta(k) - lu(kk, ii) * beta(kk)
  END DO
  beta(k) = beta(k) / lu(k, ii)
END DO
DO ii = 1, m1
  k1 = m - ii
  k = indx(k1)
  DO i = 1, ii
    kk = m - i + 1
    k2 = indx(kk)
    beta(k) = beta(k) - lu(k2, k1) * beta(k2)
  END DO
END DO

!     Calculate residuals

DO i = 1, n
  yest = DOT_PRODUCT( x(i, 1:m), beta(1:m) )
  dev1 = ABS(y(i) - yest) - z
  IF (dev1 <= deviat) CYCLE
  ydev = yest - y(i)
  deviat = dev1
  sss = i
END DO

!     Check if at optimum

IF (sss == 0) RETURN

!     Set up information on the S-th constraint

sigs = SIGN(one, ydev)
k = indx(1)
xsxf(1) = x(sss, k)
DO ii = 2, m
  k = indx(ii)
  xsxf(ii) = x(sss, k)
  ii1 = ii - 1
  DO i = 1, ii1
    xsxf(ii) = xsxf(ii) - lu(k, i) * xsxf(i)
  END DO
END DO
k = indx(m)
xsxf(m) = xsxf(m) / lu(k, m)
sumxs = -sigs + hilo(m) * xsxf(m)
DO ii = 1, m1
  k1 = m - ii
  k = indx(k1)
  DO i = 1, ii
    k2 = m - i + 1
    xsxf(k1) = xsxf(k1) - lu(k, k2) * xsxf(k2)
  END DO
  xsxf(k1) = xsxf(k1) / lu(k, k1)
  sumxs = sumxs + hilo(k1) * xsxf(k1)
END DO

!     Search for minimum ratio

210 kkk = 0
ratio = big
DO i = 1, m
  IF (sigs * SIGN(one, xsxf(i)) /= hilo(i) .OR. ABS(xsxf(i)) < acu) CYCLE
  test = ABS(xrxf(i) / xsxf(i))
  IF (test >= ratio) CYCLE
  ratio = test
  kkk = i
END DO

!     Check if R-th constraint moves interior

IF (kkk /= 0) GO TO 260

!     Process the movement of the R-th constraint

delta = ABS(deviat / sumxs)

!     Calculate the largest tolerable delta

div = ABS(sumxr) - two
IF (div < acu) GO TO 240
swing = two * z / div
IF (swing >= delta) GO TO 240

!     Switch R and S constraint indicators

SAVE = sumxs
sumxs = -sumxr + sigr + sigr
sumxr = -SAVE
SAVE = sigr
sigr = sigs
sigs = -SAVE
deviat = ABS(sumxs * delta) - two * z
z = z + delta
DO i = 1, m
  SAVE = xsxf(i)
  xsxf(i) = xrxf(i)
  xrxf(i) = SAVE
END DO
i = rrr
rrr = sss
sss = i
GO TO 210

!     Replace the R-th constraint with the S-th constraint

240 sigr = sigs
xrxf(1:m) = xsxf(1:m)
sumxr = -sumxs
z = z + delta
rrr = sss
GO TO 110

!     Process the movement of the K-th constraint

260 delta = ABS(xrxf(kkk) * deviat / (xrxf(kkk) * sumxs + xsxf(kkk) * sumxr))
top = -two * z * xrxf(kkk)
div = xrxf(kkk) * xrxf(kkk) + hilo(kkk) * sumxr
IF (SIGN(one, top) /= SIGN(one, div)) GO TO 270
IF (ABS(div) < acu) GO TO 270
swing = top / div

!     Check to see if the K-th constraint swings across

IF (swing >= delta) GO TO 270
z = z + swing
deviat = deviat - swing * ABS(sumxs + xsxf(kkk) * sumxr / xrxf(kkk))
sumxr = sumxr + two * hilo(kkk) * xrxf(kkk)
sumxs = sumxs - two * hilo(kkk) * xsxf(kkk)
hilo(kkk) = -hilo(kkk)
GO TO 210

!     Update XRXF and the LU of the current basis

270 hilo(kkk) = sigs
sumxr = sigr
xrxf(kkk) = xrxf(kkk) / xsxf(kkk)
sumxr = sumxr - hilo(kkk) * xrxf(kkk)
DO i = 1, m
  IF (i == kkk) CYCLE
  xrxf(i) = xrxf(i) - xsxf(i) * xrxf(kkk)
  sumxr = sumxr - hilo(i) * xrxf(i)
END DO
ibase(kkk) = sss

!     Update LU decomposition

CALL update(kkk, x, lu, ibase, indx, intl, n, m, ifault)
IF (ifault /= 0) RETURN
z = z + delta
ky = ky + 1
GO TO 110

CONTAINS


SUBROUTINE update(kkk, x, lu, ibase, indx, intl, n, m, ifault)

!     ALGORITHM AS 135.1  APPL. STATIST. (1979) VOL.28, NO. 1

!     Update LU decomposition matrix

INTEGER, INTENT(IN OUT)   :: kkk
REAL (dp), INTENT(IN)     :: x(:, :)   ! x(ndim, mdim)
REAL (dp), INTENT(IN OUT) :: lu(:, :)
INTEGER, INTENT(IN OUT)   :: ibase(:)
INTEGER, INTENT(IN OUT)   :: indx(:)
LOGICAL, INTENT(IN OUT)   :: intl
INTEGER, INTENT(IN)       :: n
INTEGER, INTENT(IN)       :: m
INTEGER, INTENT(OUT)      :: ifault

! Local variables

REAL (dp) :: subt, pivot
INTEGER   :: i, icol, ii, ii1, irow, isave, j, k, kk

irow = 0
DO ii = kkk, m
  IF (intl) GO TO 10
  irow = ibase(ii)
  GO TO 20

  10 irow = irow + 1
  IF (irow > n) THEN
    ifault = 1
    RETURN
  END IF
  20 lu(1:m, ii) = x(irow, 1:m)
  
!     Set up representation of incoming row
  
  IF (ii == 1) GO TO 60
  ii1 = ii - 1
  DO icol = 1, ii1
    k = indx(icol)
    subt = lu(k, ii)
    j = icol + 1
    DO i = j, m
      k = indx(i)
      lu(k, ii) = lu(k, ii) - subt * lu(k, icol)
    END DO
  END DO
  
!     Find maximum entry
  
  60 pivot = acu
  kk = 0
  DO i = ii, m
    k = indx(i)
    IF (ABS(lu(k, ii)) <= pivot) CYCLE
    pivot = ABS(lu(k, ii))
    kk = i
  END DO
  IF (kk == 0) GO TO 10
  
!     Switch order
  
  isave = indx(kk)
  indx(kk) = indx(ii)
  indx(ii) = isave
  
!     Put into columns of LU one at a time
  
  IF (intl) ibase(ii) = irow
  IF (ii == m) CYCLE
  j = ii + 1
  DO i = j, m
    k = indx(i)
    lu(k, ii) = lu(k, ii) / lu(isave, ii)
  END DO
END DO
kkk = irow

RETURN
END SUBROUTINE update

END SUBROUTINE lfnorm



PROGRAM test_lfnorm
IMPLICIT NONE

INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER                :: i, ier, niter, m, n
REAL (dp), ALLOCATABLE :: x(:,:), y(:), beta(:)
REAL (dp)              :: z

INTERFACE
  SUBROUTINE lfnorm(n, m, x, y, beta, z, ky, ifault)
    IMPLICIT NONE
    INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
    INTEGER, INTENT(IN)     :: n
    INTEGER, INTENT(IN)     :: m
    REAL (dp), INTENT(IN)   :: x(:, :)   ! x(ndim, mdim)
    REAL (dp), INTENT(IN)   :: y(:)
    REAL (dp), INTENT(OUT)  :: beta(:)
    REAL (dp), INTENT(OUT)  :: z
    INTEGER, INTENT(OUT)    :: ky
    INTEGER, INTENT(OUT)    :: ifault
  END SUBROUTINE lfnorm
END INTERFACE

DO
  WRITE(*, '(a)', ADVANCE='NO') ' Enter number of predictor variables: '
  READ(*, *) m
  IF (m < 1) THEN
    WRITE(*, *) '** Must be greater than zero **'
    CYCLE
  ELSE
    EXIT
  END IF
END DO

n = 3 * m

ALLOCATE( x(n,m), y(n), beta(m) )

! Generate artificial data with coefficients 1, 2, ..., m.

beta = (/ (DBLE(i),i=1,m) /)
DO i = 1, n
  CALL RANDOM_NUMBER( x(i,2:) )
  CALL RANDOM_NUMBER( z )
  x(i,1) = 1.0_dp
  y(i) = DOT_PRODUCT( x(i,:), beta ) + z - 0.5_dp
END DO

CALL lfnorm(n, m, x, y, beta, z, niter, ier)
IF (ier > 0) THEN
  WRITE(*, *) '** X-matrix has rank < m **'
ELSE
  WRITE(*, *) 'Regression coefficients'
  WRITE(*, '(t2, 10f7.3)') beta
  WRITE(*, '(a, i6)') ' No. of iterations = ', niter
  WRITE(*, '(a, f9.5)') ' Least max. abs. deviation = ', z
END IF

STOP
END PROGRAM test_lfnorm