C September 15, 1998
C ChE 515 illustration of solution of a system of 
C         non-linear algebraic equations.
C
C
      implicit double precision (a-h,o-z)
      parameter (neq=2) 
C
C neq = number of equations 
C
      dimension djacob(neq,neq) !Jacobian matrix 
      dimension x(neq) ! for unknown vector 
      dimension func(neq) ! function vector
      integer ipvt(neq) !Required for solver to store pivoting information.
C
C Specify initial guess
C
      x(1) = 1.d0
      x(2) = 1.d0 
C
C Convergence criterion
C
      epsi = 1.0d-8 ! abs(maximum residual) <  epsi for convergence.
      maxiter = 10  ! maximum number of Newton iterations.
C
C Typically, it is a good practice to write separate 
C subroutines for the function and the Jacobian.
C
      do k=1,maxiter
C
      write(6,*) 'ITERATION NUMBER = ', K 
c 
      call fun(x,func) ! Evaluate vector of function values.
c
      rmax = 0.d0 
      do i=1,neq
       if (dabs(func(i)).ge.rmax) then 
         rmax = dabs(func(i))
       else
       endif 
      enddo 
c 
      if (rmax.le.epsi) then 
       write (6,*) 'Solution Reached'
       write (6,*) 'Maximum Residual = ', rmax 
       do i=1,neq
        write(6,*) 'x(',i,')' 
        write(6,200) x(i) 
200     format(f16.9)  
       enddo
        write(6,*) '*************************'
       go to 99 
      else
C
       write(6,*) 'Maximum Residual = ', rmax 
       write(6,*) '---------------------------'
c 
       call jac(x,djacob) ! Evaluate the Jacobian matrix.
c
C Call dgesv to solve the system of equations
c------------------------------------------------------
       call dgesv(neq,1,djacob,neq,ipvt,func,neq,info)
c------------------------------------------------------
c
c Update variables
c 
       do i=1,neq
        x(i) = x(i)-func(i) 
       enddo
c
      endif 
      enddo 
C 
90    write(6,*) 'No convergence after',maxiter,'iterations'  
99    stop
      end
C
C Subroutine to evaluate the function
C 
      subroutine fun (x,func)
      implicit double precision (a-h,o-z)
      parameter (neq=2)
      dimension func(neq), x(neq) 
c
      func(1) = 3.d0*x(1)*x(1)-x(2)*x(2)
      func(2) = 3.d0*x(1)*x(2)*x(2)-x(1)*x(1)*x(1)-1.d0
c
      return
      end
C
C Subroutine to evaluate the Jacobian
C 
      subroutine jac (x, djacob)
      implicit double precision (a-h,o-z)
      parameter (neq=2)
      dimension djacob(neq,neq), x(neq) 
c
      djacob(1,1) = 6.d0*x(1)
      djacob(1,2) = -2.d0*x(2)
      djacob(2,1) = 3.d0*x(2)*x(2)-3.d0*x(1)*x(1)
      djacob(2,2) = 6.d0*x(1)*x(2)
c
      return
      end 
c
c     SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
*
*  -- LAPACK driver routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     March 31, 1993
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
*     ..
*
*  Purpose
*  =======
*
*  DGESV computes the solution to a real system of linear equations
*     A * X = B,
*  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
*
*  The LU decomposition with partial pivoting and row interchanges is
*  used to factor A as
*     A = P * L * U,
*  where P is a permutation matrix, L is unit lower triangular, and U is
*  upper triangular.  The factored form of A is then used to solve the
*  system of equations A * X = B.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The number of linear equations, i.e., the order of the
*          matrix A.  N >= 0.
*
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of columns
*          of the matrix B.  NRHS >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the N-by-N coefficient matrix A.
*          On exit, the factors L and U from the factorization
*          A = P*L*U; the unit diagonal elements of L are not stored.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  IPIV    (output) INTEGER array, dimension (N)
*          The pivot indices that define the permutation matrix P;
*          row i of the matrix was interchanged with row IPIV(i).
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the N-by-NRHS matrix of right hand side matrix B.
*          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
*                has been completed, but the factor U is exactly
*                singular, so the solution could not be computed.
*
*  =====================================================================