*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     post.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 post(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
      real*8 vec,objgrd
      real*8 u,work,sum,tmp
      real*8 h,hwz,e,ewu,P,PwA,Pwu
      real*8 R0

      integer*4 mode,lvec,nstate
      integer*4 i,j,isw
      integer*4 lwork

*     Local parameters

      real*8 eps
      integer*4 mwork

      parameter (eps=1.d-300)
      parameter (mwork=500)

      dimension vec(lvec),objgrd(lvec)
      dimension u(mM),work(mwork)
      dimension hwz(mM),PwA(mtheta),Pwu(mM),ewu(mM)
      dimension R0(mM*(mM+1)/2)

*     Error checks

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

      if (lvec.ne.M) then
        write(3,3002) 
        stop
      end if

*     Copy parameters NPSOL gets to see to z.

      do i=1,M
        z(i)=vec(i)
      end do 

*     Compute the posterior objf and its gradient objgrd.

      do i=1,lR
        R0(i)=theta(ibR+i)
      end do

      call drinv(R0,M)

      do i=1,M
        u(i)=0.d0
      end do

      do j=1,M
      do i=1,j
       u(i)=u(i)+R0(j*(j-1)/2+i)*z(j)
      end do
      end do

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

      sum=0.d0

      do i=1,M
        sum=sum+u(i)**2
      end do

      e=dexp(-sum/2.d0)

      do i=1,M
        ewu(i)=-e*u(i)
      end do

      h=P*e

      do i=1,M
        hwz(i)=0.d0
      end do

      do j=1,M
      do i=1,j
        hwz(j) = hwz(j) + (e*Pwu(i) + P*ewu(i))*R0(j*(j-1)/2+i)
      end do
      end do

      isw=2
      call density(isw)

      tmp=h*d+eps
      
      objf=-dlog(tmp)

      do i=1,M
        objgrd(i) = -(d*hwz(i) + h*dwz(i))/tmp
      end do

      return
3001  format(1x,'error, post, increase mwork to at least',i5)
3002  format(1x,'error, post, M and lvec don''t agree')
      end
