Using MINPACK for curve fitting: implementation?

Question

I need to implement a non-linear fitting algorithm in Fortran and chose to use MINPACK's flavor of the Levenberg-Marquardt algorithm as a basis for the least-squares stuff. However, I seem to misunderstand something very basic about how these routines are to be used. My understanding is that I have to call the routine lmdif1 and define a function fcn to be called by the MINPACK routines that calculates the square of the difference between the data to be fitted and the model function; lmdif1 would then return the best solution for the free fitting parameters in an array x.

I have written a small test program (given below) to try this out. The program produces a set of "data" by calculating values from the model function (exp(ax) in this case) with a given parameter a (=-0.5 in this case) and changing them by a small random "error". This dataset is then fitted with the same model function, with a being the free parameter. But it doesn't really work, the fits are usually very bad, unless I almost give the "correct" solution as a starting guess.

Considering that Levenberg-Marquardt is said to be quite robust and seeing how well it works e.g. in Gnuplot, I wonder what I'm doing wrong here. Any ideas? The code below should work when linking to the MINPACK routines.

Tom

! test lmdif
program testanf
implicit none
integer, parameter :: m=10,n=1,lwa=m*n+5*n+m
integer :: info,iwa(n)
double precision :: x(n),fvec(m),wa(lwa),tol=1d-3
external :: fcn
x(1)=-0.5d0
call lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
call translate_info(info,n,tol)
write(*,*) 'a=',x(1)
end program testanf

subroutine fcn(m,n,x,fvec,iflag)
implicit none
! m: number of data points
! n: number of fitting parameters
integer, intent(in) :: m,n,iflag
integer :: i
real, allocatable, save :: t(:),fdata(:)
real :: tmin,tmax,dt,rand
double precision, intent(inout) :: x(n)
double precision, intent(out) :: fvec(m)
character(len=1) :: c
logical, save :: init=.true.
if (init) then
!  create "measured data"
   tmin=0.
   tmax=10.
   dt=(tmax-tmin)/m
   call random_seed(put=(/1,2,3,4,5,6,7,8/))
   allocate(fdata(m),t(m))
   do i=1,m
      t(i)=i*dt
      fdata(i)=fmodel(-0.5d0,t(i))
   end do
   tmax=0.1*maxval(fdata)
   do i=1,m
      call random_number(rand)
      fdata(i)=fdata(i)+(rand-0.5)*tmax ! add measurement errors
      write(10,*) t(i),fdata(i)
   end do
   init=.false.
end if
! minimize the residuals
do i=1,m
   fvec(i)=(fdata(i)-fmodel(x(1),t(i)))**2
end do
if (iflag < 0) deallocate(t,fdata)
return
contains
! model function
  real function fmodel(a,t)
  implicit none
  real, intent(in) :: t
  double precision, intent(in) :: a
  fmodel=exp(a*t)
  return
  end function fmodel
end subroutine fcn