*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     nlmix.f
*
*     Copyright (C) 1991, 1992.
*
*     Marie Davidian and A. Ronald Gallant 
*     P.O. Box 5513, Raleigh NC 27650-5513, USA   
*
*     Permission to use, copy, modify, and distribute this
*     software and its documentation for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*
*     Program control is by entries in a parmfile that looks like this:
*
*     monoexp       name    format(a10  )
*        50         ncases  format(i5,5x)
*        10    1    npt,iq  format(2i5  ) (npt MCarlo if iq=0, npt quad if iq=1)
*         2         lx      format(i5,5x)
*         3         ltau    format(i5,5x)
*         2         M       format(i5,5x)
*         2    0    Kz,Iz   format(2i5  )
*         0    0    Kx,Ix   format(2i5  )
*       400         maxitr  format(i5,5x)
*        0.1D-05    ftol    format(d10.1)
*       956         iseed   format(i5,5x)
*        0.4D+01    sfac    format(d10.1)
*        30    1    ngrd,jt format(2i5  ) (plot points for joint if jt=1)
*         0         iunit9  format(i5,5x) (0 parmfile, 1 eff, 2 sim, 4 plt)
*         1         iCvec   format(i5,5x) (compute Cvec if iCvec=1 & iq=1)
*         1         muzero  format(i5,5x) (impose Ez=0 if muzero=1)
*                   blank line
*     Start values for tau.
*           0.100000000000000000D+00 0   tau(  1)    format(d30.18,1x,i1)
*           0.500000000000000000D+00 0   tau(  2)    format(d30.18,1x,i1)
*           0.500000000000000000D-01 0   tau(  3)    format(d30.18,1x,i1)
*     SNP parameters of theta start values.
*         0    0    1                    Kz, Iz, ntermz  format(3i5)
*         0    0    1                    Kx, Ix, ntermx  format(3i5)
*     Start values for theta.
*           0.100000000000000000D+01 1   theta(  1)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  2)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  3)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  4)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  5)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  6)  format(d30.18,1x,i1)
*           0.000000000000000000D+00 1   theta(  7)  format(d30.18,1x,i1)
*           0.030000000000000000D+00 0   theta(  8)  format(d30.18,1x,i1)
*           0.025000000000000000D+00 0   theta(  9)  format(d30.18,1x,i1)
*           0.080000000000000000D+00 0   theta( 10)  format(d30.18,1x,i1)
*
*     Note that in the example above, the SNP model to be fitted is 
*     bigger than the model for which start values are given.  The 
*     program uses the values for the smaller model to impute start
*     values for the larger model.
*
*     The 1's and 0's following the start values in column 32 are used 
*     to impose the start value as constraint: 0 means free, 1 means 
*     constrained.  In the SNP part, A(1,1) is always constrained to be
*     1.0D+00 and the regression part to be 0.0D+00.

*     The format for a start value line is (d30.18,1x,i1).

*     The program reads from the parmfile designated by the first entry
*     in a line from the file control.dat.  The second entry in a line 
*     from control.dat defines the output parmfile.  The imputed start
*     values are a uniform[-1,1] random number times the third entry, 
*     fnew.  The values of theta from the smaller model are replaced by
*     a theta(i) plus uniform[-1,1] random number times theta(i) times 
*     the fourth entry, fold.  Similarly for tau.  The fifth entry is 
*     the seed for the random number generator, nseed.  If nseed=0 then 
*     fnew is used as a fill value for the imputed start values and the 
*     values of theta from the smaller model are left unchanged.   The 
*     format for a line from control.dat is (a12,3x,a12,3x,2f10.1,i10); 
*     an example is
*
*     parmfile.dat   out_a.dat      0.0       0.0       8746
*     parmfile.dat   out_b.dat      0.0001    0.0       72365
*     parmfile.dat   out_c.dat      0.0001    0.0001    08745
*     parmfile.dat   out_d.dat      0.001     0.0       357
*     parmfile.dat   out_e.dat      0.001     0.001     8054
*
*     An output parmfile can be used as an input parmfile.  The easiest 
*     way to get a custom parmfile is to let the program generate one 
*     with the desired SNP settings from the example above and edit it.
*
*     The user must write subroutine density which gives the density in 
*     terms of the parameters tau(ltau) and random effects z(M) for each
*     case id = 1, ..., ncases.  Don't forget to compute nobs.
*
*     Dependence of the SNP density on x(lx) hasn't been debugged and 
*     probably doesn't work.  Thus, always put Kx=Ix=0 in the parmfile.
*     
*     If iunit9 is set to 4 and jt to 1, a file that can be used to plot
*     joint density is written to outfil instead of a parmfile; sfac 
*     determines the scaling and ngrd the number of grid points per 
*     coordinate.  The following values are printed to this file with 
*     format e30.18.  First, E(z), M values; Second, the upper triangle 
*     of Var(z) stored columnwise, M*(M+1)/2 values; Third, the grid 
*     increments, M values; Thereafter, z and h(z), written end to end.
*     Total file length is M+M*(M+1)/2+M+M*(2*ngrd+1)**M.  Note that z
*     is not constrained to have mean zero if muzero=0, which violates 
*     the typical regression conventions.  In reporting results, one 
*     may wish to add E(z) to the appropriate coordinates of tau and 
*     subtract E(z) from the plot coordinates before plotting them.  
*     If all the mass has been captured by the grid, a Riemann sum using 
*     zinc(i) as the weights should yield 1.  
*
*     If iunit9 is set to 4 and jt to 0, then files that can be used to 
*     plot the marginal densities are written.  The names of these files 
*     and the format, which is as above were M=1, are written to outfil
*
*     If iunit9 is set to 1, emperical Bayes estimates are written to
*     outfil instead of a parmfile, format e30.18.  This assumes that
*     density.f for isw=3 is written analagously to the example supplied
*     with the distribution.
*
*     If iunit9 is set to 2, a simulation from the SNP density is written 
*     to outfil instead of a parmfile, format e30.18.  
*
      program nlmix
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

*     Global data.

      include "global.f"

*     Local data

      real*8 work
      real*8 vec
      real*8 a,bl,bu
      real*8 c,cjac,clamda
      real*8 objf,objgrd,Rsmall,Rbig
      real*8 fn,fp,bic,aic,hq

      integer*4 isw,i,j,jj,ictr,ier
      integer*4 lwork,idummy
      integer*4 iwork
      integer*4 lvec
      integer*4 istate
      integer*4 msglvl,mode

      external  mcarlo,quad,post,mufun,confun

*     Local parameters

      real*8 bigbnd,feps,difint

      integer*4 itreff,vereff,dereff
      integer*4 debug
      integer*4 miwork,mwork,needc
      integer*4 derlev
      integer*4 nclin,ncnln,nrowa,ntotal

      parameter (debug=0)
*     parameter (debug=1)
      parameter (itreff=200)
      parameter (vereff  = -1)
      parameter (dereff  = 3)
      parameter (miwork=500,mwork=10000)
      parameter (bigbnd  = 1.d+10)
      parameter (feps    = 1.d-15)
      parameter (difint  = 1.d-8)
      parameter (derlev  = 2)
      parameter (nclin=0,nrowa=1)
      parameter (ntotal=mvec+nclin+mM)

      dimension iwork(miwork),needc(mM)
      dimension work(mwork)
      dimension vec(mvec)
      dimension a(1),bl(ntotal),bu(ntotal)
      dimension istate(ntotal)
      dimension c(mM),cjac(mM*mvec),clamda(ntotal)
      dimension objgrd(ntotal),Rsmall(mM*mM),Rbig(mvec*mvec)

      open(unit=3,file='detail.dat',status='unknown',form='formatted')
      open(unit=27,file='options.dat',status='unknown',form='formatted')
      open(unit=12,file='control.dat',status='unknown',form='formatted')

      iprint=0
      lnsize=80

*     Correspondences between SNP code and nlmix.f code.

*       SNP    nlmix
*      ------  ------
*      IDEGY    Kz
*      IDEGX    Kx
*      LINTY    Iz
*      LINTX    Ix
*      KY       ntermz
*      KX       ntermx

*     The contents of units 12 and 8 control program flow.

      ictr=0
900   read(12,12000,end=910) infile,outfil,fnew,fold,nseed
      ictr=ictr+1

*     Open infile as unit 8 and read parameters.

      call sizes(1)
      if (iprint.eq.0) then
        call spacer(2)
        call header('/nlmix version 3.0//_')
        call header('/line from control.dat//_')
        call spacer(2)
        write(3,3003)
        write(3,3004) infile,outfil,fnew,fold,nseed
        call header('/Dimensions and switches//_')
        call spacer(2)
        call sizes(3)
      end if

      L=1
      do i=1,mx
        if (M*i.lt.lx) L=i
      end do

      if (miwork.lt.2*M*L) then
        write(3,3905) 2*M*L
        stop
      end if

*     Setup SNP local data.

      call setmyx(M,L,L,Kz,Kx,Iz,Ix,ntermz,ntermx,idummy,iwork)
      call setprd(M,L,ntermz,ntermx,Kz,Kx)
      call setupm(Kz,M,ntermz)

*     Compute the length of theta and its subarrays.

      lA  = ntermz*ntermx
      lB0 = M
      lB  = M*M*L
      lR  = (M*(M+1))/2

      ltheta = lA + lB0 + lB + lR

*     Compute pointers to the stacked SNP parameter vector theta.

      ibA  = 0
      ibB0 = ibA  + lA
      ibB  = ibB0 + lB0
      ibR  = ibB  + lB

      iaA  = ibA  + 1
      iaB0 = ibB0 + 1
      iaB  = ibB  + 1
      iaR  = ibR  + 1
            
*     Error checks.
            
      if (mtheta.lt.ltheta) then
        write(3,3901) ltheta
        stop
      end if

      if (mtau.lt.ltau) then
        write(3,3902) ltau
        stop
      end if

      if (mtau.lt.M) then
        write(3,3902) M
        stop
      end if

      if (mM.lt.M) then
        write(3,3903) M
        stop
      end if

      if (mx.lt.lx+M) then
        write(3,3904) lx+M
        stop
      end if

*     Read start values from infile and close unit 8.

      call values (1)

      if (iprint.eq.0) then
        call header('/start values//_')
        call header('/tau//_')
        call dgmpnt(tau,ltau,1)
        call header('/A//_')
        call dgmpnt(theta(iaA),ntermx,ntermz)
        call header('/B0//_')
        call dgmpnt(theta(iaB0),M,1)
        call header('/B//_')
        call dgmpnt(theta(iaB),M,M*L)
        call header('/R//_')
        call dgmpnt(theta(iaR),lR,1)
       end if

*     Initial call to subroutine density with isw = 1 allows user to
*     read data from a file, setup local data, etc.

      isw=1
      call density(isw)

      if (iprint.eq.0) then
        msglvl=10
      else
        msglvl=0
      end if

*     Call to NPSOL as it appears in NPSOL documentation:
*
*     call npsol(n,nclin,ncnln,nrowa,nrowj,nrowr,
*    &             a,bl,bu,
*    &             confun,objfun,
*    &             inform,iter,istate,
*    &             c,cjac,clamda,objf,objgrd,r,x,
*    &             iw,leniw,w,lenw)
*

*     Write the options file that NPSOL uses.

      rewind 27
      write(27,27000)
      write(27,27001) bigbnd
      write(27,27002) maxitr
      write(27,27004) ftol
      write(27,27005) msglvl
      if (verify.eq.1) then
        write(27,27006) 3
      else
        write(27,27006) -1
      end if
      write(27,27007) derlev
*     write(27,27008) feps
*     write(27,27009) difint
      write(27,27900)
      rewind 27
      call npfile(27,inform)

*     Store parameter values and constraint indicators end to end 
*     in vsav and ivec repectively.

      lvsav = ltau+ltheta

      nfit1=0
      do i=1,ltau
        vsav(i)=tau(i)
        ivec(i)=itau(i)
        if (itau(i).eq.0) nfit1=nfit1+1
      end do

      nfit2=0
      do i=1,ltheta
        vsav(ltau+i)=theta(i)
        ivec(ltau+i)=itheta(i)
        if (itheta(i).eq.0) nfit2=nfit2+1
      end do

      nfit = nfit1 + nfit2
      pmax = ltau + lA - 1 + lR

*     The parameters that NPSOL is to see are stored in vec.

*     If ivec(i)=1, NPSOL doesn't get to see the ith parameter which 
*     constrains the ith parameter stay at its start value.

      j=0
      do i=1,lvsav
        if (ivec(i).eq.0) then
          j=j+1
          vec(j)=vsav(i)
        end if
      end do 
      lvec=j

      if (lvec.ne.nfit) then
        write(3,3999)
        stop
       end if

*     Fill in bl and bu.

      do i=1,lvec
        bl(i)=-bigbnd
        bu(i)= bigbnd
      end do

        ncnln=M
      if ((muzero.eq.1).and.(Kz.gt.0)) then
        do i=lvec+1,lvec+ncnln
          bl(i)=0.d0
          bu(i)=0.d0
        end do
        mode=0
        do i=1,M
          needc(i)=1
        end do
        call mufun(mode,ncnln,lvec,M,needc,vec,c,cjac,1)
      else
        ncnln=0
      end if

      net  = nfit-ncnln

      if (iprint.eq.0) then
        call header('/bl, vec, bu//_')
        call spacer(2)
        do i=1,lvec
          write(3,3001) bl(i),vec(i),bu(i)
        end do
        if (ncnln.gt.0) then
          do i=lvec+1,lvec+ncnln
            write(3,3001) bl(i),c(i-lvec),bu(i)
          end do
        end if
        call spacer(2)
      end if

*     Evaluate objective function at start values.

      mode = 0

      if (iq.eq.0) then
        call mcarlo(mode,lvec,vec,objf,objgrd,1)
      else
        call quad(mode,lvec,vec,objf,objgrd,1)
      end if

      if (iprint.eq.0) then
        call header('/Objective function at start values//_')
        call spacer(2)
        write(3,3002) objf, objf
      end if

      doCvec = 0

      if ((maxitr.ne.0).and.(iq.eq.0)) then
        call npsol(lvec,nclin,ncnln,nrowa,M,mvec,
     &             a,bl,bu,
     &             mufun,mcarlo,
     &             inform,iter,istate,
     &             c,cjac,clamda,objf,objgrd,Rbig,vec,
     &             iwork,miwork,work,mwork)
      elseif ((maxitr.ne.0).and.(iq.eq.1)) then
        call npsol(lvec,nclin,ncnln,nrowa,M,mvec,
     &             a,bl,bu,
     &             mufun,quad,
     &             inform,iter,istate,
     &             c,cjac,clamda,objf,objgrd,Rbig,vec,
     &             iwork,miwork,work,mwork)
      else
        call header('/NPSOL not called because maxitr=0//_')
        iter=0
        inform=4
      end if  

*     Copy the parameters that NPSOL has changed to vsav.

      j=0
      do i=1,lvsav
        if (ivec(i).eq.0) then
          j=j+1
          vsav(i)=vec(j)
        end if
      end do 

*     Make tau and theta agree with vsav.

      do i=1,ltau
        tau(i)=vsav(i)
      end do

      do i=1,ltheta
        theta(i)=vsav(ltau+i)
      end do

*     Compute some statistics.

      fn  = dfloat(nobs)
      fp  = dfloat(net)

      sn  = objf
      bic = sn + (fp/fn)*dlog(fn)/2.d0
      hq  = sn + (fp/fn)*dlog(dlog(fn))
      aic = sn + (fp/fn)

      lwork=(M+1)*ltheta+M+max0(2*ntermz,lR)
      if (lwork.gt.mwork) then
        write (3,3906) lwork
        stop
      endif

      do j=1,ltau
      do i=1,M
        dEz(M*(j-1)+i)=0.d0
      end do
      end do

      call mean(z,x,Ez,dEz(M*ltau+1),
     &  theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &  M,L,ntermz,ntermx,Kz,Kx,
     &  work,iwork)

      lwork=M*M*(ltheta+1)+ltheta+max0(2*ntermz,lR)+lR
      if (lwork.gt.mwork) then
        write (3,3906) lwork
        stop
      end if

      do j=1,ltau
      do i=1,lR
        dCz(lR*(j-1)+i)=0.d0
      end do
      end do

      call varcov(z,x,Cz,dCz(lR*ltau+1),
     &  theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &  M,L,ntermz,ntermx,Kz,Kx,
     &  work,iwork)

      if (lvec.ne.nfit) then
        write(3,3999)
        stop
       end if

      if ((iCvec.eq.1).and.(iq.eq.1)) then

        j=0
        do jj=1,lvsav
          if (ivec(jj).eq.0) then
            j=j+1
            do i=1,M
              dEz(M*(j-1)+i)=dEz(M*(jj-1)+i)
            end do  
            do i=1,lR
              dCz(lR*(j-1)+i)=dCz(lR*(jj-1)+i)
            end do  
          end if
        end do 

        doCvec = 1
        call quad(mode,lvec,vec,objf,objgrd,1)
        doCvec = 0

        do i=1,nfit1
        do j=1,nfit1
          Ctau(nfit1*(j-1)+i)=Cvec(nfit*(j-1)+i)
        end do
        end do

        call dsweep(Ctau,nfit1,1.d-13,ier)
        if (ier.ne.0) write(3,3908) outfil

        call dsweep(Cvec,nfit,1.d-13,ier)
        if (ier.ne.0) write(3,3908) outfil

        if (ncnln.gt.0) then
          lwork=max0(lR,M)*nfit
          if (lwork.gt.mwork) then
            write (3,3906) lwork
            stop
          endif
          call mufun(mode,ncnln,lvec,M,needc,vec,c,cjac,1)
          call dgmtra(cjac,work,M,nfit)
          call dgmaba(work,Cvec,Rsmall,nfit,M)
          call dsweep(Rsmall,M,1.d-13,ier)
          if (ier.ne.0) write(3,3908) outfil
          call dgmprd(cjac,Cvec,work,M,nfit,nfit)
          call dgmaba(work,Rsmall,Rbig,M,nfit)
          call dgmsub(Cvec,Rbig,Cvec,nfit,nfit)
          do i=1,nfit
            if (Cvec(nfit*(i-1)+i).lt.0.d0) then
              do j=1,nfit
                Cvec(nfit*(j-1)+i)=0.d0
                Cvec(nfit*(i-1)+j)=0.d0
              end do
            end if
          end do
        end if

        call dgmtra(dEz,work,M,nfit)
        call dgmaba(work,Cvec,CovEz,nfit,M)
        call dgmtra(dCz,work,lR,nfit)
        call dgmaba(work,Cvec,CovCz,nfit,lR)

      end if

      if (iprint.eq.0) then
        call header('/estimated values//_')
        call header('/tau//_')
        call dgmpnt(tau,ltau,1)
        call header('/A//_')
        call dgmpnt(theta(iaA),ntermx,ntermz)
        call header('/B0//_')
        call dgmpnt(theta(iaB0),M,1)
        call header('/B//_')
        call dgmpnt(theta(iaB),M,M*L)
        call header('/R//_')
        call dgmpnt(theta(iaR),lR,1)
        call header('/Ez//_')
        call dgmpnt(Ez,M,1)
        if ((iCvec.eq.1).and.(iq.eq.1)) then
          call header('/CovEz//_')
          call dgmpnt(CovEz,M,M)
        end if
        call header('/Cz//_')
        call dgmpnt(Cz,lR,1)
        if ((iCvec.eq.1).and.(iq.eq.1)) then
          call header('/CovCz//_')
          call dgmpnt(CovCz,lR,lR)
        end if
      end if

      if (debug.eq.1) then

        write(*,'(1x,''Lines printed to summary.dat'')')

        write(*,10001) infile,outfil
        write(*,10002) iter,maxitr,inform,Kz,Iz,net,pmax,sn,hq

        call derchk()

      end if

*     Write results to outfil.

      call sizes(2)
      call values (2)

*     If iunit9=4 overwrite the file written to unit 9 by subroutines
*     sizes and subroutine values by a file written by plot. 

      if (iunit9.eq.4) call plot()

*     If iunit9=2 overwrite the file written to unit 9 by subroutines
*     sizes and subroutine values by a file written by simlat. 

      if (iunit9.eq.2) call simlat()

*     If iunit9=1 overwrite the file written to unit 9 by subroutines
*     sizes and values by a file written by density called with isw=3. 

      if (iunit9.eq.1) then

        rewind 27
        write(27,27000)
        write(27,27001) bigbnd
        write(27,27002) itreff
        write(27,27004) ftol
        write(27,27005) msglvl
        write(27,27006) vereff
        write(27,27007) dereff
*       write(27,27008) feps
*       write(27,27009) difint
        write(27,27900)
        rewind 27
        call npfile(27,inform)

        lvec=M

        do i=1,lvec
          bl(i)=-bigbnd
          bu(i)= bigbnd
        end do

        open(unit=9,file=outfil,status='unknown',form='formatted')

        do id=1,ncases

          call pstart(vec,lvec)

          call npsol(lvec,nclin,ncnln,nrowa,1,mM,
     &               a,bl,bu,
     &               confun,post,
     &               inform,iter,istate,
     &               c,cjac,clamda,objf,objgrd,Rsmall,vec,
     &               iwork,miwork,work,mwork)

          if (inform.ne.0) write(3,3907) id

          do i=1,lvec
            z(i)=vec(i)
          end do
       
          isw=2
          call density(isw)
          isw=3
          call density(isw)

        end do

        close(9)

      end if

*     Write results to summary file.  The open and close stuff
*     is to allow the file to be read during computations (that is,
*     defeat buffering) and to increase the odds that the file is 
*     closed if the program aborts or is killed.

      open(unit=10,file='summary.dat',status='unknown',form='formatted')
      j=0
200   read(10,10000,end=210)
      j=j+1
      go to 200
210   continue
      rewind 10
      do i=1,j
        read(10,10000)
      end do
      if (ictr.eq.1) write(10,10003)
      write(10,10001) infile,outfil
      write(10,10002) iter,maxitr,inform,Kz,Iz,net,pmax,sn,hq
      close(10)

      if (iprint.eq.0) then
        call header('/summary of results//_')
        call spacer(2)
        write(3,10001) infile,outfil
        write(3,10002) iter,maxitr,inform,Kz,Iz,net,pmax,sn,hq
      end if

*     Call to subroutine density with isw = 4 allows user to
*     write data to a file to save global or local data, etc.

      isw=4
      call density(isw)

      go to 900
910   continue      

      close(3)
      close(12)
      close(27)
      stop
3001  format(1x,3d25.16)
3002  format(1x,'Starting value =',f20.8,5x,d30.18)
3003  format(6x,'infile',9x,'outfil',20x,'fnew',11x,'fold',5x,'nseed')
3004  format(6x,a12,3x,a12,3x,2f15.10,i10)
3901  format(1x,'error, main, increase mtheta to at least',i5)
3902  format(1x,'error, main, increase mtau to at least',i5)
3903  format(1x,'error, main, increase mM to at least',i5)
3904  format(1x,'error, main, increase mx to at least',i5)
3905  format(1x,'error, main, increase miwork to at least',i5)
3906  format(1x,'error, main, increase mwork to at least',i5)
3907  format(1x,'error, main, inform.ne.0 for id = ,',i5,
     &/'              parameter itreff should probably be increased.')
3908  format(1x,'error, main, outfil=',a12,', error return from dsweep')
3999  format(1x,'error, main, this should never happen')
10000 format(i5)
10001 format(1x,'infile=',a12,' outfil=',a12)
10002 format(1x,'iter=',i3,'/',i3,' inform=',i2,' Kz=',i2,' Iz=',i2,
     &  ' net=',i3,'/',i3,' sn=',f12.7,' hq=',f12.7)
10003 format(1x,'Program restart.')
12000 format(a12,3x,a12,3x,2f10.1,i10)
27000 format('Begin',74x,' ')
27001 format('   Infinite Bound            ',d15.3,35x,' ')
27002 format('   Major Iteration Limit     ',i15  ,35x,' ')
27003 format('   Function Precision        ',d15.3,35x,' ')
27004 format('   Optimality Tolerance      ',d15.3,35x,' ')
27005 format('   Major Print Level         ',i15  ,35x,' ')
27006 format('   Verify Level              ',i15  ,35x,' ')
27007 format('   Derivative Level          ',i15  ,35x,' ')
27008 format('   Function Precision        ',d15.3,35x,' ')
27009 format('   Difference Interval       ',d15.3,35x,' ')
27900 format('end  ',74x,' ')
      end


