! Steepest Descent Methods with a Hessian Modification from Pounds ! This program solves 10.4 5c in BandF 10th ed. program nsystems integer, parameter :: N = 3 real (kind=8), dimension(N) :: F, X, Xold, Y real (kind=8), dimension(N,N) :: H real (kind=8), parameter :: TOL = 1.0E-12; real (kind=8) :: tolcheck ! Needed for Lapack integer, dimension(N) :: IPIV integer :: INFO ! Set the initial solution vector Xold(1) = 0.0 Xold(2) = 0.0 Xold(3) = 0.0 do i=1, 10 do k=1,N X(k) = Xold(k) enddo ! Build F vector call buildf( F, X) ! Build Hessian call buildh(H, X, N ) ! Solve Linear System using LAPACK call dgesv(N,1,H,N,IPIV,F,N,INFO) ! Solve for New X do k=1,N X(k) = X(k) - F(k) enddo !Check for Tolerance tolcheck = 0.0D0 do k=1,N tolcheck = tolcheck + abs(X(k) - Xold(k)) enddo print *, X(1), X(2), X(3), tolcheck if ( tolcheck > TOL ) then do k=1,N Xold(k) = X(k) enddo endif enddo end program nsystems subroutine buildf ( F, X ) real (kind=8), dimension(*) :: F, X real (kind=8) :: x1, x2, x3 x1 = X(1) x2 = X(2) x3 = X(3) ! F(1,2,3) are the partial derivatives of the squared g function given in the ! book with respect to x1, x2, and x3. ! I used Mathematica to generate the following Fortran code F(1) = 4*(1 + x1 - x2)*(2 + x1**2 - 2*x1*(-1 + x2) - & 2.5*x2 + 2*x2**2 - x3 + x3**2) F(2) = 2*(-2.5 - 2*x1 + 4*x2)* & (2 + x1**2 - 2*x1*(-1 + x2) - 2.5*x2 + 2*x2**2 - & x3 + x3**2) F(3) = 2*(-1 + 2*x3)*(2 + x1**2 - 2*x1*(-1 + x2) - 2.5*x2 + & 2*x2**2 - x3 + x3**2) end subroutine buildf subroutine buildh ( H, X, N ) integer :: N real (kind=8), dimension(*) :: X real (kind=8), dimension(N,*) :: H real (kind=8) :: x1, x2, x3 x1 = X(1) x2 = X(2) x3 = X(3) ! Using the square of the g function as above, take the partial derivatives ! to build the Hessian matrix. Save time by remembering the rules related ! to mixed partial derivatives -- i.e. -- H_xy = H_yx. ! I used Mathematica to generate the following Fortran code H(1,1) = 12.*(1.3333333333333333 + x1**2 + x1*(2. - 2.*x2) - & 2.1666666666666665*x2 + 1.3333333333333333*x2**2 - & 0.3333333333333333*x3 + 0.3333333333333333*x3**2) H(2,2) = 16.*(1.78125 + x1**2 + x1*(2.25 - 3.*x2) - 3.75*x2 + & 3.*x2**2 - 0.5*x3 + 0.5*x3**2) H(3,3) = 12.*(0.8333333333333334 + 0.3333333333333333*x1**2 + & x1*(0.6666666666666666 - 0.6666666666666666*x2) - & 0.8333333333333334*x2 + 0.6666666666666666*x2**2 - & 1.*x3 + x3**2) H(1,2) = -12.*(1.5 + x1**2 + & x1*(2.1666666666666665 - 2.6666666666666665*x2) - & 3.*x2 + 2.*x2**2 - 0.3333333333333333*x3 + & 0.3333333333333333*x3**2) H(2,1) = H(1,2) H(1,3) = 4*(1 + x1 - x2)*(-1 + 2*x3) H(3,1) = H(1,3) H(2,3) = -4.*(1.25 + x1 - 2.*x2)*(-1. + 2.*x3) H(3,2) = H(2,3) end subroutine buildh