program ernewcox implicit double precision(a-h, o-z) parameter (kmax=20, nbmax=20, nmax=90000, nqmax=5) dimension x(kmax,nbmax), a(kmax,kmax) dimension adot(nqmax,kmax,kmax) dimension b(kmax,kmax), bdot(nqmax,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), cov(nbmax,nbmax) dimension biggam(nqmax,nqmax) dimension xnmat(kmax,nbmax), xdmat(nqmax,kmax,nbmax) dimension csc(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, isol, 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 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 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 call estsol(isol,beta,ifail) 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 (ifail .eq. 1) then write (3,*) 'Method failed' write(3,*) write(3,*) endif end C********************************************************************* C********************************************************************* subroutine estsol(isol,beta,ifail) implicit double precision(a-h, o-z) parameter (kmax=20, nbmax=20, nmax=90000) 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) 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 ifail = 0 do l = 1, nb beta(l) = 0.0d0 enddo write(3,*) 'Solution method:' write(3,*) write(3,*) 'NEWT: Simple Newton-Raphson iteration ', $ 'with analytic Jacobian' write(3,*) write(3,*) call newt(beta,ifail) 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 C do i = 1, n C indx(i) = i C enddo C call dsvrgp(n,t,t,indx) 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 do i = 1, n c write(4,*) indx(i), t(i) c 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) do 5 l = 1, nb cscor(l) = 0.0d0 5 continue do 15 i = 1, n phival(i) = phi(i,beta) do 8 l = 1, nb psi1v(i,l) = psi1(i,l,beta) 8 continue 15 continue do 90 m = 1, ndist if (evcnt(m) .eq. 0.0d0) goto 90 tcurr = t(ilow(m)) dden = 0.0d0 do 55 l = 1, nb dnum(l) = 0.0d0 55 continue do 65 j = ilow(m), n if (t0(j) .gt. tcurr) goto 65 dden = dden + phival(j) do 62 l = 1, nb dnum(l) = dnum(l) + psi1v(j,l) 62 continue 65 continue do 70 l = 1, nb addtrm = evcnt(m)*dnum(l)/dden cscor(l) = cscor(l) - addtrm 70 continue do 80 i = ilow(m), iup(m) if (del(i) .eq. 1.0d0) then do 75 l = 1, nb1 cscor(l) = cscor(l) + xnmat(kv(i),l) 75 continue do 76 l = 1, nb2 cscor(l+nb1) = cscor(l+nb1) + x2(i,l) 76 continue endif 80 continue 90 continue do 95 l = 1, nb cscor(l) = -cscor(l)/dn 95 continue 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) do 5 l1 = 1, nb do 5 l2 = 1, nb djac(l1,l2) = 0.0d0 5 continue do 10 i = 1, n phival(i) = phi(i,beta) do 8 l1 = 1, nb psi1v(i,l1) = psi1(i,l1,beta) do 7 l2 = l1, nb psi2v(i,l1,l2) = psi2(i,l1,l2,beta) 7 continue 8 continue 10 continue do 90 m = 1, ndist if (evcnt(m) .eq. 0.0d0) goto 90 tcurr = t(ilow(m)) dden = 0.0d0 do 55 l1 = 1, nb dnum1(l1) = 0.0d0 do 50 l2 = l1, nb dnum2(l1,l2) = 0.0d0 50 continue 55 continue do 65 j = ilow(m), n if (t0(j) .gt. tcurr) goto 65 dden = dden + phival(j) do 62 l1 = 1, nb dnum1(l1) = dnum1(l1) + psi1v(j,l1) do 60 l2 = l1, nb c dnum2(l1,l2) = dnum2(l1,l2) + psi2(j,l1,l2,beta) dnum2(l1,l2) = dnum2(l1,l2) + psi2v(j,l1,l2) 60 continue 62 continue 65 continue do 70 l1 = 1, nb do 69 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 69 continue 70 continue 90 continue do 95 l1 = 1, nb djac(l1,l1) = djac(l1,l1)/dn do 93 l2 = l1+1, nb djac(l1,l2) = djac(l1,l2)/dn djac(l2,l1) = djac(l1,l2) 93 continue 95 continue 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 10 l = 1, nb2 x2scor = x2scor + beta(nb1+l)*x2(i,l) 10 continue do 40 is = 1, k expterm1 = 0.0d0 do 20 l = 1, nb1 expterm1 = expterm1 + beta(l)*x(is,l) 20 continue expterm1 = expterm1 + x2scor expterm1 = dexp(expterm1) phi = phi + b(ir,is)*expterm1 40 continue 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 10 l = 1, nb2 x2scor = x2scor + beta(nb1+l)*x2(i,l) 10 continue do 40 is = 1, k expterm1 = 0.0d0 do 20 l = 1, nb1 expterm1 = expterm1 + beta(l)*x(is,l) 20 continue expterm1 = expterm1 + x2scor expterm1 = dexp(expterm1) phidot = phidot + bdot(iq,ir,is)*expterm1 40 continue 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 10 l0 = 1, nb2 x2scor = x2scor + beta(nb1+l0)*x2(i,l0) 10 continue do 40 is = 1, k expterm1 = 0.0d0 do 20 l1 = 1, nb1 expterm1 = expterm1 + beta(l1)*x(is,l1) 20 continue 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 40 continue 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 10 l0 = 1, nb2 x2scor = x2scor + beta(nb1+l0)*x2(i,l0) 10 continue do 40 is = 1, k expterm1 = 0.0d0 do 20 l1 = 1, nb1 expterm1 = expterm1 + beta(l1)*x(is,l1) 20 continue 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 40 continue 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 10 l0 = 1, nb2 x2scor = x2scor + beta(nb1+l0)*x2(i,l0) 10 continue do 40 is = 1, k expterm1 = 0.0d0 do 20 l1a = 1, nb1 expterm1 = expterm1 + beta(l1a)*x(is,l1a) 20 continue 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 40 continue 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 do 5 l1 = 1, nb do 5 l2 = 1, nb varmat(l1,l2) = 0.0d0 5 continue do 7 iq = 1, nq do 7 l = 1, nb udot(l,iq) = 0.0d0 7 continue do 10 i = 1, n phival(i) = phi(i,beta) do 8 l1 = 1, nb psi1v(i,l1) = psi1(i,l1,beta) 8 continue 10 continue do 25 iq = 1, nq do 20 i = 1, n pvdot(iq,i) = phidot(iq,i,beta) do 18 l1 = 1, nb ps1dot(iq,i,l1) = psi1dot(iq,i,l1,beta) 18 continue 20 continue 25 continue do 80 m = 1, ndist if (evcnt(m) .eq. 0.0d0) goto 80 tcurr = t(ilow(m)) dden = 0.0d0 do 30 iq = 1, nq ddendt(iq) = 0.0d0 30 enddo do 55 l1 = 1, nb dnum(l1) = 0.0d0 do 54 iq = 1, nq dnd(iq,l1) = 0.0d0 54 continue 55 continue do 65 j = ilow(m), n if (t0(j) .gt. tcurr) goto 65 dden = dden + phival(j) do 58 iq = 1, nq ddendt(iq) = ddendt(iq) + pvdot(iq,j) 58 continue do 62 l1 = 1, nb dnum(l1) = dnum(l1) + psi1v(j,l1) do 60 iq = 1, nq dnd(iq,l1) = dnd(iq,l1) + ps1dot(iq,j,l1) 60 continue 62 continue 65 continue do 67 iq = 1, nq do 67 l = 1, nb addtrm = - (dnd(iq,l)/dden) $ + (dnum(l)/dden)*(ddendt(iq)/dden) udot(l,iq) = udot(l,iq) + evcnt(m)*addtrm/dn 67 continue do 75 i = ilow(m), iup(m) if (del(i) .eq. 0.0d0) goto 75 do 72 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 69 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 69 continue if (l1 .le. nb1) then do 70 iq = 1, nq udot(l1,iq) = udot(l1,iq) + xdmat(iq,kv(i),l1)/dn 70 continue endif 72 continue 75 continue 80 continue do 90 l1 = 1, nb do 90 l2 = l1, nb do 88 iq1 = 1, nq do 88 iq2 = 1, nq addtrm = (dn**2)*udot(l1,iq1)*biggam(iq1,iq2)*udot(l2,iq2) varmat(l1,l2) = varmat(l1,l2) + addtrm 88 continue 90 continue do 95 l1 = 1, nb varmat(l1,l1) = varmat(l1,l1)/dn do 92 l2 = l1+1, nb varmat(l1,l2) = varmat(l1,l2)/dn varmat(l2,l1) = varmat(l1,l2) 92 continue 95 continue 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 105 l1 = 1, nb cov(l1,l1) = cov(l1,l1)/dn do 103 l2 = l1+1, nb cov(l1,l2) = cov(l1,l2)/dn cov(l2,l1) = cov(l1,l2) 103 continue 105 continue 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********************************************************************* subroutine newt(beta,ifail) implicit double precision (a-h,o-z) parameter (kmax=20, nbmax=20, nmax=90000) parameter (mxitr=100, tol=1.0d-3) external cscomp, jaccmp 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), bdiff(nbmax) dimension cscor(nbmax), djac(nbmax,nbmax) dimension djinv(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 C WE USE HERE A NAIVE IMPLEMENTATION OF NEWTON-RAPHSON C IF THERE ARE PROBLEMS WITH CONVERGENCE, A MORE SOPHISTICATED C ALGORITHM MAY BE NEEDED nb = nb1 + nb2 itr = 0 icflg = 0 do while ((itr .lt. mxitr) .and. (icflg .eq. 0)) itr = itr + 1 ipth = 1 call cscomp(beta,cscor,nb) call jaccmp(nb,beta,djac) eps = 1.0d-6 call gjr(djac,djinv,nb,eps,msing) do i = 1, nb bdiff(i) = 0.0d0 do l = 1, nb bdiff(i) = bdiff(i) + djinv(i,l)*cscor(l) enddo enddo icflg = 1 do ib = 1, nb beta(ib) = beta(ib) - bdiff(ib) rdiff = dabs(bdiff(ib))/dmax1(1.0d0,dabs(beta(ib))) if (rdiff .gt. tol) icflg=0 enddo enddo ifail = 1-icflg return end 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