      implicit none
      integer n, nrhs, lwork, i, j
      parameter (n = 2, nrhs = 1, lwork = n*n)
      character uplo
      integer info, ipiv(n)
      double precision x(n), fx, gx(n), hx(n,n)
      double precision xdiff, xtol, work(lwork)
      data x / -1.2, 1.0 /, uplo / 'u' /, xdiff / 1.0 /, xtol / 10e-8 /
c
      j = 0
      do while (xdiff .gt. xtol)
         call derivfn(x, n, fx, gx, hx)
c                     Ratkaistaan Newton-yhtalo Lapackilla
         call dsysv(uplo, n, nrhs, hx, n, ipiv, gx, n, work, lwork, info)
         if (info .ne. 0) then
            if (info .gt. 0) print *, 'Hessen matriisi singulaarinen'
            stop
         end if
         xdiff = 0
         do i = 1, n
            x(i) = x(i) - gx(i)
            xdiff = xdiff + gx(i)**2
         end do
         xdiff = sqrt(xdiff)
         j = j + 1
      end do
c
      write(6, 999) (x(i), i=1,n), fx, j
 999  format('x = (', 2(2x,f8.6), ' ), fx = ', e12.6, ', j = ', i2)
      end
c                     Minimoitava kohdefunktio derivaattoineen
      subroutine derivfn(x, n, fx, gx, hx)
      implicit none
      integer n
      double precision x(n), fx, gx(n), hx(n,*)
c                     Funktion arvo
      fx = 100*(x(2) - x(1)**2)**2 + (1 - x(1))**2
c                     Gradienttivektori
      gx(1) = - 400*x(1)*(x(2) - x(1)**2) + 2*(x(1) - 1)
      gx(2) = 200*(x(2) - x(1)**2)
c                     Hessen matriisin ylakolmio
      hx(1,1) = 1200*x(1)**2 - 400*x(2) + 2
      hx(1,2) = - 400*x(1)
      hx(2,2) = 200
c
      end
