        program ernewcox

c       program to implement method of Zucker and Spiegelman (2008, Stat Med)
c       version of 21 April 2010

        implicit double precision(a-h, o-z)
        parameter (kmax=20, nbmax=20, nmax=90000, nqmax=5, nsimul=100)

        dimension x(kmax,nbmax), a(kmax,kmax)
        dimension adot(nqmax,kmax,kmax)
        dimension b(kmax,kmax), bdot(nqmax,kmax,kmax)
        dimension borig(kmax,kmax), bidnt(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), id(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(nbmax), cov(nbmax,nbmax)
        dimension biggam(nqmax,nqmax)
        dimension xnmat(kmax,nbmax), xnmor(kmax,nbmax)
        dimension xdmat(nqmax,kmax,nbmax)
        dimension csc(nbmax), djac(nbmax,nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist
        common /data4/ xnmat
        common /data5/ xdmat
        common /data9/ bdot, biggam

        open(1,file='cntlr.crd',status='old')
        open(2,file='input.dat',status='old')
        open(3,file='ercox.out',status='unknown',access='append')

C       READ IN PROBLEM DIMENSIONS AND AUXILIARY INPUTS
        read(1,*) n, k, nb1, nb2, nq, icpr, ibpr
        nb = nb1 + nb2
        do l = 1, k
          read(1,*) (x(l,j), j = 1, nb1)
         enddo
        do ir = 1, k
          read(1,*) (a(ir,is), is = 1, k)
         enddo
        do iq = 1, nq
        do ir = 1, k
          read(1,*) (adot(iq,ir,is), is = 1, k)
         enddo
         enddo
        do iu = 1, nq
          read(1,*) (biggam(iu,iv), iv = 1, nq)
         enddo

C       PRINT HEADER
        write(3,*) 'COX MODEL WITH DISCRETE ERROR-PRONE COVARIATES ',
     $    'AND ARBITRARY ERROR-FREE COVARIATES'
        write(3,*)
        write(3,*) 'Number of records = ', n
        write(3,*) 'Number of error-prone covariates = ', nb1
        write(3,*)
     $    'Number of error-prone covariate configurations = ', k
        write(3,*) 'Number of error-free covariates = ', nb2
        write(3,*) 'Number of misclassication parameters = ', nq
        write(3,*)
        write(3,*) 'X MATRIX'
        do l = 1, k
          write(3,170) (x(l,j), j = 1, nb1)
         enddo
        write(3,*)
        write(3,*) 'A MATRIX'
        do ir = 1, k
          write(3,170) (a(ir,is), is = 1, k)
         enddo
        write(3,*)

C       COMPUTE B (INVERSE OF A)
        eps = 1.0d-6
        call gjr(a,b,k,eps,msing)
        if (msing .eq. 2) then
          write(3,*) 'Singular misclassification matrix'
          stop
         endif

C       CONTINUE INITIAL PRINTOUTS
        write(3,*) 'B MATRIX'
        do ir = 1, k
          write(3,170) (B(ir,is), is = 1, k)
         enddo
        write(3,*)
        write(3,*) 'ADOT MATRIX'
        do iq = 1, nq
        do ir = 1, k
          write(3,170) (adot(iq,ir,is), is = 1, k)
         enddo
         enddo
        write(3,*)
        write(3,*) 'BIGGAM MATRIX'
        do iu = 1, nq
          write(3,171) (biggam(iu,iv), iv = 1, nq)
         enddo
        write(3,*)
  170   format(10f8.4)
  171   format(7f11.7)
        write(3,*)

C       COMPUTE BDOT
        if (ibpr .eq. 1) write(3,*) 'BDOT MATRIX'
        do iq = 1, nq
          do ik1 = 1, k
          do ik2 = 1, k
            bdot(iq,ik1,ik2) = 0.0d0
            do ik3 = 1, k
            do ik4 = 1, k
              addtrm = -b(ik1,ik3)*adot(iq,ik3,ik4)*b(ik4,ik2)
              bdot(iq,ik1,ik2) = bdot(iq,ik1,ik2) + addtrm
             enddo
             enddo
            if (ibpr .eq. 1) write(3,10) iq, ik1, ik2, bdot(iq,ik1,ik2)
   10       format(3i3,f10.4)
           enddo
           enddo
         enddo
        if (ibpr .eq. 1) write(3,*)

C       COMPUTE SURROGATE X VALUES
        do ir = 1, k
          do l = 1, nb1
            xnmat(ir,l) = 0.0d0
            do is = 1, k
              xnmat(ir,l) = xnmat(ir,l) + b(ir,is)*x(is,l)
             enddo
            enddo
          enddo
        xnmor = xnmat

C       COMPUTE XDOT VALUES
        do iq = 1, nq
          do ir = 1, k
            do l = 1, nb1
              xdmat(iq,ir,l) = 0.0d0
              do is = 1, k
                addtrm = bdot(iq,ir,is)*x(is,l)
                xdmat(iq,ir,l) = xdmat(iq,ir,l) + addtrm
               enddo
              enddo
            enddo
         enddo

        do ir = 1, k
        do is = 1, k
          bidnt(ir,is) = 0.0d0
          bidnt(ir,ir) = 1.0d0
          borig(ir,is) = b(ir,is)
         enddo
         enddo

C       READ IN INPUT DATA:
C       SURV/CENS STATUS, ENTRY TIME, EXIT TIME, CONFIG, COVARIATES
        do i = 1, n
          read(2,*) del(i), t0(i), t(i), kv(i), (x2(i,l), l=1, nb2)
         enddo

C       SORT THE DATA AND IDENTIFY UNIQUE F-U TIMES
        call sort

C       CALL SUBROUTINE TO SOLVE SCORE EQUATIONS
        beta = 0.0d0
        b = bidnt
        xnmat = x
        call estsol(isol,beta,info1)
        call cscomp(beta,csc,nb)
        bet0 = beta(1)
        csc0 = csc(1)
        b = borig
        xnmat = xnmor
        call estsol(isol,beta,info2)

C       PRINT OUT RESULTS
        call cscomp(beta,csc,nb)
        call covcomp(beta,cov)
        write(3,*) 'BETA VECTOR AND STANDARD ERRORS'
        write(3,*)
        do ir = 1, nb
          sd = dsqrt(cov(ir,ir))
          write (3,200) ir, beta(ir), sd
  200     format(i4, 2F10.4)
         enddo
        write(3,*)
        write(3,*)
        write(3,*) 'FINAL ESTIMATION FUNCTION VALUES'
        write(3,*)
        write(3,*) (csc(l), l = 1, nb)
        write(3,*)
        write(3,*)
        if (icpr .eq. 1) then
          write(3,*) 'COVARIANCE MATRIX'
          write(3,*)
          do ir = 1, nb
          do is = 1, nb
            write(3,250) ir, is, cov(ir,is)
  250       format(2i4, F10.4)
           enddo
           enddo
           write(3,*)
           write(3,*)
         endif
        if (info2 .ne. 1) then
          write (3,*) 'Method failed'
          write(3,*)
          write(3,*)
         endif

        end

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

        subroutine estsol(isol,beta,info)

        implicit double precision(a-h, o-z)
        parameter (kmax=20, nbmax=20, nmax=90000, nrmax=300)
        external cjcomp

        dimension x(kmax,nbmax), b(kmax,kmax), x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(nbmax)
        dimension xnmat(kmax,nbmax)
        dimension cscor(nbmax), dqrjac(nbmax,nbmax), djac(nbmax,nbmax)
        dimension diag(nbmax), rmat(nrmax), qtf(nbmax)
        dimension wa1(nbmax), wa2(nbmax), wa3(nbmax), wa4(nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist
        common /data4/ xnmat

        nb = nb1 + nb2
        ldj = nbmax
        xtol = 0.01d0
        maxfev = 100
        mode = 1
        factor = 100.0d0
        nprint = -1
        lr = nrmax
        call hybrj(cjcomp,nb,beta,cscor,dqrjac,ldj,xtol,maxfev,diag,
     $    mode,factor,nprint,info,nfev,njev,rmat,lr,qtf,
     $    wa1,wa2,wa3,wa4)
        return
        end

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

        subroutine cjcomp(nvars,beta,cscor,djac,ldj,isw)
        implicit double precision (a-h,o-z)
        parameter (nbmax=20)
        dimension beta(nbmax), cscor(nbmax), djac(nbmax,nbmax)
        if (isw .eq. 1) call cscomp(beta,cscor,nvars)
        if (isw .eq. 2) call jaccmp(nvars,beta,djac)
        return
        end

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

        subroutine sort

        implicit double precision (a-h,o-z)
        parameter (nbmax=20,nmax=90000)

        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension indx(nmax), wksp(nmax), iwksp(nmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist

C       SORT DATA
        call qsort(indx,n,t)
        do i = 1, n
          wksp(i) = t(i)
         enddo
        do i = 1, n
          t(i) = wksp(indx(i))
         enddo
        do i = 1, n
          wksp(i) = del(i)
         enddo
        do i = 1, n
          del(i) = wksp(indx(i))
         enddo
        do i = 1, n
          iwksp(i) = kv(i)
         enddo
        do i = 1, n
          kv(i) = iwksp(indx(i))
         enddo
        do i = 1, n
          wksp(i) = t0(i)
         enddo
        do i = 1, n
          t0(i) = wksp(indx(i))
         enddo
        do j = 1, nb2
          do i = 1, n
            wksp(i) = x2(i,j)
           enddo
          do i = 1, n
            x2(i,j) = wksp(indx(i))
           enddo
         enddo

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

        return
        end

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

        subroutine cscomp(beta,cscor,nvars)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax), xnmat(kmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(nbmax), cscor(nbmax)
        dimension phival(nmax), psi1v(nmax,nbmax), dnum(nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist
        common /data4/ xnmat

        dn = n
        nb = nb1 + nb2

C        write(17,*) 'CSCOMP START'
C        write(17,*) (beta(l), l = 1, nb)

        cscor = 0.0d0

        do i = 1, n
          phival(i) = phi(i,beta)
          do l = 1, nb
            psi1v(i,l) = psi1(i,l,beta)
           enddo
          enddo

        do 90  m = 1, ndist

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

          tcurr = t(ilow(m))
          dden = 0.0d0
          dnum = 0.0d0

          do 65  j = ilow(m), n
            if (t0(j) .gt. tcurr) goto 65
            dden = dden + phival(j)
            do l = 1, nb
              dnum(l) = dnum(l) + psi1v(j,l)
             enddo
   65      continue

           do l = 1, nb
             addtrm = evcnt(m)*dnum(l)/dden
             cscor(l) = cscor(l) - addtrm
            enddo

           do i = ilow(m), iup(m)
             if (del(i) .eq. 1.0d0) then
               do l = 1, nb1
                 cscor(l) = cscor(l) + xnmat(kv(i),l)
                enddo
               do l = 1, nb2
                 cscor(l+nb1) = cscor(l+nb1) + x2(i,l)
                enddo
              endif
            enddo

   90    continue

        do l = 1, nb
          cscor(l) = -cscor(l)/dn
         enddo

C        write(17,*) 'CSCOMP END'
C        write(17,*) (cscor(l), l = 1, nb)

        return
        end

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

        subroutine jaccmp(nvars,beta,djac)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(nbmax)
        dimension phival(nmax), psi1v(nmax,nbmax)
        dimension psi2v(nmax,nbmax,nbmax)
        dimension dnum1(nbmax), dnum2(nbmax,nbmax)
        dimension djac(nbmax,nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist

        dn = n
        nb = nb1 + nb2

C        write(17,*) 'JACCMP START'
C        write(17,*) (beta(l), l = 1, nb)

        djac = 0.0d0


        do i = 1, n
          phival(i) = phi(i,beta)
          do l1 = 1, nb
            psi1v(i,l1) = psi1(i,l1,beta)
            do l2 = l1, nb
              psi2v(i,l1,l2) = psi2(i,l1,l2,beta)
             enddo
           enddo
         enddo

        do 90  m = 1, ndist

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

          tcurr = t(ilow(m))
          dden = 0.0d0
          dnum1 = 0.0d0
          dnum2 = 0.0d0

          do 65  j = ilow(m), n
            if (t0(j) .gt. tcurr) goto 65
            dden = dden + phival(j)
            do l1 = 1, nb
              dnum1(l1) = dnum1(l1) + psi1v(j,l1)
              do l2 = l1, nb
                dnum2(l1,l2) = dnum2(l1,l2) + psi2v(j,l1,l2)
               enddo
             enddo
   65      continue

           do l1 = 1, nb
           do l2 = l1, nb
             addtrm = (dnum2(l1,l2)/dden)
     $         - ((dnum1(l1)/dden)*(dnum1(l2)/dden))
             addtrm = evcnt(m)*addtrm
             djac(l1,l2) = djac(l1,l2) + addtrm
            enddo
            enddo

   90    continue

        do l1 = 1, nb
          djac(l1,l1) = djac(l1,l1)/dn
          do l2 = l1+1, nb
            djac(l1,l2) = djac(l1,l2)/dn
            djac(l2,l1) = djac(l1,l2)
           enddo
         enddo

C        write(17,*) 'JACCMP END'
C        do  l1 = 1, nb
C          write(17,*) l1, (djac(l1,l2), l2 = 1, nb)
C         enddo

        return
        end

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

        double precision function phi(i,beta)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension beta(nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2

c        write(6,*) 'PHI START', i, beta(1), beta(2)

        ir = kv(i)
        phi = 0.0d0

        x2scor = 0.0d0
        do l = 1, nb2
          x2scor = x2scor + beta(nb1+l)*x2(i,l)
         enddo

        do is = 1, k

          expterm1 = 0.0d0

          do l = 1, nb1
            expterm1 = expterm1 + beta(l)*x(is,l)
           enddo
          expterm1 = expterm1 + x2scor
          expterm1 = dexp(expterm1)

          phi = phi + b(ir,is)*expterm1

         enddo

c        write(6,*) 'PHI END', phi

        return
        end

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

        double precision function phidot(iq,i,beta)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000,nqmax=5)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension beta(nbmax)
        dimension bdot(nqmax,kmax,kmax), biggam(nqmax,nqmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data9/ bdot, biggam

        ir = kv(i)
        phidot = 0.0d0

        x2scor = 0.0d0
        do l = 1, nb2
          x2scor = x2scor + beta(nb1+l)*x2(i,l)
         enddo

        do is = 1, k

          expterm1 = 0.0d0

          do l = 1, nb1
            expterm1 = expterm1 + beta(l)*x(is,l)
           enddo
          expterm1 = expterm1 + x2scor
          expterm1 = dexp(expterm1)

          phidot = phidot + bdot(iq,ir,is)*expterm1

         enddo

        return
        end

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

        double precision function psi1(i,l,beta)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension beta(nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2

c        write(6,*) 'PSI1 START', i, beta(1), beta(2)

        ir = kv(i)
        psi1 = 0.0d0

        x2scor = 0.0d0
        do l0 = 1, nb2
          x2scor = x2scor + beta(nb1+l0)*x2(i,l0)
         enddo

        do is = 1, k

          expterm1 = 0.0d0
          do l1 = 1, nb1
            expterm1 = expterm1 + beta(l1)*x(is,l1)
           enddo
          expterm1 = expterm1 + x2scor
          expterm1 = dexp(expterm1)

          if (l .le. nb1)
     $      then
              xval = x(is,l)
            else
              xval = x2(i,l-nb1)
            endif

          psi1 = psi1 + b(ir,is)*xval*expterm1

         enddo

c        write(6,*) 'PSI1 END', l, psi1

        return
        end

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

        double precision function psi1dot(iq,i,l,beta)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000, nqmax=5)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension beta(nbmax), bdot(nqmax,kmax,kmax)
        dimension biggam(nqmax,nqmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data9/ bdot, biggam

        ir = kv(i)
        psi1dot = 0.0d0

        x2scor = 0.0d0
        do l0 = 1, nb2
          x2scor = x2scor + beta(nb1+l0)*x2(i,l0)
         enddo

        do is = 1, k

          expterm1 = 0.0d0
          do l1 = 1, nb1
            expterm1 = expterm1 + beta(l1)*x(is,l1)
           enddo
          expterm1 = expterm1 + x2scor
          expterm1 = dexp(expterm1)

          if (l .le. nb1)
     $      then
              xval = x(is,l)
            else
              xval = x2(i,l-nb1)
            endif

          psi1dot = psi1dot + bdot(iq,ir,is)*xval*expterm1

         enddo

        return
        end

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

        double precision function psi2(i,l1,l2,beta)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension beta(nbmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2

c        write(6,*) 'PSI2 START', i, beta(1), beta(2)

        ir = kv(i)
        psi2 = 0.0d0

        x2scor = 0.0d0
        do l0 = 1, nb2
          x2scor = x2scor + beta(nb1+l0)*x2(i,l0)
         enddo

        do is = 1, k

          expterm1 = 0.0d0
          do l1a = 1, nb1
            expterm1 = expterm1 + beta(l1a)*x(is,l1a)
           enddo
          expterm1 = expterm1 + x2scor
          expterm1 = dexp(expterm1)

          if (l1 .le. nb1)
     $      then
              xval1 = x(is,l1)
            else
              xval1 = x2(i,l1-nb1)
            endif

          if (l2 .le. nb1)
     $      then
              xval2 = x(is,l2)
            else
              xval2 = x2(i,l2-nb1)
            endif

          psi2 = psi2 + b(ir,is)*xval1*xval2*expterm1

         enddo

c         write(6,*) 'PSI2 END', l1, l2, psi2

         return
         end

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

        subroutine covcomp(beta,cov)

        implicit double precision (a-h,o-z)
        parameter (kmax=20, nbmax=20, nmax=90000, nqmax=5)

        dimension x(kmax,nbmax), b(kmax,kmax)
        dimension x2(nmax,nbmax), xnmat(kmax,nbmax)
        dimension kv(nmax), del(nmax), t0(nmax), t(nmax)
        dimension evcnt(nmax), ilow(nmax), iup(nmax)
        dimension beta(nbmax), cov(nbmax,nbmax)
        dimension varmat(nbmax,nbmax), plmat1(nbmax,nbmax)
        dimension dnum(nbmax), phival(nmax), psi1v(nmax,nbmax)
        dimension djac(nbmax,nbmax), djinv(nbmax,nbmax)
        dimension xdmat(nqmax,kmax,nbmax), bdot(nqmax,kmax,kmax)
        dimension biggam(nqmax,nqmax)
        dimension pvdot(nqmax,nmax), ps1dot(nqmax,nmax,nbmax)
        dimension udot(nbmax,nqmax), ddendt(nqmax), dnd(nqmax,nbmax)
        dimension abder(nbmax,nqmax)

        common /data0/ n, k, nb1, nb2, nq
        common /data1/ x, b
        common /data2/ kv, del, t0, t, x2
        common /data3/ evcnt, ilow, iup, ndist
        common /data4/ xnmat
        common /data5/ xdmat
        common /data9/ bdot, biggam

        dn = n
        nb = nb1 + nb2
        varmat = 0.0d0
        udot = 0.0d0

        do i = 1, n
          phival(i) = phi(i,beta)
        do l1 = 1, nb
          psi1v(i,l1) = psi1(i,l1,beta)
         enddo
         enddo

        do iq = 1, nq
          do i = 1, n
            pvdot(iq,i) = phidot(iq,i,beta)
          do l1 = 1, nb
            ps1dot(iq,i,l1) = psi1dot(iq,i,l1,beta)
           enddo
           enddo
         enddo

        do 80  m = 1, ndist

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

          tcurr = t(ilow(m))
          dden = 0.0d0
          ddendt = 0.0d0
          dnum = 0.0d0
          dnd = 0.0d0

          do 65  j = ilow(m), n
            if (t0(j) .gt. tcurr) goto 65
            dden = dden + phival(j)
            do iq = 1, nq
              ddendt(iq) = ddendt(iq) + pvdot(iq,j)
             enddo
            do l1 = 1, nb
              dnum(l1) = dnum(l1) + psi1v(j,l1)
              do iq = 1, nq
                dnd(iq,l1) = dnd(iq,l1) + ps1dot(iq,j,l1)
               enddo
             enddo
   65      continue


          do iq = 1, nq
          do l = 1, nb
            addtrm = - (dnd(iq,l)/dden)
     $        + (dnum(l)/dden)*(ddendt(iq)/dden)
            udot(l,iq) = udot(l,iq) + evcnt(m)*addtrm/dn
           enddo
           enddo

           do 75  i = ilow(m), iup(m)
             if (del(i) .eq. 0.0d0) goto 75
             do l1 = 1, nb
               drat1 = dnum(l1)/dden
               if (l1 .le. nb1)
     $           then
                   xnew1 = xnmat(kv(i),l1)
                 else
                   xnew1 = x2(i,l1-nb1)
                 endif
               dif1 = xnew1 - drat1
               do l2 = l1, nb
                 drat2 = dnum(l2)/dden
                 if (l2 .le. nb1)
     $             then
                     xnew2 = xnmat(kv(i),l2)
                   else
                     xnew2 = x2(i,l2-nb1)
                   endif
                 dif2 = xnew2 - drat2
                 varmat(l1,l2) = varmat(l1,l2) + dif1*dif2
                enddo
               if (l1 .le. nb1) then
                 do iq = 1, nq
                   udot(l1,iq) = udot(l1,iq) + xdmat(iq,kv(i),l1)/dn
                  enddo
                endif
              enddo
   75       continue

   80    continue

        do l1 = 1, nb
        do l2 = l1, nb
          do iq1 = 1, nq
          do iq2 = 1, nq
            addtrm = (dn**2)*udot(l1,iq1)*biggam(iq1,iq2)*udot(l2,iq2)
            varmat(l1,l2) = varmat(l1,l2) + addtrm
           enddo
           enddo
          enddo
          enddo

        do l1 = 1, nb
          varmat(l1,l1) = varmat(l1,l1)/dn
          do l2 = l1+1, nb
            varmat(l1,l2) = varmat(l1,l2)/dn
            varmat(l2,l1) = varmat(l1,l2)
           enddo
         enddo

        call jaccmp(nb,beta,djac)
        eps = 1.0d-6
        call gjr(djac,djinv,nb,eps,msing)
        call matmul(nb,nb,nb,djinv,varmat,plmat1)
        call matmul(nb,nb,nb,plmat1,djinv,cov)

        do l1 = 1, nb
          cov(l1,l1) = cov(l1,l1)/dn
          do l2 = l1+1, nb
            cov(l1,l2) = cov(l1,l2)/dn
            cov(l2,l1) = cov(l1,l2)
           enddo
         enddo

c       COMPUTE DERIVATIVES OF BETA-HAT W.R.T. A
c        if (nq .eq. 0) return
c        write(3,*) 'UDOT'
c        do l = 1, nb
c          write(3,110) (udot(l,iq), iq=1,nq)
c         enddo
c        call dmrrrr(nb,nb,djinv,nbmax,nb,nq,udot,nqmax,
c     $    nb,nq,abder,nbmax)
c        write(3,*) 'ABDER'
c        do l = 1, nb
c          write(3,110) (abder(l,iq), iq=1,nq)
c  110     format(2f15.4)
c         enddo
c        write(3,*)

        return
        end

C*********************************************************************
C
C       END OF MAIN SUBROUTINES
C       THE REST ARE UTILITY SUBROUTINES, AS FOLLOWS
C       matmul: matrix multiplication
C       gjr: matrix inversion
C       qsort: quicksort
C       hybrj: solve a nonlinear system via Powell's hybrid method
C       the rest (dogleg, dpmpar, enorm, qform,qrfac,r1mpyq,r1updt):
C         minpack routines used in hybrj
C
C*********************************************************************

        subroutine matmul(nr1,nc1,nc2,a1,a2,a3)
        parameter (nr1m=20, nc1m=20, nc2m=20)
        double precision a1(nr1m,nc1m),a2(nc1m,nc2m),a3(nr1m,nc2m)
        do i = 1, nr1
        do j = 1, nc2
          a3(i,j) = 0.0d0
          do l = 1, nc1
            a3(i,j) = a3(i,j) + a1(i,l)*a2(l,j)
           enddo
         enddo
         enddo
        return
        end

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

      SUBROUTINE  GJR(AIN,A,N,EPS,MSING)
      parameter (MM=20, NN=20)
C     GAUSS-JORDAN-RUTISHAUSER MATRIX INVERSION WITH DOUBLE PIVOTING.
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      DIMENSION AIN(MM,NN),A(MM,NN),B(50),C(50),P(50),Q(50)
      INTEGER P,Q
C      IF(N.GT.100) stop
C     COPY INPUT MATRIX INTO WORKING/OUTPUT MATRIX
      DO 5 I=1,N
      DO 5 J=1,N
 5    A(I,J)=AIN(I,J)
      MSING=1
      DO 10 K=1,N
C     DETERMINATION OF THE PIVOT ELEMENT
      PIVOT=0.
      DO 20 I=K,N
      DO 20 J=K,N
      IF(ABS(A(I,J))-ABS(PIVOT))20,20,30
   30 PIVOT=A(I,J)
      P(K)=I
      Q(K)=J
   20 CONTINUE
      IF(ABS(PIVOT)-EPS)40,40,50
C     EXCHANGE OF THE PIVOTAL ROW WITH THE KTH ROW
   50 IF(P(K)-K)60,80,60
   60 DO 70 J=1,N
      L=P(K)
      Z=A(L,J)
      A(L,J)=A(K,J)
   70 A(K,J)=Z
C     EXCHANGE OF THE PIVOTAL COLUMN WITH THE KTH COLUMN
   80 IF(Q(K)-K)85,90,85
   85 DO 100 I=1,N
      L=Q(K)
      Z=A(I,L)
      A(I,L)=A(I,K)
  100 A(I,K)=Z
  90  CONTINUE
C     JORDAN STEP
      DO 110 J=1,N
      IF(J-K)130,120,130
  120 B(J)=1./PIVOT
      C(J)=1.
      GO TO 140
  130 B(J)=-A(K,J)/PIVOT
      C(J)=A(J,K)
  140 A(K,J)=0.
  110 A(J,K)=0.
      DO 10 I=1,N
      DO 10 J=1,N
   10 A(I,J)=A(I,J)+C(I)*B(J)
C     REORDERING THE MATRIX
      DO 155 M=1,N
      K=N-M+1
      IF(P(K)-K)160,170,160
  160 DO 180 I=1,N
      L=P(K)
      Z=A(I,L)
      A(I,L)=A(I,K)
  180 A(I,K)=Z
  170 IF(Q(K)-K)190,155,190
  190 DO 150 J=1,N
      L=Q(K)
      Z=A(L,J)
      A(L,J)=A(K,J)
  150 A(K,J)=Z
  155 CONTINUE
      RETURN
   40 CONTINUE
      MSING=2
      RETURN
      END

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

C From HDK@psuvm.psu.edu Thu Dec  8 15:27:16 MST 1994
C
C The following was converted from Algol recursive to Fortran iterative
C by a colleague at Penn State (a long time ago - Fortran 66, please
C excuse the GoTo's). The following code also corrects a bug in the
C Quicksort algorithm published in the ACM (see Algorithm 402, CACM,
C Sept. 1970, pp 563-567; also you younger folks who weren't born at
C that time might find interesting the history of the Quicksort
C algorithm beginning with the original published in CACM, July 1961,
C pp 321-322, Algorithm 64). Note that the following algorithm sorts
C integer data; actual data is not moved but sort is affected by sorting
C a companion index array (see leading comments). The data type being
C sorted can be changed by changing one line; see comments after
C declarations and subsequent one regarding comparisons(Fortran
C 77 takes care of character comparisons of course, so that comment
C is merely historical from the days when we had to write character
C compare subprograms, usually in assembler language for a specific
C mainframe platform at that time). But the following algorithm is
C good, still one of the best available.


      SUBROUTINE QSORT(ORD,N,A)
C
C==============SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
C   ASCENDING ORDER VECTOR IN ORD.  THAT IS ASCENDING ORDERED A
C   IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
C   I=1,2,...,N .  THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .
C
C
C     ACM QUICKSORT - ALGORITHM #402 - IMPLEMENTED IN FORTRAN 66 BY
C                                 WILLIAM H. VERITY, WHV@PSUVM.PSU.EDU
C                                 CENTER FOR ACADEMIC COMPUTING
C                                 THE PENNSYLVANIA STATE UNIVERSITY
C                                 UNIVERSITY PARK, PA.  16802
C
      IMPLICIT INTEGER (A-Z)
C
      DIMENSION ORD(N),POPLST(2,20)
      double precision X,XX,Z,ZZ,Y
C
C     TO SORT DIFFERENT INPUT TYPES, CHANGE THE FOLLOWING
C     SPECIFICATION STATEMENTS; FOR EXAMPLE, FOR FORTRAN CHARACTER
C     USE THE FOLLOWING:  CHARACTER *(*) A(N)
C
      double precision A(N)
C
      NDEEP=0
      U1=N
      L1=1
      DO 1  I=1,N
    1 ORD(I)=I
    2 IF (U1.LE.L1) RETURN
C
    3 L=L1
      U=U1
C
C PART
C
    4 P=L
      Q=U
C     FOR CHARACTER SORTS, THE FOLLOWING 3 STATEMENTS WOULD BECOME
C     X = ORD(P)
C     Z = ORD(Q)
C     IF (A(X) .LE. A(Z)) GO TO 2
C
C     WHERE "CLE" IS A LOGICAL FUNCTION WHICH RETURNS "TRUE" IF THE
C     FIRST ARGUMENT IS LESS THAN OR EQUAL TO THE SECOND, BASED ON "LEN"
C     CHARACTERS.
C
      X=A(ORD(P))
      Z=A(ORD(Q))
      IF (X.LE.Z) GO TO 5
      Y=X
      X=Z
      Z=Y
      YP=ORD(P)
      ORD(P)=ORD(Q)
      ORD(Q)=YP
    5 IF (U-L.LE.1) GO TO 15
      XX=X
      IX=P
      ZZ=Z
      IZ=Q
C
C LEFT
C
    6 P=P+1
      IF (P.GE.Q) GO TO 7
      X=A(ORD(P))
      IF (X.GE.XX) GO TO 8
      GO TO 6
    7 P=Q-1
      GO TO 13
C
C RIGHT
C
    8 Q=Q-1
      IF (Q.LE.P) GO TO 9
      Z=A(ORD(Q))
      IF (Z.LE.ZZ) GO TO 10
      GO TO 8
    9 Q=P
      P=P-1
      Z=X
      X=A(ORD(P))
C
C DIST
C
   10 IF (X.LE.Z) GO TO 11
      Y=X
      X=Z
      Z=Y
      IP=ORD(P)
      ORD(P)=ORD(Q)
      ORD(Q)=IP
   11 IF (X.LE.XX) GO TO 12
      XX=X
      IX=P
   12 IF (Z.GE.ZZ) GO TO 6
      ZZ=Z
      IZ=Q
      GO TO 6
C
C OUT
C
   13 CONTINUE
      IF (.NOT.(P.NE.IX.AND.X.NE.XX)) GO TO 14
      IP=ORD(P)
      ORD(P)=ORD(IX)
      ORD(IX)=IP
   14 CONTINUE
      IF (.NOT.(Q.NE.IZ.AND.Z.NE.ZZ)) GO TO 15
      IQ=ORD(Q)
      ORD(Q)=ORD(IZ)
      ORD(IZ)=IQ
   15 CONTINUE
      IF (U-Q.LE.P-L) GO TO 16
      L1=L
      U1=P-1
      L=Q+1
      GO TO 17
   16 U1=U
      L1=Q+1
      U=P-1
   17 CONTINUE
      IF (U1.LE.L1) GO TO 18
C
C START RECURSIVE CALL
C
      NDEEP=NDEEP+1
      POPLST(1,NDEEP)=U
      POPLST(2,NDEEP)=L
      GO TO 3
   18 IF (U.GT.L) GO TO 4
C
C POP BACK UP IN THE RECURSION LIST
C
      IF (NDEEP.EQ.0) GO TO 2
      U=POPLST(1,NDEEP)
      L=POPLST(2,NDEEP)
      NDEEP=NDEEP-1
      GO TO 18
C
C END SORT
C END QSORT
C
      END

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

      subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
      integer n,lr
      double precision delta
      double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
c     **********
c
c     subroutine dogleg
c
c     given an m by n matrix a, an n by n nonsingular diagonal
c     matrix d, an m-vector b, and a positive number delta, the
c     problem is to determine the convex combination x of the
c     gauss-newton and scaled gradient directions that minimizes
c     (a*x - b) in the least squares sense, subject to the
c     restriction that the euclidean norm of d*x be at most delta.
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization of a. that is, if a = q*r, where q has
c     orthogonal columns and r is an upper triangular matrix,
c     then dogleg expects the full upper triangle of r and
c     the first n components of (q transpose)*b.
c
c     the subroutine statement is
c
c       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an input array of length lr which must contain the upper
c         triangular matrix r stored by rows.
c
c       lr is a positive integer input variable not less than
c         (n*(n+1))/2.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       delta is a positive input variable which specifies an upper
c         bound on the euclidean norm of d*x.
c
c       x is an output array of length n which contains the desired
c         convex combination of the gauss-newton direction and the
c         scaled gradient direction.
c
c       wa1 and wa2 are work arrays of length n.
c
c     subprograms called
c
c       minpack-supplied ... dpmpar,enorm
c
c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,j,jj,jp1,k,l
      double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
     *                 temp,zero
      double precision dpmpar,enorm
      data one,zero /1.0d0,0.0d0/
c
c     epsmch is the machine precision.
c
      epsmch = dpmpar(1)
c
c     first, calculate the gauss-newton direction.
c
      jj = (n*(n + 1))/2 + 1
      do 50 k = 1, n
         j = n - k + 1
         jp1 = j + 1
         jj = jj - k
         l = jj + 1
         sum = zero
         if (n .lt. jp1) go to 20
         do 10 i = jp1, n
            sum = sum + r(l)*x(i)
            l = l + 1
   10       continue
   20    continue
         temp = r(jj)
         if (temp .ne. zero) go to 40
         l = j
         do 30 i = 1, j
            temp = dmax1(temp,dabs(r(l)))
            l = l + n - i
   30       continue
         temp = epsmch*temp
         if (temp .eq. zero) temp = epsmch
   40    continue
         x(j) = (qtb(j) - sum)/temp
   50    continue
c
c     test whether the gauss-newton direction is acceptable.
c
      do 60 j = 1, n
         wa1(j) = zero
         wa2(j) = diag(j)*x(j)
   60    continue
      qnorm = enorm(n,wa2)
      if (qnorm .le. delta) go to 140
c
c     the gauss-newton direction is not acceptable.
c     next, calculate the scaled gradient direction.
c
      l = 1
      do 80 j = 1, n
         temp = qtb(j)
         do 70 i = j, n
            wa1(i) = wa1(i) + r(l)*temp
            l = l + 1
   70       continue
         wa1(j) = wa1(j)/diag(j)
   80    continue
c
c     calculate the norm of the scaled gradient and test for
c     the special case in which the scaled gradient is zero.
c
      gnorm = enorm(n,wa1)
      sgnorm = zero
      alpha = delta/qnorm
      if (gnorm .eq. zero) go to 120
c
c     calculate the point along the scaled gradient
c     at which the quadratic is minimized.
c
      do 90 j = 1, n
         wa1(j) = (wa1(j)/gnorm)/diag(j)
   90    continue
      l = 1
      do 110 j = 1, n
         sum = zero
         do 100 i = j, n
            sum = sum + r(l)*wa1(i)
            l = l + 1
  100       continue
         wa2(j) = sum
  110    continue
      temp = enorm(n,wa2)
      sgnorm = (gnorm/temp)/temp
c
c     test whether the scaled gradient direction is acceptable.
c
      alpha = zero
      if (sgnorm .ge. delta) go to 120
c
c     the scaled gradient direction is not acceptable.
c     finally, calculate the point along the dogleg
c     at which the quadratic is minimized.
c
      bnorm = enorm(n,qtb)
      temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
      temp = temp - (delta/qnorm)*(sgnorm/delta)**2
     *       + dsqrt((temp-(delta/qnorm))**2
     *               +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
      alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
  120 continue
c
c     form appropriate convex combination of the gauss-newton
c     direction and the scaled gradient direction.
c
      temp = (one - alpha)*dmin1(sgnorm,delta)
      do 130 j = 1, n
         x(j) = temp*wa1(j) + alpha*x(j)
  130    continue
  140 continue
      return
c
c     last card of subroutine dogleg.
c
      end

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

      double precision function dpmpar(i)
      integer i
c     **********
c
c     Function dpmpar
c
c     This function provides double precision machine parameters
c     when the appropriate set of data statements is activated (by
c     removing the c from column 1) and all other data statements are
c     rendered inactive. Most of the parameter values were obtained
c     from the corresponding Bell Laboratories Port Library function.
c
c     The function statement is
c
c       double precision function dpmpar(i)
c
c     where
c
c       i is an integer input variable set to 1, 2, or 3 which
c         selects the desired machine parameter. If the machine has
c         t base b digits and its smallest and largest exponents are
c         emin and emax, respectively, then these parameters are
c
c         dpmpar(1) = b**(1 - t), the machine precision,
c
c         dpmpar(2) = b**(emin - 1), the smallest magnitude,
c
c         dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
c
c     Argonne National Laboratory. MINPACK Project. November 1996.
c     Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
c
c     **********
      integer mcheps(4)
      integer minmag(4)
      integer maxmag(4)
      double precision dmach(3)
      equivalence (dmach(1),mcheps(1))
      equivalence (dmach(2),minmag(1))
      equivalence (dmach(3),maxmag(1))
c
c     Machine constants for the IBM 360/370 series,
c     the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
c     the Xerox Sigma 5/7/9 and the Sel systems 85/86.
c
c     data mcheps(1),mcheps(2) / z34100000, z00000000 /
c     data minmag(1),minmag(2) / z00100000, z00000000 /
c     data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
c
c     Machine constants for the Honeywell 600/6000 series.
c
c     data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
c     data minmag(1),minmag(2) / o402400000000, o000000000000 /
c     data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
c
c     Machine constants for the CDC 6000/7000 series.
c
c     data mcheps(1) / 15614000000000000000b /
c     data mcheps(2) / 15010000000000000000b /
c
c     data minmag(1) / 00604000000000000000b /
c     data minmag(2) / 00000000000000000000b /
c
c     data maxmag(1) / 37767777777777777777b /
c     data maxmag(2) / 37167777777777777777b /
c
c     Machine constants for the PDP-10 (KA processor).
c
c     data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
c     data minmag(1),minmag(2) / "033400000000, "000000000000 /
c     data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
c
c     Machine constants for the PDP-10 (KI processor).
c
c     data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
c     data minmag(1),minmag(2) / "000400000000, "000000000000 /
c     data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
c
c     Machine constants for the PDP-11.
c
c     data mcheps(1),mcheps(2) /   9472,      0 /
c     data mcheps(3),mcheps(4) /      0,      0 /
c
c     data minmag(1),minmag(2) /    128,      0 /
c     data minmag(3),minmag(4) /      0,      0 /
c
c     data maxmag(1),maxmag(2) /  32767,     -1 /
c     data maxmag(3),maxmag(4) /     -1,     -1 /
c
c     Machine constants for the Burroughs 6700/7700 systems.
c
c     data mcheps(1) / o1451000000000000 /
c     data mcheps(2) / o0000000000000000 /
c
c     data minmag(1) / o1771000000000000 /
c     data minmag(2) / o7770000000000000 /
c
c     data maxmag(1) / o0777777777777777 /
c     data maxmag(2) / o7777777777777777 /
c
c     Machine constants for the Burroughs 5700 system.
c
c     data mcheps(1) / o1451000000000000 /
c     data mcheps(2) / o0000000000000000 /
c
c     data minmag(1) / o1771000000000000 /
c     data minmag(2) / o0000000000000000 /
c
c     data maxmag(1) / o0777777777777777 /
c     data maxmag(2) / o0007777777777777 /
c
c     Machine constants for the Burroughs 1700 system.
c
c     data mcheps(1) / zcc6800000 /
c     data mcheps(2) / z000000000 /
c
c     data minmag(1) / zc00800000 /
c     data minmag(2) / z000000000 /
c
c     data maxmag(1) / zdffffffff /
c     data maxmag(2) / zfffffffff /
c
c     Machine constants for the Univac 1100 series.
c
c     data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
c     data minmag(1),minmag(2) / o000040000000, o000000000000 /
c     data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
c
c     Machine constants for the Data General Eclipse S/200.
c
c     Note - it may be appropriate to include the following card -
c     static dmach(3)
c
c     data minmag/20k,3*0/,maxmag/77777k,3*177777k/
c     data mcheps/32020k,3*0/
c
c     Machine constants for the Harris 220.
c
c     data mcheps(1),mcheps(2) / '20000000, '00000334 /
c     data minmag(1),minmag(2) / '20000000, '00000201 /
c     data maxmag(1),maxmag(2) / '37777777, '37777577 /
c
c     Machine constants for the Cray-1.
c
c     data mcheps(1) / 0376424000000000000000b /
c     data mcheps(2) / 0000000000000000000000b /
c
c     data minmag(1) / 0200034000000000000000b /
c     data minmag(2) / 0000000000000000000000b /
c
c     data maxmag(1) / 0577777777777777777777b /
c     data maxmag(2) / 0000007777777777777776b /
c
c     Machine constants for the Prime 400.
c
c     data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
c     data minmag(1),minmag(2) / :10000000000, :00000100000 /
c     data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
c
c     Machine constants for the VAX-11.
c
c     data mcheps(1),mcheps(2) /   9472,  0 /
c     data minmag(1),minmag(2) /    128,  0 /
c     data maxmag(1),maxmag(2) / -32769, -1 /
c
c     Machine constants for IEEE machines.
c
      data dmach(1) /2.22044604926d-16/
      data dmach(2) /2.22507385852d-308/
      data dmach(3) /1.79769313485d+308/
c
      dpmpar = dmach(i)
      return
c
c     Last card of function dpmpar.
c
      end

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

      double precision function enorm(n,x)
      integer n
      double precision x(n)
c     **********
c
c     function enorm
c
c     given an n-vector x, this function calculates the
c     euclidean norm of x.
c
c     the euclidean norm is computed by accumulating the sum of
c     squares in three different sums. the sums of squares for the
c     small and large components are scaled so that no overflows
c     occur. non-destructive underflows are permitted. underflows
c     and overflows do not occur in the computation of the unscaled
c     sum of squares for the intermediate components.
c     the definitions of small, intermediate and large components
c     depend on two constants, rdwarf and rgiant. the main
c     restrictions on these constants are that rdwarf**2 not
c     underflow and rgiant**2 not overflow. the constants
c     given here are suitable for every known computer.
c
c     the function statement is
c
c       double precision function enorm(n,x)
c
c     where
c
c       n is a positive integer input variable.
c
c       x is an input array of length n.
c
c     subprograms called
c
c       fortran-supplied ... dabs,dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i
      double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,
     *                 x1max,x3max,zero
      data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
      s1 = zero
      s2 = zero
      s3 = zero
      x1max = zero
      x3max = zero
      floatn = n
      agiant = rgiant/floatn
      do 90 i = 1, n
         xabs = dabs(x(i))
         if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
            if (xabs .le. rdwarf) go to 30
c
c              sum for large components.
c
               if (xabs .le. x1max) go to 10
                  s1 = one + s1*(x1max/xabs)**2
                  x1max = xabs
                  go to 20
   10          continue
                  s1 = s1 + (xabs/x1max)**2
   20          continue
               go to 60
   30       continue
c
c              sum for small components.
c
               if (xabs .le. x3max) go to 40
                  s3 = one + s3*(x3max/xabs)**2
                  x3max = xabs
                  go to 50
   40          continue
                  if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
   50          continue
   60       continue
            go to 80
   70    continue
c
c           sum for intermediate components.
c
            s2 = s2 + xabs**2
   80    continue
   90    continue
c
c     calculation of norm.
c
      if (s1 .eq. zero) go to 100
         enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
         go to 130
  100 continue
         if (s2 .eq. zero) go to 110
            if (s2 .ge. x3max)
     *         enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
            if (s2 .lt. x3max)
     *         enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
            go to 120
  110    continue
            enorm = x3max*dsqrt(s3)
  120    continue
  130 continue
      return
c
c     last card of function enorm.
c
      end

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

      subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode,
     *                 factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,
     *                 wa3,wa4)
      integer n,ldfjac,maxfev,mode,nprint,info,nfev,njev,lr
      double precision xtol,factor
      double precision x(n),fvec(n),fjac(ldfjac,n),diag(n),r(lr),
     *                 qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
c     **********
c
c     subroutine hybrj
c
c     the purpose of hybrj is to find a zero of a system of
c     n nonlinear functions in n variables by a modification
c     of the powell hybrid method. the user must provide a
c     subroutine which calculates the functions and the jacobian.
c
c     the subroutine statement is
c
c       subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,
c                        mode,factor,nprint,info,nfev,njev,r,lr,qtf,
c                        wa1,wa2,wa3,wa4)
c
c     where
c
c       fcn is the name of the user-supplied subroutine which
c         calculates the functions and the jacobian. fcn must
c         be declared in an external statement in the user
c         calling program, and should be written as follows.
c
c         subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
c         integer n,ldfjac,iflag
c         double precision x(n),fvec(n),fjac(ldfjac,n)
c         ----------
c         if iflag = 1 calculate the functions at x and
c         return this vector in fvec. do not alter fjac.
c         if iflag = 2 calculate the jacobian at x and
c         return this matrix in fjac. do not alter fvec.
c         ---------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless
c         the user wants to terminate execution of hybrj.
c         in this case set iflag to a negative integer.
c
c       n is a positive integer input variable set to the number
c         of functions and variables.
c
c       x is an array of length n. on input x must contain
c         an initial estimate of the solution vector. on output x
c         contains the final estimate of the solution vector.
c
c       fvec is an output array of length n which contains
c         the functions evaluated at the output x.
c
c       fjac is an output n by n array which contains the
c         orthogonal matrix q produced by the qr factorization
c         of the final approximate jacobian.
c
c       ldfjac is a positive integer input variable not less than n
c         which specifies the leading dimension of the array fjac.
c
c       xtol is a nonnegative input variable. termination
c         occurs when the relative error between two consecutive
c         iterates is at most xtol.
c
c       maxfev is a positive integer input variable. termination
c         occurs when the number of calls to fcn with iflag = 1
c         has reached maxfev.
c
c       diag is an array of length n. if mode = 1 (see
c         below), diag is internally set. if mode = 2, diag
c         must contain positive entries that serve as
c         multiplicative scale factors for the variables.
c
c       mode is an integer input variable. if mode = 1, the
c         variables will be scaled internally. if mode = 2,
c         the scaling is specified by the input diag. other
c         values of mode are equivalent to mode = 1.
c
c       factor is a positive input variable used in determining the
c         initial step bound. this bound is set to the product of
c         factor and the euclidean norm of diag*x if nonzero, or else
c         to factor itself. in most cases factor should lie in the
c         interval (.1,100.). 100. is a generally recommended value.
c
c       nprint is an integer input variable that enables controlled
c         printing of iterates if it is positive. in this case,
c         fcn is called with iflag = 0 at the beginning of the first
c         iteration and every nprint iterations thereafter and
c         immediately prior to return, with x and fvec available
c         for printing. fvec and fjac should not be altered.
c         if nprint is not positive, no special calls of fcn
c         with iflag = 0 are made.
c
c       info is an integer output variable. if the user has
c         terminated execution, info is set to the (negative)
c         value of iflag. see description of fcn. otherwise,
c         info is set as follows.
c
c         info = 0   improper input parameters.
c
c         info = 1   relative error between two consecutive iterates
c                    is at most xtol.
c
c         info = 2   number of calls to fcn with iflag = 1 has
c                    reached maxfev.
c
c         info = 3   xtol is too small. no further improvement in
c                    the approximate solution x is possible.
c
c         info = 4   iteration is not making good progress, as
c                    measured by the improvement from the last
c                    five jacobian evaluations.
c
c         info = 5   iteration is not making good progress, as
c                    measured by the improvement from the last
c                    ten iterations.
c
c       nfev is an integer output variable set to the number of
c         calls to fcn with iflag = 1.
c
c       njev is an integer output variable set to the number of
c         calls to fcn with iflag = 2.
c
c       r is an output array of length lr which contains the
c         upper triangular matrix produced by the qr factorization
c         of the final approximate jacobian, stored rowwise.
c
c       lr is a positive integer input variable not less than
c         (n*(n+1))/2.
c
c       qtf is an output array of length n which contains
c         the vector (q transpose)*fvec.
c
c       wa1, wa2, wa3, and wa4 are work arrays of length n.
c
c     subprograms called
c
c       user-supplied ...... fcn
c
c       minpack-supplied ... dogleg,dpmpar,enorm,
c                            qform,qrfac,r1mpyq,r1updt
c
c       fortran-supplied ... dabs,dmax1,dmin1,mod
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,iflag,iter,j,jm1,l,ncfail,ncsuc,nslow1,nslow2
      integer iwa(1)
      logical jeval,sing
      double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
     *                 prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
     *                 zero
      double precision dpmpar,enorm
      data one,p1,p5,p001,p0001,zero
     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
c
c     epsmch is the machine precision.
c
      epsmch = dpmpar(1)
c
      info = 0
      iflag = 0
      nfev = 0
      njev = 0
c
c     check the input parameters for errors.
c
      if (n .le. 0 .or. ldfjac .lt. n .or. xtol .lt. zero
     *    .or. maxfev .le. 0 .or. factor .le. zero
     *    .or. lr .lt. (n*(n + 1))/2) go to 300
      if (mode .ne. 2) go to 20
      do 10 j = 1, n
         if (diag(j) .le. zero) go to 300
   10    continue
   20 continue
c
c     evaluate the function at the starting point
c     and calculate its norm.
c
      iflag = 1
      call fcn(n,x,fvec,fjac,ldfjac,iflag)
      nfev = 1
      if (iflag .lt. 0) go to 300
      fnorm = enorm(n,fvec)
c
c     initialize iteration counter and monitors.
c
      iter = 1
      ncsuc = 0
      ncfail = 0
      nslow1 = 0
      nslow2 = 0
c
c     beginning of the outer loop.
c
   30 continue
         jeval = .true.
c
c        calculate the jacobian matrix.
c
         iflag = 2
         call fcn(n,x,fvec,fjac,ldfjac,iflag)
         njev = njev + 1
         if (iflag .lt. 0) go to 300
c
c        compute the qr factorization of the jacobian.
c
         call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
c
c        on the first iteration and if mode is 1, scale according
c        to the norms of the columns of the initial jacobian.
c
         if (iter .ne. 1) go to 70
         if (mode .eq. 2) go to 50
         do 40 j = 1, n
            diag(j) = wa2(j)
            if (wa2(j) .eq. zero) diag(j) = one
   40       continue
   50    continue
c
c        on the first iteration, calculate the norm of the scaled x
c        and initialize the step bound delta.
c
         do 60 j = 1, n
            wa3(j) = diag(j)*x(j)
   60       continue
         xnorm = enorm(n,wa3)
         delta = factor*xnorm
         if (delta .eq. zero) delta = factor
   70    continue
c
c        form (q transpose)*fvec and store in qtf.
c
         do 80 i = 1, n
            qtf(i) = fvec(i)
   80       continue
         do 120 j = 1, n
            if (fjac(j,j) .eq. zero) go to 110
            sum = zero
            do 90 i = j, n
               sum = sum + fjac(i,j)*qtf(i)
   90          continue
            temp = -sum/fjac(j,j)
            do 100 i = j, n
               qtf(i) = qtf(i) + fjac(i,j)*temp
  100          continue
  110       continue
  120       continue
c
c        copy the triangular factor of the qr factorization into r.
c
         sing = .false.
         do 150 j = 1, n
            l = j
            jm1 = j - 1
            if (jm1 .lt. 1) go to 140
            do 130 i = 1, jm1
               r(l) = fjac(i,j)
               l = l + n - i
  130          continue
  140       continue
            r(l) = wa1(j)
            if (wa1(j) .eq. zero) sing = .true.
  150       continue
c
c        accumulate the orthogonal factor in fjac.
c
         call qform(n,n,fjac,ldfjac,wa1)
c
c        rescale if necessary.
c
         if (mode .eq. 2) go to 170
         do 160 j = 1, n
            diag(j) = dmax1(diag(j),wa2(j))
  160       continue
  170    continue
c
c        beginning of the inner loop.
c
  180    continue
c
c           if requested, call fcn to enable printing of iterates.
c
            if (nprint .le. 0) go to 190
            iflag = 0
            if (mod(iter-1,nprint) .eq. 0)
     *         call fcn(n,x,fvec,fjac,ldfjac,iflag)
            if (iflag .lt. 0) go to 300
  190       continue
c
c           determine the direction p.
c
            call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
c
c           store the direction p and x + p. calculate the norm of p.
c
            do 200 j = 1, n
               wa1(j) = -wa1(j)
               wa2(j) = x(j) + wa1(j)
               wa3(j) = diag(j)*wa1(j)
  200          continue
            pnorm = enorm(n,wa3)
c
c           on the first iteration, adjust the initial step bound.
c
            if (iter .eq. 1) delta = dmin1(delta,pnorm)
c
c           evaluate the function at x + p and calculate its norm.
c
            iflag = 1
            call fcn(n,wa2,wa4,fjac,ldfjac,iflag)
            nfev = nfev + 1
            if (iflag .lt. 0) go to 300
            fnorm1 = enorm(n,wa4)
c
c           compute the scaled actual reduction.
c
            actred = -one
            if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
c
c           compute the scaled predicted reduction.
c
            l = 1
            do 220 i = 1, n
               sum = zero
               do 210 j = i, n
                  sum = sum + r(l)*wa1(j)
                  l = l + 1
  210             continue
               wa3(i) = qtf(i) + sum
  220          continue
            temp = enorm(n,wa3)
            prered = zero
            if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
c
c           compute the ratio of the actual to the predicted
c           reduction.
c
            ratio = zero
            if (prered .gt. zero) ratio = actred/prered
c
c           update the step bound.
c
            if (ratio .ge. p1) go to 230
               ncsuc = 0
               ncfail = ncfail + 1
               delta = p5*delta
               go to 240
  230       continue
               ncfail = 0
               ncsuc = ncsuc + 1
               if (ratio .ge. p5 .or. ncsuc .gt. 1)
     *            delta = dmax1(delta,pnorm/p5)
               if (dabs(ratio-one) .le. p1) delta = pnorm/p5
  240       continue
c
c           test for successful iteration.
c
            if (ratio .lt. p0001) go to 260
c
c           successful iteration. update x, fvec, and their norms.
c
            do 250 j = 1, n
               x(j) = wa2(j)
               wa2(j) = diag(j)*x(j)
               fvec(j) = wa4(j)
  250          continue
            xnorm = enorm(n,wa2)
            fnorm = fnorm1
            iter = iter + 1
  260       continue
c
c           determine the progress of the iteration.
c
            nslow1 = nslow1 + 1
            if (actred .ge. p001) nslow1 = 0
            if (jeval) nslow2 = nslow2 + 1
            if (actred .ge. p1) nslow2 = 0
c
c           test for convergence.
c
            if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
            if (info .ne. 0) go to 300
c
c           tests for termination and stringent tolerances.
c
            if (nfev .ge. maxfev) info = 2
            if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
            if (nslow2 .eq. 5) info = 4
            if (nslow1 .eq. 10) info = 5
            if (info .ne. 0) go to 300
c
c           criterion for recalculating jacobian.
c
            if (ncfail .eq. 2) go to 290
c
c           calculate the rank one modification to the jacobian
c           and update qtf if necessary.
c
            do 280 j = 1, n
               sum = zero
               do 270 i = 1, n
                  sum = sum + fjac(i,j)*wa4(i)
  270             continue
               wa2(j) = (sum - wa3(j))/pnorm
               wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
               if (ratio .ge. p0001) qtf(j) = sum
  280          continue
c
c           compute the qr factorization of the updated jacobian.
c
            call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
            call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
            call r1mpyq(1,n,qtf,1,wa2,wa3)
c
c           end of the inner loop.
c
            jeval = .false.
            go to 180
  290    continue
c
c        end of the outer loop.
c
         go to 30
  300 continue
c
c     termination, either normal or user imposed.
c
      if (iflag .lt. 0) info = iflag
      iflag = 0
      if (nprint .gt. 0) call fcn(n,x,fvec,fjac,ldfjac,iflag)
      return
c
c     last card of subroutine hybrj.
c
      end

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

      subroutine qform(m,n,q,ldq,wa)
      integer m,n,ldq
      double precision q(ldq,m),wa(m)
c     **********
c
c     subroutine qform
c
c     this subroutine proceeds from the computed qr factorization of
c     an m by n matrix a to accumulate the m by m orthogonal matrix
c     q from its factored form.
c
c     the subroutine statement is
c
c       subroutine qform(m,n,q,ldq,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a and the order of q.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       q is an m by m array. on input the full lower trapezoid in
c         the first min(m,n) columns of q contains the factored form.
c         on output q has been accumulated into a square matrix.
c
c       ldq is a positive integer input variable not less than m
c         which specifies the leading dimension of the array q.
c
c       wa is a work array of length m.
c
c     subprograms called
c
c       fortran-supplied ... min0
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,j,jm1,k,l,minmn,np1
      double precision one,sum,temp,zero
      data one,zero /1.0d0,0.0d0/
c
c     zero out upper triangle of q in the first min(m,n) columns.
c
      minmn = min0(m,n)
      if (minmn .lt. 2) go to 30
      do 20 j = 2, minmn
         jm1 = j - 1
         do 10 i = 1, jm1
            q(i,j) = zero
   10       continue
   20    continue
   30 continue
c
c     initialize remaining columns to those of the identity matrix.
c
      np1 = n + 1
      if (m .lt. np1) go to 60
      do 50 j = np1, m
         do 40 i = 1, m
            q(i,j) = zero
   40       continue
         q(j,j) = one
   50    continue
   60 continue
c
c     accumulate q from its factored form.
c
      do 120 l = 1, minmn
         k = minmn - l + 1
         do 70 i = k, m
            wa(i) = q(i,k)
            q(i,k) = zero
   70       continue
         q(k,k) = one
         if (wa(k) .eq. zero) go to 110
         do 100 j = k, m
            sum = zero
            do 80 i = k, m
               sum = sum + q(i,j)*wa(i)
   80          continue
            temp = sum/wa(k)
            do 90 i = k, m
               q(i,j) = q(i,j) - temp*wa(i)
   90          continue
  100       continue
  110    continue
  120    continue
      return
c
c     last card of subroutine qform.
c
      end

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

      subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
      integer m,n,lda,lipvt
      integer ipvt(lipvt)
      logical pivot
      double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
c     **********
c
c     subroutine qrfac
c
c     this subroutine uses householder transformations with column
c     pivoting (optional) to compute a qr factorization of the
c     m by n matrix a. that is, qrfac determines an orthogonal
c     matrix q, a permutation matrix p, and an upper trapezoidal
c     matrix r with diagonal elements of nonincreasing magnitude,
c     such that a*p = q*r. the householder transformation for
c     column k, k = 1,2,...,min(m,n), is of the form
c
c                           t
c           i - (1/u(k))*u*u
c
c     where u has zeros in the first k-1 positions. the form of
c     this transformation and the method of pivoting first
c     appeared in the corresponding linpack subroutine.
c
c     the subroutine statement is
c
c       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       a is an m by n array. on input a contains the matrix for
c         which the qr factorization is to be computed. on output
c         the strict upper trapezoidal part of a contains the strict
c         upper trapezoidal part of r, and the lower trapezoidal
c         part of a contains a factored form of q (the non-trivial
c         elements of the u vectors described above).
c
c       lda is a positive integer input variable not less than m
c         which specifies the leading dimension of the array a.
c
c       pivot is a logical input variable. if pivot is set true,
c         then column pivoting is enforced. if pivot is set false,
c         then no column pivoting is done.
c
c       ipvt is an integer output array of length lipvt. ipvt
c         defines the permutation matrix p such that a*p = q*r.
c         column j of p is column ipvt(j) of the identity matrix.
c         if pivot is false, ipvt is not referenced.
c
c       lipvt is a positive integer input variable. if pivot is false,
c         then lipvt may be as small as 1. if pivot is true, then
c         lipvt must be at least n.
c
c       rdiag is an output array of length n which contains the
c         diagonal elements of r.
c
c       acnorm is an output array of length n which contains the
c         norms of the corresponding columns of the input matrix a.
c         if this information is not needed, then acnorm can coincide
c         with rdiag.
c
c       wa is a work array of length n. if pivot is false, then wa
c         can coincide with rdiag.
c
c     subprograms called
c
c       minpack-supplied ... dpmpar,enorm
c
c       fortran-supplied ... dmax1,dsqrt,min0
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,j,jp1,k,kmax,minmn
      double precision ajnorm,epsmch,one,p05,sum,temp,zero
      double precision dpmpar,enorm
      data one,p05,zero /1.0d0,5.0d-2,0.0d0/
c
c     epsmch is the machine precision.
c
      epsmch = dpmpar(1)
c
c     compute the initial column norms and initialize several arrays.
c
      do 10 j = 1, n
         acnorm(j) = enorm(m,a(1,j))
         rdiag(j) = acnorm(j)
         wa(j) = rdiag(j)
         if (pivot) ipvt(j) = j
   10    continue
c
c     reduce a to r with householder transformations.
c
      minmn = min0(m,n)
      do 110 j = 1, minmn
         if (.not.pivot) go to 40
c
c        bring the column of largest norm into the pivot position.
c
         kmax = j
         do 20 k = j, n
            if (rdiag(k) .gt. rdiag(kmax)) kmax = k
   20       continue
         if (kmax .eq. j) go to 40
         do 30 i = 1, m
            temp = a(i,j)
            a(i,j) = a(i,kmax)
            a(i,kmax) = temp
   30       continue
         rdiag(kmax) = rdiag(j)
         wa(kmax) = wa(j)
         k = ipvt(j)
         ipvt(j) = ipvt(kmax)
         ipvt(kmax) = k
   40    continue
c
c        compute the householder transformation to reduce the
c        j-th column of a to a multiple of the j-th unit vector.
c
         ajnorm = enorm(m-j+1,a(j,j))
         if (ajnorm .eq. zero) go to 100
         if (a(j,j) .lt. zero) ajnorm = -ajnorm
         do 50 i = j, m
            a(i,j) = a(i,j)/ajnorm
   50       continue
         a(j,j) = a(j,j) + one
c
c        apply the transformation to the remaining columns
c        and update the norms.
c
         jp1 = j + 1
         if (n .lt. jp1) go to 100
         do 90 k = jp1, n
            sum = zero
            do 60 i = j, m
               sum = sum + a(i,j)*a(i,k)
   60          continue
            temp = sum/a(j,j)
            do 70 i = j, m
               a(i,k) = a(i,k) - temp*a(i,j)
   70          continue
            if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
            temp = a(j,k)/rdiag(k)
            rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
            if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
            rdiag(k) = enorm(m-j,a(jp1,k))
            wa(k) = rdiag(k)
   80       continue
   90       continue
  100    continue
         rdiag(j) = -ajnorm
  110    continue
      return
c
c     last card of subroutine qrfac.
c
      end

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

      subroutine r1mpyq(m,n,a,lda,v,w)
      integer m,n,lda
      double precision a(lda,n),v(n),w(n)
c     **********
c
c     subroutine r1mpyq
c
c     given an m by n matrix a, this subroutine computes a*q where
c     q is the product of 2*(n - 1) transformations
c
c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c     and gv(i), gw(i) are givens rotations in the (i,n) plane which
c     eliminate elements in the i-th and n-th planes, respectively.
c     q itself is not given, rather the information to recover the
c     gv, gw rotations is supplied.
c
c     the subroutine statement is
c
c       subroutine r1mpyq(m,n,a,lda,v,w)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       a is an m by n array. on input a must contain the matrix
c         to be postmultiplied by the orthogonal matrix q
c         described above. on output a*q has replaced a.
c
c       lda is a positive integer input variable not less than m
c         which specifies the leading dimension of the array a.
c
c       v is an input array of length n. v(i) must contain the
c         information necessary to recover the givens rotation gv(i)
c         described above.
c
c       w is an input array of length n. w(i) must contain the
c         information necessary to recover the givens rotation gw(i)
c         described above.
c
c     subroutines called
c
c       fortran-supplied ... dabs,dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,j,nmj,nm1
      double precision cos,one,sin,temp
      data one /1.0d0/
c
c     apply the first set of givens rotations to a.
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 50
      do 20 nmj = 1, nm1
         j = n - nmj
         if (dabs(v(j)) .gt. one) cos = one/v(j)
         if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
         if (dabs(v(j)) .le. one) sin = v(j)
         if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
         do 10 i = 1, m
            temp = cos*a(i,j) - sin*a(i,n)
            a(i,n) = sin*a(i,j) + cos*a(i,n)
            a(i,j) = temp
   10       continue
   20    continue
c
c     apply the second set of givens rotations to a.
c
      do 40 j = 1, nm1
         if (dabs(w(j)) .gt. one) cos = one/w(j)
         if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
         if (dabs(w(j)) .le. one) sin = w(j)
         if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
         do 30 i = 1, m
            temp = cos*a(i,j) + sin*a(i,n)
            a(i,n) = -sin*a(i,j) + cos*a(i,n)
            a(i,j) = temp
   30       continue
   40    continue
   50 continue
      return
c
c     last card of subroutine r1mpyq.
c
      end

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

      subroutine r1updt(m,n,s,ls,u,v,w,sing)
      integer m,n,ls
      logical sing
      double precision s(ls),u(m),v(n),w(m)
c     **********
c
c     subroutine r1updt
c
c     given an m by n lower trapezoidal matrix s, an m-vector u,
c     and an n-vector v, the problem is to determine an
c     orthogonal matrix q such that
c
c                   t
c           (s + u*v )*q
c
c     is again lower trapezoidal.
c
c     this subroutine determines q as the product of 2*(n - 1)
c     transformations
c
c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c     where gv(i), gw(i) are givens rotations in the (i,n) plane
c     which eliminate elements in the i-th and n-th planes,
c     respectively. q itself is not accumulated, rather the
c     information to recover the gv, gw rotations is returned.
c
c     the subroutine statement is
c
c       subroutine r1updt(m,n,s,ls,u,v,w,sing)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of s.
c
c       n is a positive integer input variable set to the number
c         of columns of s. n must not exceed m.
c
c       s is an array of length ls. on input s must contain the lower
c         trapezoidal matrix s stored by columns. on output s contains
c         the lower trapezoidal matrix produced as described above.
c
c       ls is a positive integer input variable not less than
c         (n*(2*m-n+1))/2.
c
c       u is an input array of length m which must contain the
c         vector u.
c
c       v is an array of length n. on input v must contain the vector
c         v. on output v(i) contains the information necessary to
c         recover the givens rotation gv(i) described above.
c
c       w is an output array of length m. w(i) contains information
c         necessary to recover the givens rotation gw(i) described
c         above.
c
c       sing is a logical output variable. sing is set true if any
c         of the diagonal elements of the output s are zero. otherwise
c         sing is set false.
c
c     subprograms called
c
c       minpack-supplied ... dpmpar
c
c       fortran-supplied ... dabs,dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more,
c     john l. nazareth
c
c     **********
      integer i,j,jj,l,nmj,nm1
      double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,
     *                 zero
      double precision dpmpar
      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
c
c     giant is the largest magnitude.
c
      giant = dpmpar(3)
c
c     initialize the diagonal element pointer.
c
      jj = (n*(2*m - n + 1))/2 - (m - n)
c
c     move the nontrivial part of the last column of s into w.
c
      l = jj
      do 10 i = n, m
         w(i) = s(l)
         l = l + 1
   10    continue
c
c     rotate the vector v into a multiple of the n-th unit vector
c     in such a way that a spike is introduced into w.
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 nmj = 1, nm1
         j = n - nmj
         jj = jj - (m - j + 1)
         w(j) = zero
         if (v(j) .eq. zero) go to 50
c
c        determine a givens rotation which eliminates the
c        j-th element of v.
c
         if (dabs(v(n)) .ge. dabs(v(j))) go to 20
            cotan = v(n)/v(j)
            sin = p5/dsqrt(p25+p25*cotan**2)
            cos = sin*cotan
            tau = one
            if (dabs(cos)*giant .gt. one) tau = one/cos
            go to 30
   20    continue
            tan = v(j)/v(n)
            cos = p5/dsqrt(p25+p25*tan**2)
            sin = cos*tan
            tau = sin
   30    continue
c
c        apply the transformation to v and store the information
c        necessary to recover the givens rotation.
c
         v(n) = sin*v(j) + cos*v(n)
         v(j) = tau
c
c        apply the transformation to s and extend the spike in w.
c
         l = jj
         do 40 i = j, m
            temp = cos*s(l) - sin*w(i)
            w(i) = sin*s(l) + cos*w(i)
            s(l) = temp
            l = l + 1
   40       continue
   50    continue
   60    continue
   70 continue
c
c     add the spike from the rank 1 update to w.
c
      do 80 i = 1, m
         w(i) = w(i) + v(n)*u(i)
   80    continue
c
c     eliminate the spike.
c
      sing = .false.
      if (nm1 .lt. 1) go to 140
      do 130 j = 1, nm1
         if (w(j) .eq. zero) go to 120
c
c        determine a givens rotation which eliminates the
c        j-th element of the spike.
c
         if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
            cotan = s(jj)/w(j)
            sin = p5/dsqrt(p25+p25*cotan**2)
            cos = sin*cotan
            tau = one
            if (dabs(cos)*giant .gt. one) tau = one/cos
            go to 100
   90    continue
            tan = w(j)/s(jj)
            cos = p5/dsqrt(p25+p25*tan**2)
            sin = cos*tan
            tau = sin
  100    continue
c
c        apply the transformation to s and reduce the spike in w.
c
         l = jj
         do 110 i = j, m
            temp = cos*s(l) + sin*w(i)
            w(i) = -sin*s(l) + cos*w(i)
            s(l) = temp
            l = l + 1
  110       continue
c
c        store the information necessary to recover the
c        givens rotation.
c
         w(j) = tau
  120    continue
c
c        test for zero diagonal elements in the output s.
c
         if (s(jj) .eq. zero) sing = .true.
         jj = jj + (m - j + 1)
  130    continue
  140 continue
c
c     move w back into the last column of the output s.
c
      l = jj
      do 150 i = n, m
         s(l) = w(i)
         l = l + 1
  150    continue
      if (s(jj) .eq. zero) sing = .true.
      return
c
c     last card of subroutine r1updt.
c
      end
