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

*     Global data.

      include "global.f"

*     Local data

      real*8 vec
      real*8 objf,objgrd,c
      real*8 objsav,vecsav,del,delsav

      integer*4 i,ii,j
      integer*4 lvec
      integer*4 mode,nstate,needc

*     Local parameters

      real*8 delta

      parameter (delta=1.d-8)

      dimension vec(mvec)
      dimension objgrd(mM*mvec),c(mM)
      dimension delsav(mvec)
      dimension needc(mM)

      mode=2
      nstate=1

      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

      write(*,'(1x,''Begin objfun derivative checks'')')

      if (iq.eq.0) then
        write(*,'(1x,''Computed by mcarlo.f'')')
        call mcarlo(mode,lvec,vec,objf,objgrd,nstate)
      else
      write(*,'(1x,''Computed by quad.f'')')
        call quad(mode,lvec,vec,objf,objgrd,nstate)
      end if

      write(*,'(1x,''objf'',d25.16)') objf

      objsav=objf

      do i=1,lvec
        delsav(i)=objgrd(i)
      end do

      do i=1,lvec
        vecsav=vec(i)
        vec(i)=vec(i)+delta
        if (iq.eq.0) then
          call mcarlo(mode,lvec,vec,objf,objgrd,nstate)
        else
          call quad(mode,lvec,vec,objf,objgrd,nstate)
        end if
        vec(i)=vecsav
        del=(objf-objsav)/delta
        write(*,'(1x,''i, objgrd, del'',i5,2d25.16)') i,delsav(i),del
      end do

      write(*,'(1x,''Begin mufun.f derivative checks'')')

      do ii=1,M

        do i=1,M
          needc(i)=0
        end do
        needc(ii)=1

        call mufun(mode,M,lvec,M,needc,vec,c,objgrd,nstate)

*       call dgmpnt(objgrd,M,lvec)

        objf=c(ii)
        do j=1,lvec
          objgrd(j)=objgrd(M*(j-1)+ii)
        end do 

        write(*,'(1x,''Eyl('',i2,'') = '',d25.16)') ii, objf

        objsav=objf

        do i=1,lvec
          delsav(i)=objgrd(i)
        end do

        do i=1,lvec
          vecsav=vec(i)
          vec(i)=vec(i)+delta

          call mufun(mode,M,lvec,M,needc,vec,c,objgrd,nstate)
          objf=c(ii)
          do j=1,lvec
            objgrd(j)=objgrd(M*(j-1)+i)
          end do 

          vec(i)=vecsav
          del=(objf-objsav)/delta
          write(*,'(1x,''i, objgrd, del'',i5,2d25.16)') i,delsav(i),del
        end do

      end do

      write(*,'(1x,''Begin post.f derivative checks'')')

      lvec=M
      id=1

      do i=1,lvec
        vec(i)=Ez(i)+dsqrt(Cz(i*(i+1)/2))
      end do

      call post(mode,lvec,vec,objf,objgrd,nstate)

      write(*,'(1x,''objf'',d25.16)') objf

      objsav=objf

      do i=1,lvec
        delsav(i)=objgrd(i)
      end do

      do i=1,lvec
        vecsav=vec(i)
        vec(i)=vec(i)+delta*dsqrt(Cz(i*(i+1)/2))
        call post(mode,lvec,vec,objf,objgrd,nstate)
        vec(i)=vecsav
        del=(objf-objsav)/(delta*dsqrt(Cz(i*(i+1)/2)))
        write(*,'(1x,''i, objgrd, del'',i5,2d25.16)') i,delsav(i),del
      end do

      write(*,'(1x,''End derivative checks'')')

      return
3999  format(1x,'error, derchk, this should never happen')
      end
