*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     quad.f
*
*     Copyright (C) 1991.
*
*     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 quad(mode,lvec,vec,objf,objgrd,nstate)
      implicit  real*8 (a-h,o-z)
      implicit  integer*4 (i-n)
      save

*     Global data

      include "global.f"

*     Local data

      real*8 objf,tmp
      real*8 vec,objgrd
      real*8 T,W,u,work
      real*8 endpts
      real*8 sum,sumdel
      real*8 fobs,factor
      real*8 h,hdel,P,PwA,Pwz,Eh,Ehdel
      real*8 R0
      real*8 weight

      integer*4 mode,lvec,nstate
      integer*4 i,j,ii,jj,isw
      integer*4 izz
      integer*4 lwork

*     Local parameters

      real*8 eps
      real*8 rootpi,root2

      integer*4 mT,mwork

      parameter (eps=1.d-300)
      parameter (rootpi=1.7724538509055160272d0)
      parameter (root2 =1.4142135623730950488d0)
      parameter (mT=130)
      parameter (mwork=500)

      dimension vec(lvec),objgrd(lvec)
      dimension sumdel(mtau+mtheta)
      dimension T(mT),W(mT),u(mM),work(mwork)
      dimension endpts(2)
      dimension izz(mM)
      dimension hdel(mtheta),PwA(mtheta),Pwz(mM),Ehdel(mtheta)
      dimension R0((mM*(mM+1))/2)

*     Error checks

      lwork=M
      if (lwork.gt.mwork) then
        write(3,3001) lwork
        stop
      end if

      lwork=2*ntermz
      if (lwork.gt.mwork) then
        write(3,3001) lwork
        stop
      end if

      lwork=npt
      if (lwork.gt.mwork) then
        write(3,3001) lwork
        stop
      end if

      if (npt.gt.mT) then
        write(3,3002) mT
        stop
      end if

*     Copy parameters NPSOL gets to see 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

*     This is an identity to pass to EHYX instead of R.

      do j=1,M
      do i=1,j
        if (i.eq.j) then
          R0(j*(j-1)/2+i)=1.d0
        else
          R0(j*(j-1)/2+i)=0.d0
        end if
      end do
      end do

      if (doCvec.eq.1) then
        do i=1,lvec*lvec
          Cvec(i)=0.d0
        end do
      end if

*     Get Gauss-Hermite quadrature points.

      call gaussq(4,npt,0.d0,0.d0,0,endpts,work,T,W)

      do i=1,npt
        T(i)=root2*T(i)
      end do

      factor=1.d0/(rootpi**M)

*     Compute log likelihood objf and its gradient objgrd.

      fobs=dfloat(nobs)

      objf=0.d0

      do i=1,lvec
        objgrd(i)=0.d0
      end do

      do id=1,ncases

*     Integrate the d density with respect to the SNP density and
*     and store it in sum. 
*     Compute derivatives with respect to the parameters tau and 
*     theta and store them in sumdel.

        sum=0.d0

        do j=1,lvsav
          sumdel(j)=0.d0
        end do

        do i=1,M
          izz(i)=1
        end do
        izz(1)=izz(1)-1

        do j=1,npt**M

          do jj=1,M
            if (izz(jj).lt.npt) then
              izz(jj)=izz(jj)+1
              go to 10
            else
              izz(jj)=1
            end if
          end do
 10       continue

          do jj=1,M
            z(jj)=0.d0
          end do

          do jj=1,M
            u(jj)=T(izz(jj))
            do ii=1,jj
              z(ii)=z(ii)+theta(ibR+jj*(jj-1)/2+ii)*u(jj)
            end do
          end do

          isw=2
          call density(isw)

          call pzx(u,x,P,PwA,Pwz,theta(iaA),
     &      M,L,ntermz,ntermx,work)

          call ehyx(u,x,Eh,Ehdel,
     &      theta(iaA),theta(iaB0),theta(iaB),R0,
     &      M,L,ntermz,ntermx,Kz,Kx,
     &      work)

          h=factor*P/Eh

          do jj=1,lA
            hdel(ibA+jj)
     &        = factor*(PwA(ibA+jj)/Eh - (P/Eh**2)*Ehdel(ibA+jj))
          end do

          weight=1.d0
          do jj=1,M
            weight=weight*W(izz(jj))
          end do

          sum=sum+h*d*weight

          do jj=1,ltau
            sumdel(jj)=sumdel(jj)+h*dwtau(jj)*weight
          end do

          do jj=1,lA
            sumdel(ltau+ibA+jj)=sumdel(ltau+ibA+jj)
     &        + hdel(ibA+jj)*d*weight
          end do

          do jj=1,M
          do ii=1,jj
            sumdel(ltau+ibR+jj*(jj-1)/2+ii)
     &        = sumdel(ltau+ibR+jj*(jj-1)/2+ii)
     &         + h*dwz(ii)*T(izz(jj))*weight
          end do
          end do

        end do

        sum=sum+eps

        objf=objf-dlog(sum)/fobs

        tmp=sum*fobs

        j=0
        do i=1,lvsav
          if (ivec(i).eq.0) then
            j=j+1
            objgrd(j)=objgrd(j)-sumdel(i)/tmp
          end if
        end do

        if (doCvec.eq.1) then

          ii=0
          do i=1,lvsav
          if (ivec(i).eq.0) ii=ii+1
          jj=0
          do j=1,lvsav
          if (ivec(j).eq.0) jj=jj+1

            if ((ivec(i).eq.0).and.(ivec(j).eq.0)) then
              Cvec(lvec*(jj-1)+ii)=Cvec(lvec*(jj-1)+ii)
     &          + (sumdel(i)/sum)*(sumdel(j)/sum)
            end if

          end do
          end do

        end if
  
      end do

      return
3001  format(1x,'error, quad, increase mwork to at least',i5)
3002  format(1x,'error, quad, npt must be less than',i5)
      end

