subroutine sort

        implicit double precision (a-h,o-z)
        parameter (nmax=5000,kbmax=20)
        external dsvrgp

        dimension t(nmax), del(nmax), z(nmax,kbmax)
        dimension indx(nmax), wksp(nmax)
        dimension ilow(nmax), iup(nmax), rsk(nmax), evcnt(nmax)

        common /data1/  n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist

C       SORT DATA
        do i = 1, n
          indx(i) = i
         enddo
        call dsvrgp(n,t,t,indx)
        do i = 1, n
          wksp(i) = del(i)
         enddo
        do i = 1, n
          del(i) = wksp(indx(i))
         enddo
        do j = 1, kb
          do i = 1, n
            wksp(i) = z(i,j)
           enddo
          do i = 1, n
            z(i,j) = wksp(indx(i))
           enddo
         enddo

C       IDENTIFY UNIQUE F-U TIMES
        ix = 1
        told = t(1)
        ilow(1) = 1
        rsknum = n
        rsk(1) = rsknum
        rm = 1.0d0
        evcnt(1) = del(1)
        do i  = 2, n
          if (t(i) .ne. told) then
            iup(ix) = i-1
            rsk(ix) = rsknum
            rsknum = rsknum - rm
            ix = ix+1
            ilow(ix) = i
            rm = 0.0d0
            evcnt(ix) = 0.0d0
            told = t(i)
           endif
          rm = rm + 1.0d0
          evcnt(ix) = evcnt(ix) + del(i)
          if (i .eq. n)  iup(ix) = n
         enddo
        ndist = ix

        return
        end


        subroutine rcomp(bta,rho,ifail)

        implicit double precision (a-h,o-z)
        parameter(kbmax=20,nmax=5000)
        parameter(didl=1.0d-2,rskmin=20.0d0)
        parameter(slbak=10.0d0,slbmin=10.0d0)

        dimension t(nmax), del(nmax), z(nmax,kbmax), bta(kbmax)
        dimension r(nmax), p(nmax), dlamgam(nmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data4/ r

        pp = 1.0d0
        ndmx = ndist
        icut = 0
        ifail = 0
        evtt = 0.0d0

        m = 1

        do 100 while ((m .le. ndist) .and. (icut .eq. 0)) 

          if (rsk(m) .lt. rskmin) icut=1

          pfac = 1.0d0 - rho*evcnt(m)/rsk(m)

          if (pfac .le. didl) icut=1

          if (icut .eq. 1) then
            ndmx = m-1 
            goto 90
           endif

          pp = pp*pfac
          p(m) = pp

          hg = 0.0d0
          if (evcnt(m) .eq. 0.0d0) goto 85
          evtt = evtt + 1.0d0
          do 30  i = ilow(m), iup(m)
            if (del(i) .eq. 0.0d0)  goto 29
            zsum = 0.0d0
            do 10  jj = 1, kb
              zsum = zsum + z(i,jj)*bta(jj)
   10        continue
            gami = dexp(-zsum)
            hg = hg + gami
   29      continue
   30      continue

   85     continue

          dlamgam(m) = hg/rsk(m)

  90      continue

          m = m + 1

 100     continue

        if (evtt .le. slbmin) then
          ifail = 1
          return
         endif

        cumsum = dlamgam(1)
        r(1) = cumsum/p(1)

        do 200  m = 2, ndmx
          cumsum = cumsum + p(m-1)*dlamgam(m)
          r(m) = cumsum/p(m)
 200     continue

        if (ndmx .eq. ndist) return

        tmax = t(ndmx)
        rmax = r(ndmx)
        t2sum = 0.0d0
        rtsum = 0.0d0
        m = ndmx-1
        slb = dmin1(evtt,slbak)
        evtim = 0.0d0
        do 280 while (evtim .lt. slb) 
          if (evcnt(m) .gt. 0.0d0) then
            evtim = evtim + 1.0d0
            td = t(ilow(m)) - tmax
            t2sum = t2sum + td**2
            rtsum = rtsum + td*(r(m)-rmax)
           endif
          m = m-1
  280    continue
        slope = rtsum/t2sum
        if (slope .lt. 0.0d0) slope=0.0d0
        do 290  m = ndmx+1, ndist
          td = t(ilow(m)) - tmax
          r(m) = rmax + slope*td
  290    continue

        return
        end


        subroutine pslik(np,pars,func)

        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000,didl=1.d-8)

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(kbmax), r(nmax)

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data4/ r

        dn = n

        do 10  jj = 1, kb
          beta(jj) = pars(jj)
   10    continue
        rho = pars(kb+1)+didl

        fllik = 0.d0

        do 50  m = 1, ndist
        do 50  i = ilow(m), iup(m)

          zsum = 0.d0
          do 30  jj = 1, kb
            zsum = zsum + z(i,jj)*beta(jj)
   30      continue
          gam = dexp(zsum)

          extrml = dlog(1.0d0 + rho*r(m)*gam)

          if (del(i) .eq. 1.0d0) then
            fllik = fllik + zsum - extrml
           endif
          fllik = fllik - extrml/rho

   50    continue

        func = -fllik/dn

        return
        end


        subroutine psgrad(np,pars,grad)

        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000,didl=1.d-8)

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension beta(kbmax), grad(kmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)
        dimension r(nmax)

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data4/ r

        dn = n

        do 10  jj = 1, kb
          beta(jj) = pars(jj)
   10    continue
        rho = pars(kb+1)+didl

        do 20  l = 1, k
          grad(l) = 0.0d0
   20    continue

        do 70  m = 1, ndist
        do 70  i = ilow(m), iup(m)

          zsum = 0.d0
          do 30  jj = 1, kb
            zsum = zsum + z(i,jj)*beta(jj)
   30      continue
          gam = dexp(zsum)

          rre = rho*r(m)*gam
          extrm = 1.0d0 + rre
          extrml = dlog(extrm)
          exrat = rre/extrm
          exrm = r(m)*gam/extrm

          if (del(i) .eq. 1.0d0) then
            do 35  l = 1, kb
              grad(l) = grad(l) - z(i,l)/extrm
   35        continue
            grad(kb+1) = grad(kb+1) + exrm
           endif
          do 40  l = 1, kb
            grad(l) = grad(l) + z(i,l)*exrm
   40      continue
          grad(kb+1) = grad(kb+1) + (exrat-extrml)/(rho**2)

   70    continue

        do 80  l = 1, k
          grad(l) = grad(l)/dn
   80    continue

        return
        end


        subroutine rtilcmp(pars)

        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000)

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension r(nmax), rtil(nmax), gam(nmax)
        dimension rsk(nmax), ilow(nmax), iup(nmax), evcnt(nmax)        

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist        
        common /data4/ r
        common /data5/ rtil

        do 10  m = 1, ndist
          rtil(m) = r(m)
   10    continue
        return
        end

c        rho = pars(kb+1)
c
c        do 20  i = 1, n
c          zsum = 0.0d0
c          do 15 l = 1, kb
c            zsum = zsum + z(i,l)*pars(l)
c   15      continue
c          gam(i) = dexp(zsum)
c   20    continue
c
c        rtcum = 0.0d0
c
c        do 70  m = 1, ndist
c
c          rden = 0.0d0
c
c          do 40  j = ilow(m), n
c            rden = rden + gam(j)/(1.0d0+rho*r(m)*gam(j))
c   40      continue
c
c          rtcum = rtcum + evcnt(m)/rden
c          rtil(m) = rtcum
c
c   70    continue
c
c         return 
c         end


        subroutine marlik(np,pars,func)

        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000,didl=1.d-8)

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension rtil(nmax), zsu(nmax), gam(nmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)   

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data5/ rtil

        dn = n

        rho = pars(kb+1)+didl

        do 20  i = 1, n
          zsum = 0.0d0
          do 15 l = 1, kb
            zsum = zsum + z(i,l)*pars(l)
   15      continue
          zsu(i) = zsum
          gam(i) = dexp(zsum)
   20    continue

        fllik = 0.0d0

        do 55  m = 1, ndist
          if (evcnt(m) .eq. 0.0d0)  goto 54
          do 40  i = ilow(m), iup(m)
            if (del(i) .ne. 0.0d0) then
              extrml = dlog(1.0d0 + rho*rtil(m)*gam(i))
              fllik = fllik + zsu(i) - extrml
             endif
   40      continue

          cyrr = 0.0d0
          do 50  j = ilow(m), n
            cyrr = cyrr + gam(j)/(1.0d0+rho*rtil(m)*gam(j))
   50      continue

          fllik = fllik - evcnt(m)*dlog(cyrr)

   54    continue
   55    continue

        func = -fllik/dn

        return
        end


        subroutine margrd(np,pars,grad)

        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000,didl=1.d-8)

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension gam(nmax), grad(kmax)
        dimension expnum(kmax), rtil(nmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data5/ rtil

        dn = n

        rho = pars(kb+1)+didl

        do 20  i = 1, n
          zsum = 0.0d0
          do 15 l = 1, kb
            zsum = zsum + z(i,l)*pars(l)
   15      continue
          gam(i) = dexp(zsum)
   20    continue

        do 30  l = 1, k
          grad(l) = 0.0d0
   30    continue

        do 70  m = 1, ndist

          if (evcnt(m) .eq. 0.0d0) goto 69

          do 50  i = ilow(m), iup(m)
            if (del(i) .eq. 0.0d0)  goto 49
            exti = 1.0d0 + rho*rtil(m)*gam(i)
            do 35  l = 1, kb
              grad(l) = grad(l) - z(i,l)/exti
   35        continue
            grad(kb+1) = grad(kb+1) + rtil(m)*gam(i)/exti
   49      continue
   50      continue

          expden = 0.0
          do 55  l = 1, k
            expnum(l) = 0.0d0
   55      continue

          do 65  j = ilow(m), n
            extj = 1.0d0+rho*rtil(m)*gam(j)
            wj = gam(j)/extj
            expden = expden + wj
            do 62  l = 1, kb
              expnum(l) = expnum(l) + wj*z(j,l)/extj
   62        continue
            expnum(kb+1) = expnum(kb+1) + (wj**2)*rtil(m)
   65      continue

          do 68  l = 1, kb
            grad(l) = grad(l) + evcnt(m)*expnum(l)/expden
   68      continue
          grad(kb+1) = grad(kb+1) - evcnt(m)*expnum(kb+1)/expden

   69    continue
   70    continue

        do 80  l = 1, k
          grad(l) = grad(l)/dn
   80    continue

        return
        end