        subroutine erckcmp(iopt,iropt,zval,nref1,nref2)

        implicit double precision(a-h, o-z)
        parameter (kmax=50, nbmax=5, nqmax=5, nmax=90000)
        parameter (mxitr=100, tol=1.0d-5,didl=1.0d-5)
        parameter (hllow=1.0d-8)
        external dlinrg, dlinds, dmxtyf, dmrrrr, dmurrv, dmxytf

        dimension x(kmax,nbmax), a(kmax,kmax), f(nqmax,kmax,kmax)
        dimension biggam(nqmax,nqmax), kv(nmax), del(nmax), t(nmax)
        dimension b(kmax,kmax), pih(kmax), cnt(nmax), iwt(nmax)
        dimension xtil(kmax,nbmax), xtilex(kmax,nbmax)
        dimension ilow(nmax), iup(nmax)
        dimension sst(kmax), dcl(nmax,kmax)
        dimension rsk(kmax), ev(kmax), grn(kmax)
        dimension sstar(nmax,kmax), green(nmax,kmax), gstar(nmax,kmax)
        dimension s(nmax,kmax), svar(nmax,kmax), hl(kmax)
        dimension q(nmax,kmax,kmax)
        dimension omega(kmax,kmax), phi(kmax,nqmax)
        dimension ominv(kmax,kmax), xttomi(nbmax,kmax)
        dimension desmat(nbmax,nbmax), v(nbmax,nbmax), rhs(nbmax)
        dimension beta(nbmax), diff(nbmax), bdiff(nbmax), evec(nmax)
        dimension dlhl(kmax), omnew(kmax,kmax), xtinv(nbmax,nbmax)
        dimension xtinvom(nbmax,nbmax), omni(kmax,kmax)
        dimension xttomni(nbmax,kmax)

        common /data0/ n, k, nb, nq
        common /data1/ kv, del, t
        common /data2/ cnt, ilow, iup
        common /data3/ x, a, f, biggam
        common /data4/ beta, v, ifail

        dn = n
        nb1 = nb + 1
        ipth = 1
        ifail=0
        nref1 = 0
        nref2 = 0

C       COMPUTE B (INVERSE OF A)
        call dlinrg(k,a,kmax,b,kmax)

C       SET UP EXTENDED DESIGN MATRIX
        do l = 1, k
          xtil(l,1) = 1.0d0
          do ir = 1, nb
            xtil(l,ir+1) = x(l,ir)
           enddo
         enddo

C       TALLY COUNTS IN EACH CONFIG
        do l = 1, k
          pih(l) = 0.0d0
         enddo
        do i = 1, n
          pih(kv(i)) = pih(kv(i))+1.0d0
         enddo

C       SORT THE DATA AND IDENTIFY UNIQUE F-U TIMES
        call sort(n,ndist)

C       INITIALIZATIONS
        do l = 1, k
          sst(l) = 1.0d0
          grn(l) = 0.0d0
          hl(l) = 0.0d0
         enddo

C       LOOP OVER F-U TIMES

        do m = 1, ndist
         
C         COMPUTATION OF KAPLAN-MEIER AND CUML HAZ ESTIMATES

          do l = 1, k
            rsk(l) = 0.0d0
            ev(l) = 0.0d0
           enddo

          do i = ilow(m), n
            rsk(kv(i)) = rsk(kv(i)) + 1.0d0
           enddo

          do i = ilow(m), iup(m)
            ev(kv(i)) = ev(kv(i)) + del(i)
           enddo

          do l = 1, k
            p = ev(l)/rsk(l)
            sst(l) = sst(l)*(1.0d0-p)
            grn(l) = grn(l) + p/(rsk(l)*(1.0d0-p))
            sstar(m,l) = sst(l)
            green(m,l) = grn(l)
            dcl(m,l) = p
            gstar(m,l) = rsk(l)/pih(l)
           enddo
 
C         END OF COMPUTATION OF K-M AND CUML HAZ ESTIMATES

C         COMPUTE S = B*SSTAR AND ITS VARIANCE
          do l1 = 1, k
            s(m,l1) = 0.0d0
            svar(m,l1) = 0.0d0
            do l2 = 1, k
              addtrm = b(l1,l2)*sstar(m,l2)
              s(m,l1) = s(m,l1) + addtrm
              svar(m,l1) = svar(m,l1) + (addtrm**2)*green(m,l2)
             enddo
           enddo

         enddo

C       END OF LOOP OVER F-U TIMES

C       IDENTIFY FIRST AND LAST TIMES AT WHICH BOUNDS ON S ARE MET
        istflg=0
        do m = 1, ndist
          iwflg = 1
          do l = 1, k
            chwid = zval*dsqrt(svar(m,l))
            shi = s(m,l) + chwid
            slo = s(m,l) - chwid
            if ((iropt .eq. 1) .and. (shi .gt. 1.0d0)) then
              iwflg=0
              nref1 = nref1 + 1
             endif
            if (slo .le. 0.0d0) then
              iwflg=0
              nref2 = nref2 + 1
             endif
           enddo
          iwt(m) = iwflg
          if ((istflg .eq. 0) .and. (iwflg .eq. 1)) then
            mlow=m
            istflg=1
           endif
          if (iwflg .eq. 1) mhi=m
         enddo
        if (istflg .eq. 0) then
          write(10,*) 'Computation Failed - S Bounds Violated'
          ifail=1
          return
         endif

C       COMPUTE L'S
        totcnt = 0.0d0
        do m = mlow, mhi
          if (iwt(m) .eq. 1) then
            totcnt = totcnt + cnt(m)
            do l = 1, k
              hlkm = -dlog(s(m,l))
              hl(l) = hl(l) + cnt(m)*hlkm
             enddo
           endif
         enddo
        do l = 1, k
          hl(l) = hl(l)/totcnt
         enddo 

C       COMPUTE Q

        do ir = 1, k
        do is = 1, k

          do m1 = 1, mhi
            q(m1,ir,is) = 0.0d0
            mstart = max(m1,mlow)
            do m2 = mstart, mhi
              if (iwt(m2) .eq. 1) then
                addtrm = cnt(m2)*sstar(m2,is)/s(m2,ir)
                q(m1,ir,is) = q(m1,ir,is) + addtrm
               endif
             enddo
            q(m1,ir,is) = q(m1,ir,is)/totcnt
           enddo

         enddo
         enddo

C       COMPUTE OMEGA ASSUMING KNOWN CLASSIFICATION RATES

        do ir = 1, k
        do is = 1, k

          omega(ir,is) = 0.0d0

          do m = 1, mhi
            do l = 1, k
              addtrm = q(m,ir,l)*q(m,is,l)/gstar(m,l)
              addtrm = addtrm*b(ir,l)*b(is,l)*dcl(m,l)/(pih(l)/dn)
              omega(ir,is) = omega(ir,is) + addtrm/dn
             enddo
           enddo

          write(3,*) ir,is,omega(ir,is)            
   
         enddo
         enddo

C       ADJUST OMEGA FOR ESTIMATION OF CLASSIFICATION RATES

        if (nq .gt. 0) then

          do l = 1, k
          do iu = 1, nq
            phi(l,iu) = 0.0d0
           enddo
           enddo

          do m = mlow, mhi
            if (iwt(m) .eq. 1) then
              do l = 1, k
              do iu = 1, nq
              do ir = 1, k
              do is = 1, k
                addtrm = (1.0d0/s(m,l))*b(l,ir)
     $            *f(iu,ir,is)*s(m,is)*cnt(m)
                phi(l,iu) = phi(l,iu) + addtrm/totcnt
               enddo
               enddo
               enddo
               enddo
             endif
           enddo

          write(3,*) 'PHI'
          do l = 1, k
          do iu = 1, nq
            write(3,*) l, iu, pih(l), phi(l,iu)
           enddo
           enddo   
   

          do ir = 1, k
          do is = 1, k
          do it = 1, nq
          do iu = 1, nq
            addtrm = phi(ir,it)*phi(is,iu)*biggam(it,iu)
            omega(ir,is) = omega(ir,is) + addtrm
           enddo
           enddo

          write(3,*) ir,is,omega(ir,is)            

           enddo
           enddo

         endif

C       COMPUTE WEIGHTED LEAST SQUARES ESTIMATES

C       SETUP UP LOG-TRANSFORMED RESPONSE VECTOR 
C       AND CORRESPONDING COVARIANCE MATRIX
        do l = 1, k
          if (hl(l) .gt. hllow)
     $      then
              dlhl(l) = dlog(hl(l))
            else
              ifail=1
              write(10,*) 'Computation Failed - L Value Too Near 1'
              return
            endif
         enddo
        do ir = 1, k
          omnew(ir,ir) = omega(ir,ir)/(hl(ir)**2)
          do is = ir+1, k
            omnew(ir,is) = omega(ir,is)/(hl(ir)*hl(is))
            omnew(is,ir) = omnew(ir,is)
           enddo
          enddo

C       SQUARE CASE 
        if (nb1 .eq. k) then
          call dlinrg(nb1,xtil,kmax,xtinv,nbmax)
          call dmurrv(nb1,nb1,xtinv,nbmax,nb1,dlhl,ipth,k,beta)
          call dmrrrr(nb1,nb1,xtinv,nbmax,nb1,nb1,omnew,kmax,nb1,nb1,
     $      xtinvom,nbmax)
          call dmxytf(nb1,nb1,xtinvom,nbmax,nb1,nb1,xtinv,nbmax,
     $      nb1,nb1,v,nbmax)
          goto 100
         endif

C       NON-ITERATIVE ESTIMATE
        call dlinds(k,omnew,kmax,omni,kmax)
        call dmxtyf(k,nb1,xtil,kmax,k,k,omni,kmax,nb1,k,
     $    xttomni,nbmax)
        call dmrrrr(nb1,k,xttomni,nbmax,k,nb1,xtil,kmax,
     $    nb1,nb1,desmat,nbmax)
        call dlinrg(nb1,desmat,nbmax,v,nbmax)
        call dmurrv(nb1,k,xttomni,nbmax,k,dlhl,ipth,nb1,rhs)
        call dmurrv(nb1,nb1,v,nbmax,nb1,rhs,ipth,nb1,beta)

C       IF REQUESTED, PROCEED TO ITERATIVE ESTIMATE. OTHERWISE, FINISH.
        if (iopt .ne. 1) goto 100         
        call dlinds(k,omega,kmax,ominv,kmax)

C       WE USE HERE A NAIVE IMPLEMENTATION OF GAUSS-NEWTON FITTING
C       IF THERE ARE PROBLEMS WITH CONVERGENCE, A MORE SOPHISTICATED
C       ALGORITHM MAY BE NEEDED

        itr = 0
        icflg = 0

        do while ((itr .lt. mxitr) .and. (icflg .eq. 0))

          itr = itr + 1

          do l = 1, k
            evec(l) = 0.0d0
            do ib = 1, nb1
              evec(l) = evec(l) + xtil(l,ib)*beta(ib)
             enddo
            evec(l) = dexp(evec(l))
           enddo

           do l = 1, k
             diff(l) = hl(l) - evec(l)
           do ib = 1, nb1
             xtilex(l,ib) = evec(l)*xtil(l,ib)
            enddo
            enddo
 
          call dmxtyf(k,nb1,xtilex,kmax,k,k,ominv,kmax,nb1,k,
     $      xttomi,nbmax)
          call dmrrrr(nb1,k,xttomi,nbmax,k,nb1,xtilex,kmax,
     $      nb1,nb1,desmat,nbmax)
          do ib = 1, nb1
            desmat(ib,ib) = desmat(ib,ib)+didl
           enddo
          call dlinds(nb1,desmat,nbmax,v,nbmax)
          call dmurrv(nb1,k,xttomi,nbmax,k,diff,ipth,nb1,rhs)
          call dmurrv(nb1,nb1,v,nbmax,nb1,rhs,ipth,nb1,bdiff)

          icflg = 1
          do ib = 1, nb1
            beta(ib) = beta(ib) + bdiff(ib)
            rdiff = dabs(bdiff(ib))/dmax1(1.0d0,dabs(beta(ib)))
            if (rdiff .gt. tol) icflg=0
           enddo

         enddo

        ifail = 1-icflg

100     continue 

        if (ifail .ne. 0) 
     $    then
            write(10,*) 'Computation Failed - G-N Did Not Converge'
c          else
c            do ib1 = 1, nb1
c              v(ib1,ib1) = v(ib1,ib1)/dn
c              do ib2 = ib1+1, nb1
c                v(ib1,ib2) = v(ib1,ib2)/dn
c                v(ib2,ib1) = v(ib1,ib2)
c               enddo
c             enddo
          endif

        return
        end

C**********************************************************************

        subroutine sort(nobs,ndis)

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

        dimension t(nmax), del(nmax), kv(nmax), indx(nmax), wksp(nmax)
        dimension ilow(nmax), iup(nmax), cnt(nmax)

        common /data1/ kv, del, t
        common /data2/ cnt, ilow, iup

C       SORT DATA
        do i = 1, nobs
          indx(i) = i
         enddo
        call dsvrgp(nobs,t,t,indx)
        do i = 1, nobs
          wksp(i) = del(i)
         enddo
        do i = 1, nobs
          del(i) = wksp(indx(i))
         enddo
        do i = 1, nobs
          wksp(i) = kv(i)
         enddo
        do i = 1, nobs
          kv(i) = wksp(indx(i))
         enddo

C       IDENTIFY UNIQUE F-U TIMES
        ix = 1
        told = t(1)
        ilow(1) = 1
        do i  = 2, nobs
          if (t(i) .ne. told) then
            iup(ix) = i-1
            ix = ix+1
            ilow(ix) = i
            told = t(i)
           endif
          if (i .eq. nobs)  iup(ix) = nobs
         enddo
        ndis = ix
        do m = 1, ndis
          ncnt = iup(m) - ilow(m) + 1
          cnt(m) = ncnt
         enddo

        return
        end
