*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     values.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
      subroutine values(isw)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      real*4 ran
      
      real*8 fran
      real*8 rhoij,cii,fii,cjj,fjj
      real*8 fn,fp,bic,aic,hq

      integer*4 mseed
      integer*4 isw,i,j,ij,iprsav
      integer*4 Kz0,Iz0,nz0
      integer*4 Kx0,Ix0,nx0
      integer*4 lthet0
      integer*4 iwork,liwork
      integer*4 alpha,beta
      integer*4 ialpha,ibeta

      
*     Uncomment the following two lines, and those right before
*     the format statements, on a Sun running under SunOS 4.1.1 
*     and perhaps other Sun OS's.

*     integer*4 ieeer,ieee_flags
*     character*32 out

      integer*4 miwork

      parameter (miwork=200)

      dimension iwork(miwork)
      dimension alpha(100),beta(100)


      if (isw.eq.1) then
        read(8,8000)

        do i=1,ltau
          read(8,8001) tau(i),itau(i)
        end do

        read(8,8000)

        read(8,8002) Kz0,Iz0,nz0
        read(8,8002) Kx0,Ix0,nx0

        read(8,8000)

        lthet0 = nz0*nx0 + M + M*M*L + lR

        do i=1,lthet0
          read(8,8001) theta(i),itheta(i)
          vsav(i)=dble(itheta(i))
        end do

        liwork = 2*M*L + M*nz0 + M*L*nx0 + nz0 + nx0

        if (liwork.gt.miwork) then
          write(3,3001) liwork
          stop
        end if

        iprsav=iprint
        iprint=1
        
        call snprst
     &   (vsav,M,0,L,Kz0,Kx0,Iz0,Ix0,
     &   0,L,Kz,Kx,Iz,Ix,0.d0,ij,iwork)

        call snprst
     &   (theta,M,0,L,Kz0,Kx0,Iz0,Ix0,
     &   0,L,Kz,Kx,Iz,Ix,fnew,ij,iwork)

        if (lthet0.ne.ltheta) then
          do i=1,ltheta
            itheta(i)=int(vsav(i))
          end do
        end if

        do i=1,lthet0
          vsav(i)=1.d0
        end do

        call snprst
     &   (vsav,M,0,L,Kz0,Kx0,Iz0,Ix0,
     &   0,L,Kz,Kx,Iz,Ix,0.d0,ij,iwork)

        iprint=iprsav

        mseed=nseed
        
        if (mseed.ne.0) then
          do i=1,ltau
            if (itau(i).ne.1) then 
              fran = -1.d0 + 2.d0*dble(ran(mseed))
              tau(i)=tau(i)+tau(i)*fold*fran
            end if
          end do
          do i=1,ltheta
            if (itheta(i).ne.1) then 
              fran = -1.d0 + 2.d0*dble(ran(mseed))
              if (vsav(i).eq.1.d0) then
                theta(i)=theta(i)+theta(i)*fold*fran
              else
                theta(i)=fnew*fran
              end if
            end if
          end do
        end if       

        
        liwork = 2*M*L + M*nz0 + M*L*nx0 + nz0 + nx0
*     The first element of the polynomial part of the SNP model
*     is always fixed at one.

        theta(ibA+1)=1.d0
        itheta(ibA+1)=1

*     The regression part of the SNP model is not used, zeroing
*     them out accomplishes this.

        do i=ibB0+1,ibB0+M
          theta(i)=0.d0
          itheta(i)=1
        end do
        do i=ibB+1,ibB+M*M*L
          theta(i)=0.d0
          itheta(i)=1
        end do

*    With these built in restrictions, pmax is the maximum number 
*    of free parameters.

        close(8)
        return
      end if

      if (isw.eq.2) then
        write(9,9000)

        do i=1,ltau
          write(9,9001) tau(i),itau(i),i
        end do

        write(9,9002)

        write(9,9003) Kz,Iz,ntermz
        write(9,9004) Kx,Ix,ntermx

        write(9,9005)

        do i=1,ltheta
          write(9,9009) theta(i),itheta(i),i
        end do

        write(9,9006)

        do i=1,M
          write(9,9018) Ez(i),i
        end do

        write(9,9007)

        do j=1,M
        do i=1,j
          write(9,9019) Cz(j*(j-1)/2+i),i,j
        end do
        end do

        write(9,9017)

        do j=1,M
          cjj = Cz(j*(j-1)/2+j)
          if (cjj.gt.0) then
            fjj = 1.d0/dsqrt(cjj)
          else
            fjj = 0.d0
          end if
        do i=1,j
          cii = Cz(i*(i-1)/2+i)
          if (cii.gt.0) then
            fii = 1.d0/dsqrt(cii)
          else
            fii = 0.d0
          end if

          rhoij = Cz(j*(j-1)/2+i)*fjj*fii

          write(9,9020) rhoij,i,j
        end do
        end do

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

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

        write(9,9101)
        write(9,9102) infile
        write(9,9103) outfil
        write(9,9104) fnew,fnew
        write(9,9114) fold,fold
        write(9,9115) nseed
        write(9,9105) iter,maxitr
        write(9,9106) inform
        write(9,9109) nobs
        write(9,9110) net,pmax
        write(9,9107) sn,sn
        write(9,9108) bic,bic
        write(9,9113) hq,hq 
        write(9,9112) aic,aic
        write(9,9111) (2*nobs)*sn,(2*nobs)*sn

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

          write(9,9010)
          j=0
          do i=1,ltau
            if (itau(i).eq.0) then
              j=j+1
              write(9,9013) tau(i),dsqrt(Cvec(nfit*(j-1)+j)),i
            else
              write(9,9015) tau(i),i
            end if
          end do
          
          write(9,9011)

          do ialpha=1,ntermz
          do ibeta=1,ntermx
            ij=ntermx*(ialpha-1)+ibeta
            call multiy(ialpha,alpha)
            call multix(ibeta,beta)
            if (itheta(ibA+ij).eq.0) then
              j=j+1
              write(9,9012) 
     &          theta(ij),dsqrt(Cvec(nfit*(j-1)+j)),ij,(alpha(i),i=1,M)
            else
              write(9,9014) theta(ij),ij,(alpha(i),i=1,M)
            end if
          end do
          end do

          do i=ibB0+1,ibB0+M
            if (itheta(i).eq.0) then
              j=j+1
              write(9,9012) 
     &          theta(i),dsqrt(Cvec(nfit*(j-1)+j)),i
            else
              write(9,9014) theta(i),i
            end if
          end do

          do i=ibB+1,ibB+M*M*L
            if (itheta(i).eq.0) then
              j=j+1
              write(9,9012) 
     &          theta(i),dsqrt(Cvec(nfit*(j-1)+j)),i
            else
              write(9,9014) theta(i),i
            end if
          end do

          do i=ibR+1,ibR+lR
            if (itheta(i).eq.0) then
              j=j+1
              write(9,9012) 
     &          theta(i),dsqrt(Cvec(nfit*(j-1)+j)),i
            else
              write(9,9014) theta(i),i
            end if
          end do

          write(9,9016)
          j=0
          do i=1,ltau
            if (itau(i).eq.0) then
              j=j+1
              write(9,9013) tau(i),dsqrt(Ctau(nfit1*(j-1)+j)),i
            else
              write(9,9015) tau(i),i
            end if
          end do

          write(9,9021)

          if (muzero.eq.1) then 
            do i=1,M
              write(9,9025) Ez(i),i
            end do
          else
            do i=1,M
              write(9,9023) Ez(i),dsqrt(CovEz(M*(i-1)+i)),i
            end do
          end if

          write(9,9022)

          do j=1,M
          do i=1,j
            ij=j*(j-1)/2+i
            write(9,9024) Cz(ij),sqrt(CovCz(lR*(ij-1)+ij)),i,j
          end do
          end do

        end if

*       Uncomment the following seven lines, and those amongst the 
*       variable declarations, on a Sun running under SunOS 4.1.1 
*       and perhaps other Sun OS's.

*       ieeer = ieee_flags('get','exception','invalid',out)
*       if (out.eq.'invalid') write(9,9901) out
*       ieeer = ieee_flags('get','exception','overflow',out)
*       if (out.eq.'overflow') write(9,9901) out
*       ieeer = ieee_flags('get','exception','division',out)
*       if (out.eq.'division') write(9,9901) out
*       ieeer = ieee_flags('clear','exception','common',out)

        close(9)
        return
      end if

3001  format(1x,'error, values, increase miwork to at least', i5)
8000  format(a10)
8001  format(d30.18,1x,i1)
8002  format(3i5)
9000  format('Start values for tau.')
9001  format(d30.18,1x,i1,3x,'tau(',i3,')    format(d30.18,1x,i1)')
9002  format('SNP parameters of theta start values.')
9003  format(3i5,20x,'Kz, Iz, ntermz  format(3i5)')
9004  format(3i5,20x,'Kx, Ix, ntermx  format(3i5)')
9005  format('Start values for theta.')
9006  format('Mean of z, Ez.')
9007  format('Variance matrix of z, Cz.')
9009  format(d30.18,1x,i1,3x,'theta(',i3,')  format(d30.18,1x,i1)')
9010  format('Values of tau (repeated) and standard errors')
9011  format
     &  ('Values of theta (repeated), standard errors,',
     &  ' and multi-indexes')
9012  format(f20.8,3x,f20.8,3x,'theta(',i3,')',3x,100i1)
9013  format(f20.8,3x,f20.8,3x,'tau(',i3,')')
9014  format(f20.8,12x,'constrained',3x,'theta(',i3,')',3x,100i1)
9015  format(f20.8,12x,'constrained',3x,'tau(',i3,')')
9016  format('Values of tau (repeated) and conditional standard errors')
9017  format
     &('Correlation matrix of z, rho.')
9018  format(d30.18,5x,'Ez(',i3,')')
9019  format(d30.18,5x,'Cz(',i3,',',i3,')')
9020  format(f23.15,12x,'rho(',i3,',',i3,')')
9021  format('Mean of z, Ez, (repeated) and standard errors.')
9022  format
     &('Variance matrix of z, Cz, (repeated) and standard errors.')
9023  format(f20.8,3x,f20.8,3x,'Ez(',i3,')')
9024  format(f20.8,3x,f20.8,3x,'Cz(',i3,',',i3,')')
9025  format(f20.8,12x,'constrained',3x,'Ez(',i3,')')
9101  format('Results from fit.')
9102  format('  infile            = ',2x,a12)
9103  format('  outfil            = ',2x,a12)
9104  format('  fnew              = ',f15.8,5x,d27.18)
9105  format('  iter/maxitr       = ',i5,' /',i5)
9106  format('  inform            = ',i5)
9107  format('  sn = objf         = ',f15.8,5x,d27.18)
9108  format('  Schwarz = bic     = ',f15.8,5x,d27.18)
9109  format('  nobs              = ',i8)
9110  format('  net/pmax          = ',i5,' /',i5)
9111  format('  -2 ln likelihood  = ',f15.8,5x,d27.18)
9112  format('  Akaike = aic      = ',f15.8,5x,d27.18)
9113  format('  Hannan&Quinn = hq = ',f15.8,5x,d27.18)
9114  format('  fold              = ',f15.8,5x,d27.18)
9115  format('  nseed             = ',i10)
9901  format('  IEEE floating-point arithmetic exception, ',a32)
      end


