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