newton_raphson_downhill Subroutine

   subroutine newton_raphson_downhill(x, fvec, eps, ndir, rel)
   implicit none
   integer(kind=I4), intent(in   )                 :: ndir
   real(kind=R8),    intent(in   )                 :: eps
   real(kind=R8),    intent(in   )                 :: rel
   real(kind=R8),    intent(inout), dimension(1:2) :: x
   real(kind=R8),    intent(out  )                 :: fvec

      real(kind=R8), dimension(1:2)    :: x0
      real(kind=R8), dimension(1:ndir) :: ct, st, ff, x1, x2

      real(kind=R8)    :: f0, fd1, fd2, ti, dfdx1, dfdx2, dx1, dx2, dt
      integer(kind=I4) :: i

      do i = 1, ndir
         ti = PI_R8*( -1. +(i -1)*2./ndir )
         ct(i) = cos(ti)
         st(i) = sin(ti)
      enddo
      x0 = x

      do
         f0  = abs_diff_sk_ku(chrom = (/x0(1)     , x0(2)     /))
         fd1 = abs_diff_sk_ku(chrom = (/x0(1) +eps, x0(2)     /))
         fd2 = abs_diff_sk_ku(chrom = (/x0(1)     , x0(2) +eps/))

         dfdx1 = (fd1 -f0)/eps
         dfdx2 = (fd2 -f0)/eps

         !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th)
         !$OMP DO PRIVATE(i, dt, dx1, dx2)
            do i = 1, ndir
               dt  = -f0/( eps*(dfdx1*ct(i) +dfdx2*st(i)) )
               dx1 = eps*ct(i)*dt
               dx2 = eps*st(i)*dt

               x1(i) = x0(1) +rel*dx1
               x2(i) = x0(2) +rel*dx2

               ff(i) = abs_diff_sk_ku(chrom = (/x1(i), x2(i)/))
            enddo
         !$OMP END DO
         !$OMP END PARALLEL

         i  = minloc(ff, 1)
         x0 = (/x1(i), x2(i)/)

         if (ff(i)<1.e-8) exit
         write(*,*) x0(1), x0(2), ff(i)
      enddo

      x    = x0
      fvec = ff(i)
   return
   endsubroutine newton_raphson_downhill