fdjac_mn_sub Subroutine

subroutine fdjac_mn_sub(funcv, x, m, fjac, epsfcn)

Arguments

Type IntentOptional Attributes Name
subroutine funcv(x, m, y)
Arguments
Type IntentOptional Attributes Name
real(kind=8), intent(in), dimension(:) :: x
integer :: m
real(kind=8), dimension(m) :: y
real(kind=8), intent(in) :: x(:)
integer :: m
real(kind=8) :: fjac(m,size(x))
real(kind=8), optional :: epsfcn

Source Code

subroutine fdjac_mn_sub(funcv,x,m,fjac,epsfcn)
  implicit none
  interface 
     subroutine funcv(x,m,y)
       implicit none
       integer                         :: m
       real(8),dimension(:),intent(in) :: x
       real(8),dimension(m)            :: y
     end subroutine funcv
  end interface
  integer          ::  n
  integer          ::  m
  real(8),intent(in) ::  x(:)
  real(8)            ::  x_(size(x))
  real(8)          ::  fvec(m)
  real(8)          ::  fjac(m,size(x))
  real(8),optional ::  epsfcn
  real(8)          ::  eps,eps_
  real(8)          ::  epsmch
  real(8)          ::  h,temp
  real(8)          ::  wa1(m)
  real(8)          ::  wa2(m)
  integer          :: i,j,k
  n=size(x)
  x_ = x
  eps_= 0.d0; if(present(epsfcn))eps_=epsfcn
  epsmch = epsilon(epsmch)
  eps    = sqrt(max(eps_,epsmch))
  call funcv(x_,m,fvec)
  do j=1,n
     temp = x_(j)
     h    = eps*abs(temp)
     if(h==0.d0) h = eps
     x_(j) = temp + h
     call funcv(x_,m,wa1)
     x_(j) = temp
     fjac(1:m,j) = (wa1(1:m) - fvec(1:m))/h
  enddo
end subroutine fdjac_mn_sub