SUBROUTINE gllm(ni, nid, nj, nk, nkp, ji, y, c, conv, e, f, cspr, cslr, ifault)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-04-26  Time: 11:05:17

! N.B. Arguments W & V have been removed - they were workspaces.

!    Algorithm AS 207 Appl. Statist. (1984) vol.33, no.3, pp. 358-362.
!    by Michael Haber

!    Fitting a generalized log-linear model
!    to fully or partially classified frequencies.

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

INTEGER, INTENT(IN)        :: ni, nid, nj, nk, nkp, ji(ni)
REAL (dp), INTENT(IN)      :: y(nj)
REAL (dp), INTENT(IN OUT)  :: c(nid,nkp)
REAL (dp), INTENT(IN)      :: conv
REAL (dp), INTENT(OUT)     :: e(ni), f(nj), cspr, cslr
INTEGER, INTENT(OUT)       :: ifault

REAL (dp)  :: w(4,ni), v(2,nkp)
REAL (dp)  :: cmin, csum, ctmax, ff, yy
INTEGER    :: i, j, k, nkk
LOGICAL    :: lconv
REAL (dp), PARAMETER  :: eps = 0.00001_dp, zero = 0.0_dp, one = 1.0_dp

ifault = 1

!       Check the ji array

DO  i = 1, ni
  IF (ji(i) < 1 .OR. ji(i) > nj) RETURN
END DO
ifault = 0

!       initialize

e(1:ni) = one
w(3,1:nj) = zero

!       standardize the C matrix

cmin = zero
DO  i = 1, ni
  DO  k = 1, nk
    cmin = MIN(cmin,c(i,k))
  END DO
END DO
IF (cmin /= zero) THEN
    c(1:ni,1:nk) = c(1:ni,1:nk) - cmin
END IF

80 ctmax = zero
DO  i = 1, ni
  csum = zero
  DO  k = 1, nk
    csum = csum + c(i,k)
  END DO
  ctmax = MAX(ctmax,csum)
  w(4,i) = csum
END DO

IF (ctmax <= eps) THEN
  ifault = 2
  RETURN
END IF
IF (ABS(ctmax-one) > eps) THEN
  DO  i = 1, ni
    DO  k = 1, nk
      c(i,k) = c(i,k) / ctmax
    END DO
  END DO
  GO TO 80
END IF
DO  i = 1, ni
  IF (ABS(w(4,i)-one) > eps) GO TO 140
END DO
nkk = nk
GO TO 160
140 nkk = nkp
DO  i = 1, ni
  c(i,nkk) = one - w(4,i)
END DO

!       Enter the EM algorithm

160 f(1:nj) = zero
DO  i = 1, ni
  j = ji(i)
  f(j) = f(j) + e(i)
END DO

!       Check for convergence

lconv = .true.
DO  j = 1, nj
  IF (ABS(f(j)-w(3,j)) > conv) lconv = .false.
  w(3,j) = f(j)
END DO
IF (.NOT.lconv) THEN
  
  DO  i = 1, ni
    j = ji(i)
    w(1,i) = y(j)
    IF (f(j) > eps) w(1,i) = e(i) * y(j) / f(j)
  END DO
  DO  k = 1, nkk
    v(1,k) = SUM( c(1:ni,k) * w(1,1:ni) )
  END DO
  
!       Enter the IPF algorithm
  
  230 w(2,1:ni) = e(1:ni)
  DO  k = 1, nkk
    v(2,k) = SUM( c(1:ni,k) * e(1:ni) )
    DO  i = 1, ni
      IF (c(i,k) > eps .AND. v(2,k) > eps) e(i) = e(i) * (v(1,k)/  &
          v(2,k)) ** c(i,k)
    END DO
  END DO
  DO  i = 1, ni
    IF (ABS(e(i)-w(2,i)) > conv) GO TO 230
  END DO
  GO TO 160
END IF

!       calculate the goodness-of-fit statistics

cspr = zero
cslr = zero
DO  j = 1, nj
  yy = y(j)
  ff = f(j)
  IF (ff > eps) THEN
    cspr = cspr + (yy-ff) ** 2 / ff
    IF (yy > eps) THEN
      cslr = cslr + yy * LOG(yy/ff)
    END IF
  END IF
END DO
cslr = 2.0D0 * cslr
RETURN
END SUBROUTINE gllm