program main

!*****************************************************************************80
!
!! MAIN is the main program for P06_OPT01.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    10 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Raman Mehra, William Kessel, James Carroll,
!    Global stability and contral analysis of aircraft at high angles of attack,
!    Technical Report CR-215-248-1, -2, -3,
!    Office of Naval Research, June 1977.
!
!    Rami Melhem, Werner Rheinboldt,
!    A Comparison of Methods for Determining Turning Points of Nonlinear Equations,
!    Computing,
!    Volume 29, Number 3, September 1982, pages 201-226.
!
!    Albert Schy, Margery Hannah,
!    Prediction of Jump Phenomena in Roll-coupled Maneuvers of Airplanes,
!    Journal of Aircraft,
!    Volume 14, Number 4, 1977,  pages 375-382.
!
!    John Young, Albert Schy, Katherine Johnson,,
!    Prediction of Jump Phenomena in Aircraft Maneuvers, Including
!    Nonlinear Aerodynamic Effects,
!    Journal of Guidance and Control,
!    Volume 1, Number 1, 1978, pages 26-31.
!
!  Local parameters:
!
!    Input, integer OPTION, the option index.
!
  implicit none

  integer ( kind = 4 ) nvar

  character ( len = 80 ) file_t_name
  character ( len = 80 ) file_x_name
  integer   ( kind = 4 ) file_t_unit
  integer   ( kind = 4 ) file_x_unit
  real    ( kind = 8 ) h
  integer ( kind = 4 ) h_reduction
  real    ( kind = 8 ) hmax
  real    ( kind = 8 ) hmin
  real    ( kind = 8 ) ht
  integer ( kind = 8 ) i
  integer ( kind = 4 ) option
  integer ( kind = 4 ) par_index
  real    ( kind = 8 ) r8_sign
  integer ( kind = 4 ) status
  integer ( kind = 4 ) step
  integer ( kind = 4 ) step_max
  real    ( kind = 8 ), allocatable, dimension ( : ) :: tan
  real    ( kind = 8 ), allocatable, dimension ( : ) :: x
  real    ( kind = 8 ), allocatable, dimension ( : ) :: x1

  option = 1

  write ( file_x_name, '(a,i1,a)' ) 'p06_opt0', option, '_x000.txt'
  write ( file_t_name, '(a,i1,a)' ) 'p06_opt0', option, '_t000.txt'

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'P06_OPT01'
  write ( *, '(a)' ) '  Compute a series of solutions for problem 6.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  The option chosen is ', option
  write ( *, '(a)' ) ' '
!
!  Get problem size.
!
  call p06_nvar ( option, nvar )

  allocate ( tan(1:nvar) )
  allocate ( x(1:nvar) )
  allocate ( x1(1:nvar) )
!
!  Get starting point.
!
  call p06_start ( option, nvar, x )
!
!  Get the tangent vector.
!
  call p06_tan ( option, nvar, x, tan )
!
!  For correction of initial point, use variable index 7.
!
  par_index = 7
!
!  Force F(X) = 0.
!
  step = -1
  write ( *, '(2x,i2,8(1x,f8.5))' ) step, x(1:nvar)

  call p06_newton ( option, nvar, x, par_index, status )

  if ( status < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Newton iteration failed on starting point.'
    return
  end if

  call p06_tan ( option, nvar, x, tan )

  step = 0
  write ( *, '(2x,i2,8(1x,f8.5))' ) step, x(1:nvar)

  call get_unit ( file_x_unit )
  open ( unit = file_x_unit, file = file_x_name, status = 'replace' )
  do i = 1, nvar
    write ( file_x_unit, '(g16.8)' ) x(i)
  end do
  close ( unit = file_x_unit )
  call file_name_inc ( file_x_name )

  call get_unit ( file_t_unit )
  open ( unit = file_t_unit, file = file_t_name, status = 'replace' )
  do i = 1, nvar
    write ( file_t_unit, '(g16.8)' ) tan(i)
  end do
  close ( unit = file_t_unit )
  call file_name_inc ( file_t_name )
!
!  Get stepsize.
!
  call p06_stepsize ( option, h, hmin, hmax )
!
!  LOOP:
!
  step_max = 30

  do step = 1, step_max
!
!  Compute the tangent vector.
!
    call p06_tan ( option, nvar, x, tan )

    h_reduction = 0
!
!  Use X + H * TAN as a starting estimate for Newton iteration.
!
    do

      if ( hmax < abs ( h ) ) then
        h = hmax * r8_sign ( h )
      end if

      if ( abs ( h ) < hmin ) then
        h = hmin * r8_sign ( h )
      end if

      x1(1:nvar) = x(1:nvar) + h * tan(1:nvar)

      par_index = 0
      call p06_newton ( option, nvar, x1, par_index, status )
!
!  If we didn't get it, can we try again?
!
      if ( status < 0 ) then

        if ( abs ( h ) <= hmin ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'P06_OPT01:'
          write ( *, '(a)' ) '  Cannot decrease stepsize any more!'
          write ( *, '(a)' ) '  Did not reach target point.'
          stop
        else
          h = h / 4.0D+00
          h_reduction = h_reduction + 1
        end if
!
!  We computed the point.
!  Should we change the stepsize?
!
      else

        if ( h_reduction == 0 ) then

          if ( status <= 1 ) then
            h = h * 4.0D+00
          else if ( status <= 3 ) then
            h = h * 2.0D+00
          else if ( 12 <= status ) then
            h = h / 4.0D+00
          else if ( 8 <= status ) then
            h = h / 2.0D+00
          end if

        end if

        exit

      end if

    end do

    call p06_tan ( option, nvar, x1, tan )
    write ( *, '(2x,i2,8(1x,f8.5))' ) step, x1(1:nvar)

    call get_unit ( file_x_unit )
    open ( unit = file_x_unit, file = file_x_name, status = 'replace' )
    do i = 1, nvar
      write ( file_x_unit, '(g16.8)' ) x1(i)
    end do
    close ( unit = file_x_unit )
    call file_name_inc ( file_x_name )

    call get_unit ( file_t_unit )
    open ( unit = file_t_unit, file = file_t_name, status = 'replace' )
    do i = 1, nvar
      write ( file_t_unit, '(g16.8)' ) tan(i)
    end do
    close ( unit = file_t_unit )
    call file_name_inc ( file_t_name )

    if ( step == step_max ) then
      exit
    end if

    x(1:nvar) = x1(1:nvar)

  end do

  deallocate ( tan )
  deallocate ( x )
  deallocate ( x1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'P06_OPT01:'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine file_name_inc ( file_name )

!*****************************************************************************80
!
!! FILE_NAME_INC increments a partially numeric filename.
!
!  Discussion:
!
!    It is assumed that the digits in the name, whether scattered or
!    connected, represent a number that is to be increased by 1 on
!    each call.  If this number is all 9's on input, the output number
!    is all 0's.  Non-numeric letters of the name are unaffected.
!
!    If the name is empty, then the routine stops.
!
!    If the name contains no digits, the empty string is returned.
!
!  Example:
!
!      Input            Output
!      -----            ------
!      'a7to11.txt'     'a7to12.txt'
!      'a7to99.txt'     'a8to00.txt'
!      'a9to99.txt'     'a0to00.txt'
!      'cat.txt'        ' '
!      ' '              STOP!
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    14 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character ( len = * ) FILE_NAME.
!    On input, a character string to be incremented.
!    On output, the incremented string.
!
  implicit none

  character c
  integer change
  integer digit
  character ( len = * ) file_name
  integer i
  integer lens

  lens = len_trim ( file_name )

  if ( lens <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_NAME_INC - Fatal error!'
    write ( *, '(a)' ) '  The input string is empty.'
    stop
  end if

  change = 0

  do i = lens, 1, -1

    c = file_name(i:i)

    if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      change = change + 1

      digit = ichar ( c ) - 48
      digit = digit + 1

      if ( digit == 10 ) then
        digit = 0
      end if

      c = char ( digit + 48 )

      file_name(i:i) = c

      if ( c /= '0' ) then
        return
      end if

    end if

  end do

  if ( change == 0 ) then
    file_name = ' '
    return
  end if

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT, the free unit number.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen
 
  iunit = 0
 
  do i = 1, 99
 
    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )
 
      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if
 
  end do

  return
end
subroutine p06_barray ( b )

!*****************************************************************************80
!
!! P06_BARRAY sets the B array.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 August 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) B(5,8), the array of coefficients for the linear
!    part of the aircraft stability function.
!
  implicit none

  real    ( kind = 8 ) b(5,8)
  real    ( kind = 8 ), parameter, dimension ( 5, 8 ) :: b_save = reshape ( (/ &
    -3.933D+00,  0.0D+00,    0.002D+00,  0.0D+00,    0.0D+00, &
     0.107D+00, -0.987D+00,  0.0D+00,    1.0D+00,    0.0D+00, &
     0.126D+00,  0.0D+00,   -0.235D+00,  0.0D+00,   -1.0D+00, &
     0.0D+00,  -22.95D+00,   0.0D+00,   -1.0D+00,    0.0D+00, &
    -9.99D+00,   0.0D+00,    5.67D+00,   0.0D+00,   -0.196D+00, &
     0.0D+00,  -28.37D+00,   0.0D+00,   -0.168D+00,  0.0D+00, &
   -45.83D+00,   0.0D+00,   -0.921D+00,  0.0D+00,   -0.0071D+00, &
    -7.64D+00,   0.0D+00,   -6.51D+00,   0.0D+00,    0.0D+00 /), (/ 5, 8 /) )

  b(1:5,1:8) = b_save(1:5,1:8)

  return
end
subroutine p06_carray ( c )

!*****************************************************************************80
!
!! P06_CARRAY sets the C array.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    17 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) C(5,8,8), the array of coefficients for the nonlinear
!    part of the aircraft stability function.
!
  implicit none

  real    ( kind = 8 ) c(5,8,8)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k

  c(1:5,1:8,1:8) = 0.0D+00

  c(1,2,3) = -   0.727D+00
  c(1,3,4) =     8.39D+00
  c(1,4,5) = - 684.4D+00
  c(1,4,7) =  + 63.5D+00

  c(2,1,3) =   + 0.949D+00
  c(2,1,5) =   + 0.173D+00

  c(3,1,2) =   - 0.716D+00
  c(3,1,4) =   - 1.578D+00
  c(3,4,7) =   + 1.132D+00

  c(4,1,5) =   - 1.0D+00

  c(5,1,4) =   + 1.0D+00

  return
end
subroutine p06_fun ( option, nvar, x, fx )

!*****************************************************************************80
!
!! P06_FUN evaluates the function for problem 6.
!
!  Title:
!
!    The aircraft stability problem.
!
!  Description:
!
!    The equations describe the behavior of an aircraft under the
!    control of a pilot.  The variables are
!
!      X(1) = roll
!      X(2) = pitch
!      X(3) = yaw
!      X(4) = angle of attack
!      X(5) = sideslip
!      X(6) = elevator
!      X(7) = aileron
!      X(8) = rudder
!
!    The function is of the following form
!
!    For indices I=1 through 5,
!
!      F(I) = SUM ( 1 <= J <= 8 ) B(I,J) * X(J)
!           + SUM ( 1 <= J <= 8, 1 <= K <= 8 ) C(I,J,K) * X(J) * X(K)
!
!    with the last two equations fixing the values of the elevator
!    and rudder:
!
!      F(6) = X(6) - value
!      F(7) = X(8)
!
!  Options:
!
!    There are five options, which vary in the value they fix the
!    elevator value in function 6:
!
!      Option   Constant Elevator value
!
!       1        -0.050
!       2        -0.008
!       3         0.0D+00
!       4         0.05
!       5         0.1
!
!  Special points:
!
!    Melhem lists the following limit points in X(7)
!    (note that melhem has B(4,1)=1.0, B(4,2)=0.0)
!
!       X(1)    X(2)     X(3)    X(4)    X(5)    X(6)   X(7)     X(8)
!
!    -2.9691  0.8307  -0.0727  0.4102 -0.2688 -0.05   0.5091   0.0
!    -2.8158 -0.1748  -0.0894  0.0263  0.0709 -0.008  0.2044   0.0
!    -3.7571 -0.6491  -0.3935  0.0918  0.1968 -0.008 -0.0038   0.0
!    -4.1637  0.0922  -0.0926  0.0224 -0.0171 -0.008  0.3782   0.0
!    -2.5839 -0.2212  -0.0540  0.0135  0.0908  0.0    0.1860   0.0
!    -3.9007 -1.1421  -0.5786  0.1328  0.3268  0.0   -0.5070   0.0
!    -2.3610 -0.7236   0.0327 -0.0391  0.2934  0.05   0.2927   0.0
!    -2.2982  1.4033   0.0632 -0.0793  0.5833  0.10   0.5833   0.0
!
!    Rheinboldt lists the following limit points in X(7), with
!    B(4,1)=0.0, B(4,2)=1.0:
!
!       X(1)    X(2)     X(3)    X(4)    X(5)    X(6)   X(7)     X(8)
!
!     2.9648  0.8255   0.0736  0.0413  0.2673 -0.050 -0.0504   0.0
!     2.8173 -0.1762   0.0899  0.0264 -0.0714 -0.008 -0.2049   0.0
!     3.7579 -0.6554   0.3865  0.0925 -0.1986 -0.008  0.0062   0.0
!     4.1638  0.0891   0.0948  0.0228  0.1623 -0.008 -0.3776   0.0
!     2.5873 -0.2235   0.0546  0.0136 -0.0916  0.000 -0.1869   0.0
!     3.9005 -1.1481   0.5815  0.1335 -0.3285  0.000  0.5101   0.0
!     2.3639 -0.7297  -0.3160 -0.0387 -0.2958  0.050 -0.2957   0.0
!     2.2992 -1.4102  -0.0618 -0.0790 -0.5862  0.100 -0.6897   0.0
!
!    Rheinboldt lists the following bifurcation points:
!
!       X(1)    X(2)     X(3)    X(4)    X(5)    X(6)    X(7)    X(8)
!
!     4.482   0.1632   0.0237  0.0062  0.0352 -0.0006 -0.3986  0.0
!     3.319  -0.1869   0.1605  0.0437 -0.0688 -0.0125 -0.2374  0.0
!     4.466   0.1467   0.0404  0.0097  0.0308 -0.0061 -0.3995  0.0
!    -3.325   0.1880  -0.1614  0.0439  0.0691 -0.0124  0.2367  0.0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    13 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Raman Mehra, William Kessel, James Carroll,
!    Global stability and contral analysis of aircraft at high angles of attack,
!    Technical Report CR-215-248-1, -2, -3,
!    Office of Naval Research, June 1977.
!
!    Rami Melhem, Werner Rheinboldt,
!    A Comparison of Methods for Determining Turning Points of Nonlinear Equations,
!    Computing,
!    Volume 29, Number 3, September 1982, pages 201-226.
!
!    Albert Schy, Margery Hannah,
!    Prediction of Jump Phenomena in Roll-coupled Maneuvers of Airplanes,
!    Journal of Aircraft,
!    Volume 14, Number 4, 1977,  pages 375-382.
!
!    John Young, Albert Schy, Katherine Johnson,,
!    Prediction of Jump Phenomena in Aircraft Maneuvers, Including
!    Nonlinear Aerodynamic Effects,
!    Journal of Guidance and Control,
!    Volume 1, Number 1, 1978, pages 26-31.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Input, real ( kind = 8 ) X(NVAR), the argument of the function.
!
!    Output, real ( kind = 8 ) FX(NVAR-1), the value of the function at X.
!
  implicit none

  integer ( kind = 4 ) nvar

  real    ( kind = 8 ) b(5,8)
  real    ( kind = 8 ) c(5,8,8)
  real    ( kind = 8 ) fx(nvar-1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) option
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  real    ( kind = 8 ) val
  real    ( kind = 8 ) x(nvar)
!
!  Compute the linear term.
!
  call p06_barray ( b )

  fx(1:5) = matmul ( b(1:5,1:8), x(1:8) )
!
!  Compute the nonlinear terms.
!
  call p06_carray ( c )

  do i = 1, 5
    do j = 1, 8
      do k = 1, 8
        fx(i) = fx(i) + c(i,j,k) * x(j) * x(k)
      end do
    end do
  end do
!
!  Set function values for two fixed variables.
!
  if ( option == 1 ) then
    val = - 0.050D+00
  else if ( option == 2 ) then
    val = - 0.008D+00
  else if ( option == 3 ) then
    val =   0.000D+00
  else if ( option == 4 ) then
    val =   0.050D+00
  else if ( option == 5 ) then
    val =   0.100D+00
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P06_FUN - Fatal error!'
    write ( *, '(a,i8)' ) '  Unrecognized option  = ', option
    stop
  end if

  fx(6) = x(6) - val
  fx(7) = x(8)

  return
end
subroutine p06_jac ( option, nvar, x, jac )

!*****************************************************************************80
!
!! P06_JAC evaluates the jacobian for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    13 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Input, real ( kind = 8 ) X(NVAR), the argument of the jacobian.
!
!    Output, real ( kind = 8 ) JAC(NVAR,NVAR), the jacobian matrix evaluated
!    at X.  The NVAR-th row is not set by this routine.
!
  implicit none

  integer ( kind = 4 ) nvar

  real    ( kind = 8 ) b(5,8)
  real    ( kind = 8 ) c(5,8,8)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) option
  integer ( kind = 4 ) j
  real    ( kind = 8 ) jac(nvar,nvar)
  integer ( kind = 4 ) k
  real    ( kind = 8 ) x(nvar)

  jac(1:nvar,1:nvar) = 0.0D+00
!
!  Set the jacobian to the linear coefficients.
!
  call p06_barray ( b )

  jac(1:5,1:8) = b(1:5,1:8)
!
!  Add the nonlinear terms.
!
  call p06_carray ( c )

  do i = 1, 5
    do j = 1, 8
      do k = 1, 8
        jac(i,j) = jac(i,j) + ( c(i,j,k) + c(i,k,j) ) * x(k)
      end do
    end do
  end do
!
!  Constraint equations.
!
  jac(6,6) = 1.0D+00
  jac(7,8) = 1.0D+00

  return
end
subroutine p06_newton ( option, nvar, x, par_index, status )

!*****************************************************************************80
!
!! P06_NEWTON applies Newton's method to an approximate root.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    06 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Input/output, real ( kind = 8 ) X(NVAR).
!    On input, the starting point of Newton's method.
!    On output, an improved estimate of the root of F(X)=0, if the
!    algorithm converged.
!
!    Input, integer ( kind = 4 ) PAR_INDEX, the index of the parameter 
!    to be fixed.  This variable should be between 1 and NVAR.  However, the user can
!    set it to 0, indicating that the program should make an intelligent
!    choice for the index.
!
!    Output, integer ( kind = 4 ) STATUS, the status of the iteration.
!    -3, the full number of steps was taken without convergence.
!        (however, the output X might be CLOSE to a good solution).
!    -2, the iteration seemed to be diverging, and was halted.
!    -1, the jacobian was singular, and the iteration was halted.
!     nonnegative, the convergence test was satisfied, and this is the
!        number of steps taken (possibly 0).
!
  implicit none

  integer ( kind = 4 ) nvar

  real    ( kind = 8 ) fx(nvar)
  real    ( kind = 8 ), parameter :: FX_ABS_TOL = 0.000001D+00
  real    ( kind = 8 ) fx_max
  real    ( kind = 8 ) fx_max_init
  integer ( kind = 4 ) ipar
  integer ( kind = 4 ) ipivot(nvar)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) info
  integer ( kind = 4 ) option
  integer ( kind = 4 ) it
  integer ( kind = 4 ), parameter :: IT_MAX = 20
  integer ( kind = 4 ) j
  real    ( kind = 8 ) jac(nvar,nvar)
  integer ( kind = 4 ) job
  integer ( kind = 4 ) par_index
  real    ( kind = 8 ) par_value
  integer ( kind = 4 ) status
  logical, parameter :: VERBOSE = .false.
  real    ( kind = 8 ) x(nvar)

  if ( par_index < 1 .or. nvar < par_index ) then
    call p06_par_index ( option, nvar, x, par_index )

    if ( VERBOSE ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a,i8,a)' ) '  Iteration will hold index ', par_index, ' fixed.'
    end if

  end if
  par_value = x(par_index)

  if ( VERBOSE ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_NEWTON'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Step   F(X)'
    write ( *, '(a)' ) ' '
  end if

  do it = 0, IT_MAX
!
!  Compute the function value.
!
    call p06_fun ( option, nvar, x, fx )
    fx(nvar) = x(par_index) - par_value
!
!  Compute the norm of the function value.
!
    fx_max = maxval ( abs ( fx(1:nvar) ) )

    if ( VERBOSE ) then
      write ( *, '(2x,i8,2x,g14.6)' ) it, fx_max
    end if

    if ( it == 0 ) then
      fx_max_init = fx_max
    end if
!
!  If the function norm is small enough, return.
!
    if ( abs ( fx_max ) < FX_ABS_TOL ) then
      status = it
      exit
    end if
!
!  If the function norm seems to be exploding, halt.
!
    if ( 1000.0 * fx_max_init < abs ( fx_max ) ) then
      status = -2
      exit
    end if

    if ( it == IT_MAX ) then
      status = -3
      exit
    end if
!
!  Compute the jacobian.
!
    call p06_jac ( option, nvar, x, jac )

    jac(nvar,1:nvar) = 0.0D+00
    jac(nvar,par_index) = 1.0D+00
!
!  Factor the jacobian.
!
    call sge_fa ( nvar, nvar, jac, ipivot, info )

    if ( info /= 0 ) then
      status = -1
      exit
    end if
!
!  Solve the system JAC * DX = FX
!
    job = 0
    call sge_sl ( nvar, nvar, jac, ipivot, fx, job )
!
!  Update X = X - DX.
!
    x(1:nvar) = x(1:nvar) - fx(1:nvar)

  end do

  return
end
subroutine p06_nvar ( option, nvar )

!*****************************************************************************80
!
!! P06_NVAR sets the number of variables for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Output, integer ( kind = 4 ) NVAR, the number of variables.
!
  implicit none

  integer ( kind = 4 ) option
  integer ( kind = 4 ) nvar

  nvar = 8

  return
end
subroutine p06_option_num ( option_num )

!*****************************************************************************80
!
!! P06_OPTION_NUM returns the number of options for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) OPTION_NUM, the number of options.
!
  implicit none

  integer ( kind = 4 ) option_num

  option_num = 5

  return
end
subroutine p06_par_index ( option, nvar, x, par_index )

!*****************************************************************************80
!
!! P00_PAR_INDEX chooses the index of the continuation parameter.
!
!  Discussion:
!
!    Given the NVAR-dimensional point X, the (NVAR-1)-dimensional function 
!    F(X), and the NVAR-1 by NVAR jacobian matrix, let the NVAR-dimensional 
!    vector TAN be any null vector of JAC.
!
!      JAC * TAN = 0
!
!    Choose PAR_INDEX to be the index of TAN of maximum absolute value.
!     
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 September 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer OPTION, the option index.
!
!    Input, integer NVAR, the number of variables.
!
!    Input, real X(NVAR), the starting point of Newton's method.
!
!    Output, integer PAR_INDEX, the index of the parameter to be held fixed.
!    This variable will be between 1 and NVAR.  It is the index of the variable
!    which is currently changing most rapidly along the curve F(X) = 0.
!
  implicit none

  integer ( kind = 4 ) nvar

  integer ( kind = 4 ) option
  integer ( kind = 4 ) par_index
  real    ( kind = 8 ) tan(nvar)
  real    ( kind = 8 ) x(nvar)

  call p06_tan ( option, nvar, x, tan )

  call r8vec_amax_index ( nvar, tan, par_index )

  return
end
subroutine p06_start ( option, nvar, x )

!*****************************************************************************80
!
!! P06_START returns a starting point for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    13 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Output, real ( kind = 8 ) X(NVAR), the starting point.
!
  implicit none

  integer ( kind = 4 ) nvar

  integer ( kind = 4 ) i
  integer ( kind = 4 ) option
  real    ( kind = 8 ) x(nvar)

  if ( option == 1 ) then

    x(1:nvar) = (/ &
       1.06001162985175758D-03, &
       5.12061216467178115D-02, &
       5.79953409787390485D-05, &
       5.96060845777059631D-02, &
       2.64683802731226678D-05, &
      -5.00000000000000000D-02, &
       0.00000000000000000D+00, &  
       0.00000000000000000D+00 /)

  else if ( option == 2 ) then

    x(1:nvar) = (/ &
       0.000001548268247D+00, &
       0.008192973225663D+00, &
      -0.000000682134573D+00, &
       0.009536973221178D+00, &
       0.000002896734870D+00, &
      -0.008000000000000D+00, &
       0.000018188778989D+00, &
       0.000000000000000D+00 /)

  else if ( option == 3 ) then

    x(1:nvar) = (/ &
       0.0D+00, &
       0.0D+00, &
       0.0D+00, &
       0.0D+00, &
       0.0D+00, &
       0.0D+00, &
       0.0D+00, &
       0.0D+00 /)

  else if ( option == 4 ) then

    x(1:nvar) = (/ &
      -0.000010655314069D+00, &
      -0.051206082422980D+00, &
       0.000005600187501D+00, &
      -0.059606082643400D+00, &
      -0.000020891016199D+00, &
       0.050000000000000D+00, &
      -0.000122595323216D+00, &
       0.000000000000000D+00 /)

  else if ( option == 5 ) then

    x(1:nvar) = (/ &
      -0.000027083319493D+00, &
      -0.102412164106124D+00, &
       0.000014540858026D+00, &
      -0.119212165322433D+00, &
      -0.000048014067202D+00, &
       0.100000000000000D+00, &
      -0.000267808407544D+00, &
       0.000000000000000D+00 /)

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P06_START - Fatal error!'
    write ( *, '(a,i8)' ) '  Unrecognized option  = ', option
    stop
  end if

  return
end
subroutine p06_stepsize ( option, h, hmin, hmax )

!*****************************************************************************80
!
!! P06_STEPSIZE returns step sizes for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    24 September 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Output, real ( kind = 8 ) H, HMIN, HMAX, suggested values for the
!    initial step, the minimum step, and the maximum step.
!
  implicit none

  real    ( kind = 8 ) h
  real    ( kind = 8 ) hmax
  real    ( kind = 8 ) hmin
  integer ( kind = 4 ) option

  h =   - 0.250D+00
  hmin =  0.001D+00
  hmax =  0.500D+00

  return
end
subroutine p06_tan ( option, nvar, x, tan )

!*****************************************************************************80
!
!! P06_TAN determines a tangent vector at X.
!
!  Discussion:
!
!    If X is a solution of F(Y) = 0, then the vector TAN
!    is tangent to the curve of solutions at X.
!
!    If X is not a solution of F(Y) = 0, then the vector TAN
!    is tangent to the curve F(Y) = F(X) at X.
!
!    The vector will have unit euclidean norm.
!
!    The sign of TAN will be chosen so that the determinant
!    of F'(X) augmented with a final row equal to TAN will be positive.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    03 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Input, real ( kind = 8 ) X(NVAR), the evaluation point.
!
!    Output, real ( kind = 8 ) TAN(NVAR), a tangent vector at X.
!
  implicit none

  integer ( kind = 4 ) nvar

  real    ( kind = 8 ) jac_det
  real    ( kind = 8 ) jac(nvar,nvar)
  real    ( kind = 8 ), allocatable, dimension ( :, : ) :: nullspace
  integer ( kind = 4 ) nullspace_size
  integer ( kind = 4 ) option
  real    ( kind = 8 ) tan(nvar)
  real    ( kind = 8 ) tan_norm
  real    ( kind = 8 ) x(nvar)
!
!  Compute the jacobian.
!
  call p06_jac ( option, nvar, x, jac )
!
!  Compute the nullspace size.
!
  call r8mat_nullspace_size ( nvar, nvar, jac, nullspace_size )

  if ( nullspace_size < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_TAN - Fatal error!'
    write ( *, '(a)' ) '  The matrix seems to have no nullspace.'
    write ( *, '(a)' ) '  The tangent vector could not be computed.'
    stop
  end if
!
!  Compute the nullspace.
!
  allocate ( nullspace(1:nvar,1:nullspace_size) )

  call r8mat_nullspace ( nvar, nvar, jac, nullspace_size, nullspace )

  tan(1:nvar) = nullspace(1:nvar,1)

  deallocate ( nullspace )
!
!  Choose the sign of TAN by the determinant condition.
!
  jac(nvar,1:nvar) = tan(1:nvar)

  call r8mat_det ( nvar, jac, jac_det )

  if ( jac_det < 0.0D+00 ) then
    tan(1:nvar) = - tan(1:nvar)
  end if

  tan_norm = sqrt ( sum ( tan(1:nvar)**2 ) )

  tan(1:nvar) = tan(1:nvar) / tan_norm

  return
end
subroutine p06_title ( option, title )

!*****************************************************************************80
!
!! P06_TITLE sets the title for problem 6.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!    TITLE will never be longer than 80 characters.
!
  implicit none

  integer ( kind = 4 ) option
  character ( len = * ) title

  title = '???'

  if ( option == 1 ) then
    title = 'Aircraft function, x(6) = - 0.050.'
  else if ( option == 2 ) then
    title = 'Aircraft function, x(6) = - 0.008.'
  else if ( option == 3 ) then
    title = 'Aircraft function, x(6) =   0.000.'
  else if ( option == 4 ) then
    title = 'Aircraft function, x(6) = + 0.050.'
  else if ( option == 5 ) then
    title = 'Aircraft function, x(6) = + 0.100.'
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P06_TITLE - Fatal error!'
    write ( *, '(a,i8)' ) '  Unrecognized option  = ', option
    stop
  end if

  return
end
function r8_sign ( x )

!*****************************************************************************80
!
!! R8_SIGN returns the sign of an R8.
!
!  Discussion:
!
!    value = -1 if X < 0;
!    value =  0 if X => 0.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    27 March 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the number whose sign is desired.
!
!    Output, real ( kind = 8 ) R8_SIGN, the sign of X:
!
  implicit none

  real    ( kind = 8 ) r8_sign
  real    ( kind = 8 ) x

  if ( x < 0.0D+00 ) then
    r8_sign = -1.0D+00
  else
    r8_sign = +1.0D+00
  end if

  return
end
subroutine r8_swap ( x, y )

!*****************************************************************************80
!
!! R8_SWAP swaps two R8's.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none

  real    ( kind = 8 ) x
  real    ( kind = 8 ) y
  real    ( kind = 8 ) z

  z = x
  x = y
  y = z

  return
end
subroutine r8mat_det ( n, a, det )

!*****************************************************************************80
!
!! R8MAT_DET computes the determinant of an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an array of R8 values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    Original FORTRAN77 version by Helmut Spaeth
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Helmut Spaeth,
!    Cluster Analysis Algorithms
!    for Data Reduction and Classification of Objects,
!    Ellis Horwood, 1980, page 125-127.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order of the matrix.
!
!    Input, real ( kind = 8 ) A(N,N), the matrix whose determinant is desired.
!
!    Output, real ( kind = 8 ) DET, the determinant of the matrix.
!
  implicit none

  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(n,n)
  real    ( kind = 8 ) b(n,n)
  real    ( kind = 8 ) det
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) piv(1)
  real    ( kind = 8 ) t

  b(1:n,1:n) = a(1:n,1:n)

  det = 1.0D+00

  do k = 1, n

    piv = maxloc ( abs ( b(k:n,k) ) )

    m = piv(1) + k - 1

    if ( m /= k ) then
      det = - det
      t      = b(m,k)
      b(m,k) = b(k,k)
      b(k,k) = t
    end if

    det = det * b(k,k)

    if ( b(k,k) /= 0.0D+00 ) then

      b(k+1:n,k) = -b(k+1:n,k) / b(k,k)

      do j = k + 1, n
        if ( m /= k ) then
          t      = b(m,j)
          b(m,j) = b(k,j)
          b(k,j) = t
        end if
        b(k+1:n,j) = b(k+1:n,j) + b(k+1:n,k) * b(k,j)
      end do

    end if

  end do

  return
end
subroutine r8mat_nullspace ( m, n, a, nullspace_size, nullspace )

!*****************************************************************************80
!
!! R8MAT_NULLSPACE computes the nullspace of a matrix.
!
!  Discussion:
!
!    Let A be an MxN matrix.
!
!    If X is an N-vector, and A*X = 0, then X is a null vector of A.
!
!    The set of all null vectors of A is called the nullspace of A.
!
!    The 0 vector is always in the null space.
!
!    If the 0 vector is the only vector in the nullspace of A, then A
!    is said to have maximum column rank.  (Because A*X=0 can be regarded
!    as a linear combination of the columns of A).  In particular, if A
!    is square, and has maximum column rank, it is nonsingular.
!
!    The dimension of the nullspace is the number of linearly independent
!    vectors that span the nullspace.  If A has maximum column rank,
!    its nullspace has dimension 0.  
!
!    This routine uses the reduced row echelon form of A to determine
!    a set of NULLSPACE_SIZE independent null vectors.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    02 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns of
!    the matrix A.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix to be analyzed. 
!
!    Input, integer ( kind = 4 ) NULLSPACE_SIZE, the size of the nullspace.
!
!    Output, real ( kind = 8 ) NULLSPACE(N,NULLSPACE_SIZE), vectors that
!    span the nullspace.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nullspace_size

  real    ( kind = 8 ) a(m,n)
  integer ( kind = 4 ) col(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i2
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  real    ( kind = 8 ) nullspace(n,nullspace_size)
  integer ( kind = 4 ) row(m)
  real    ( kind = 8 ) rref(m,n)
!
!  Make a copy of A.
!
  rref(1:m,1:n) = a(1:m,1:n)
!
!  Get the reduced row echelon form of A.
!
  call r8mat_rref ( m, n, rref )
!
!  Note in ROW the columns of the leading nonzeros.
!  COL(J) = +J if there is a leading 1 in that column, and -J otherwise.
!
  row(1:m) = 0

  do j = 1, n
    col(j) = - j
  end do

  do i = 1, m
    do j = 1, n
      if ( rref(i,j) == 1.0D+00 ) then
        row(i) = j
        col(j) = j
        exit
      end if
    end do
  end do

  nullspace(1:n,1:nullspace_size) = 0.0D+00

  j2 = 0
!
!  If column J does not contain a leading 1, then it contains
!  information about a null vector.
!
  do j = 1, n

    if ( col(j) < 0 ) then
 
      j2 = j2 + 1

      do i = 1, m
        if ( rref(i,j) /= 0.0D+00 ) then
          i2 = row(i)
          nullspace(i2,j2) = - rref(i,j)
        end if
      end do

      nullspace(j,j2) = 1.0D+00

    end if
    
  end do

  return
end
subroutine r8mat_nullspace_size ( m, n, a, nullspace_size )

!*****************************************************************************80
!
!! R8MAT_NULLSPACE_SIZE computes the size of the nullspace of a matrix.
!
!  Discussion:
!
!    Let A be an MxN matrix.
!
!    If X is an N-vector, and A*X = 0, then X is a null vector of A.
!
!    The set of all null vectors of A is called the nullspace of A.
!
!    The 0 vector is always in the null space.
!
!    If the 0 vector is the only vector in the nullspace of A, then A
!    is said to have maximum column rank.  (Because A*X=0 can be regarded
!    as a linear combination of the columns of A).  In particular, if A
!    is square, and has maximum column rank, it is nonsingular.
!
!    The dimension of the nullspace is the number of linearly independent
!    vectors that span the nullspace.  If A has maximum column rank,
!    its nullspace has dimension 0.  
!
!    This routine ESTIMATES the dimension of the nullspace.  Cases of
!    singularity that depend on exact arithmetic will probably be missed.
!
!    The nullspace will be estimated by counting the leading 1's in the
!    reduced row echelon form of A, and subtracting this from N.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    02 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns of
!    the matrix A.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix to be analyzed. 
!
!    Output, integer ( kind = 4 ) NULLSPACE_SIZE, the estimated size 
!    of the nullspace.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(m,n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) leading
  integer ( kind = 4 ) nullspace_size
  real    ( kind = 8 ) rref(m,n)
!
!  Get the reduced row echelon form of A.
!
  rref(1:m,1:n) = a(1:m,1:n)

  call r8mat_rref ( m, n, rref )
!
!  Count the leading 1's in A.
!
  leading = 0
  do i = 1, m
    do j = 1, n
      if ( rref(i,j) == 1.0D+00 ) then
        leading = leading + 1
        exit
      end if
    end do
  end do

  nullspace_size = n - leading

  return
end
subroutine r8mat_rref ( m, n, a )

!*****************************************************************************80
!
!! R8MAT_RREF computes the reduced row echelon form of a matrix.
!
!  Discussion:
!
!    A matrix is in row echelon form if:
!
!    * The first nonzero entry in each row is 1.
!
!    * The leading 1 in a given row occurs in a column to
!      the right of the leading 1 in the previous row.
!
!    * Rows which are entirely zero must occur last.
!
!    The matrix is in reduced row echelon form if, in addition to
!    the first three conditions, it also satisfies:
!
!    * Each column containing a leading 1 has no other nonzero entries.
!
!  Example:
!
!    Input matrix:
!
!     1.0  3.0  0.0  2.0  6.0  3.0  1.0
!    -2.0 -6.0  0.0 -2.0 -8.0  3.0  1.0
!     3.0  9.0  0.0  0.0  6.0  6.0  2.0
!    -1.0 -3.0  0.0  1.0  0.0  9.0  3.0
!
!    Output matrix:
!
!     1.0  3.0  0.0  0.0  2.0  0.0  0.0
!     0.0  0.0  0.0  1.0  2.0  0.0  0.0
!     0.0  0.0  0.0  0.0  0.0  1.0  0.3
!     0.0  0.0  0.0  0.0  0.0  0.0  0.0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns of 
!    the matrix A.
!
!    Input/output, real ( kind = 8 ) A(M,N).  On input, the matrix to be
!    analyzed.  On output, the RREF form of the matrix.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(m,n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) lead
  integer ( kind = 4 ) r
  real    ( kind = 8 ) temp

  lead = 1

  do r = 1, m

    if ( n < lead ) then
      exit
    end if

    i = r

    do while ( a(i,lead) == 0.0 )

      i = i + 1

      if ( m < i ) then
        i = r
        lead = lead + 1
        if ( n < lead ) then
          lead = -1
          exit
        end if
      end if

    end do

    if ( lead < 0 ) then
      exit
    end if

    do j = 1, n
      temp   = a(i,j)
      a(i,j) = a(r,j)
      a(r,j) = temp
    end do

    a(r,1:n) = a(r,1:n) / a(r,lead)

    do i = 1, m
      if ( i /= r ) then
        a(i,1:n) = a(i,1:n) - a(i,lead) * a(r,1:n)
      end if
    end do

    lead = lead + 1

  end do

  return
end
subroutine r8vec_amax_index ( n, a, amax_index )

!*****************************************************************************80
!
!! R8VEC_AMAX_INDEX returns the index of the maximum absolute value in an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the array.
!
!    Input, real ( kind = 8 ) A(N), the array.
!
!    Output, integer ( kind = 4 ) AMAX_INDEX, the index of the entry of largest magnitude.
!
  implicit none

  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(n)
  real    ( kind = 8 ) amax
  integer ( kind = 4 ) amax_index
  integer ( kind = 4 ) i

  if ( n <= 0 ) then

    amax_index = -1

  else

    amax_index = 1
    amax = abs ( a(1) )

    do i = 2, n
      if ( amax < abs ( a(i) ) ) then
        amax_index = i
        amax = abs ( a(i) )
      end if
    end do

  end if

  return
end
subroutine sge_check ( lda, m, n, ierror )

!*****************************************************************************80
!
!! SGE_CHECK checks the dimensions of a general matrix.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) LDA, the leading dimension of the array.
!    LDA must be at least M.
!
!    Input, integer ( kind = 4 ) M, the number of rows of the matrix.
!    M must be positive.
!
!    Input, integer ( kind = 4 ) N, the number of columns of the matrix.
!    N must be positive.
!
!    Output, integer ( kind = 4 ) IERROR, reports whether any errors were detected.
!    IERROR is set to 0 before the checks are made, and then:
!    IERROR = IERROR + 1 if LDA is illegal;
!    IERROR = IERROR + 2 if M is illegal;
!    IERROR = IERROR + 4 if N is illegal.
!
  implicit none

  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) lda
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  ierror = 0

  if ( lda < m ) then
    ierror = ierror + 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_CHECK - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal LDA = ', lda
    stop
  end if

  if ( m < 1 ) then
    ierror = ierror + 2
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_CHECK - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal M = ', m
    stop
  end if

  if ( n < 1 ) then
    ierror = ierror + 4
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_CHECK - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal N = ', n
    stop
  end if

  return
end
subroutine sge_fa ( lda, n, a, pivot, info )

!*****************************************************************************80
!
!! SGE_FA factors a general matrix.
!
!  Discussion:
!
!    SGE_FA is a simplified version of the LINPACK routine SGEFA.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    26 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) LDA, the leading dimension of the array.
!    LDA must be at least N.
!
!    Input, integer ( kind = 4 ) N, the order of the matrix.
!    N must be positive.
!
!    Input/output, real ( kind = 8 ) A(LDA,N), the matrix to be factored.
!    On output, A contains an upper triangular matrix and the multipliers
!    which were used to obtain it.  The factorization can be written
!    A = L * U, where L is a product of permutation and unit lower
!    triangular matrices and U is upper triangular.
!
!    Output, integer ( kind = 4 ) PIVOT(N), a vector of pivot indices.
!
!    Output, integer ( kind = 4 ) INFO, singularity flag.
!    0, no singularity detected.
!    nonzero, the factorization failed on the INFO-th step.
!
  implicit none

  integer ( kind = 4 ) lda
  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(lda,n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) info
  integer ( kind = 4 ) pivot(n)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  real    ( kind = 8 ) t
!
!  Check the dimensions.
!
  call sge_check ( lda, n, n, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_FA - Fatal error!'
    write ( *, '(a)' ) '  Illegal dimensions.'
    stop
  end if

  info = 0

  do k = 1, n - 1
!
!  Find L, the index of the pivot row.
!
    l = k
    do i = k + 1, n
      if ( abs ( a(l,k) ) < abs ( a(i,k) ) ) then
        l = i
      end if
    end do

    pivot(k) = l
!
!  If the pivot index is zero, the algorithm has failed.
!
    if ( a(l,k) == 0.0D+00 ) then
      info = k
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SGE_FA - Fatal error!'
      write ( *, '(a,i8)' ) '  Zero pivot on step ', info
      return
    end if
!
!  Interchange rows L and K if necessary.
!
    if ( l /= k ) then
      call r8_swap ( a(l,k), a(k,k) )
    end if
!
!  Normalize the values that lie below the pivot entry A(K,K).
!
    a(k+1:n,k) = -a(k+1:n,k) / a(k,k)
!
!  Row elimination with column indexing.
!
    do j = k + 1, n

      if ( l /= k ) then
        call r8_swap ( a(l,j), a(k,j) )
      end if

      a(k+1:n,j) = a(k+1:n,j) + a(k+1:n,k) * a(k,j)

    end do

  end do

  pivot(n) = n

  if ( a(n,n) == 0.0D+00 ) then
    info = n
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_FA - Fatal error!'
    write ( *, '(a,i8)' ) '  Zero pivot on step ', info
  end if

  return
end
subroutine sge_sl ( lda, n, a, pivot, b, job )

!*****************************************************************************80
!
!! SGE_SL solves a system factored by SGE_FA.
!
!  Discussion:
!
!    SGE_SL is a simplified version of the LINPACK routine SGESL.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    04 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) LDA, the leading dimension of the array.
!    LDA must be at least N.
!
!    Input, integer ( kind = 4 ) N, the order of the matrix.
!    N must be positive.
!
!    Input, real ( kind = 8 ) A(LDA,N), the LU factors from SGE_FA.
!
!    Input, integer ( kind = 4 ) PIVOT(N), the pivot vector from SGE_FA.
!
!    Input/output, real ( kind = 8 ) B(N).
!    On input, the right hand side vector.
!    On output, the solution vector.
!
!    Input, integer ( kind = 4 ) JOB, specifies the operation.
!    0, solve A * x = b.
!    nonzero, solve transpose ( A ) * x = b.
!
  implicit none

  integer ( kind = 4 ) lda
  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(lda,n)
  real    ( kind = 8 ) b(n)
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) pivot(n)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) job
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  real    ( kind = 8 ) t
!
!  Check the dimensions.
!
  call sge_check ( lda, n, n, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SGE_SL - Fatal error!'
    write ( *, '(a)' ) '  Illegal dimensions.'
    stop
  end if
!
!  Solve A * x = b.
!
  if ( job == 0 ) then
!
!  Solve PL * Y = B.
!
    do k = 1, n - 1

      l = pivot(k)

      if ( l /= k ) then
        call r8_swap ( b(l), b(k) )
      end if

      b(k+1:n) = b(k+1:n) + a(k+1:n,k) * b(k)

    end do
!
!  Solve U * X = Y.
!
    do k = n, 1, -1
      b(k) = b(k) / a(k,k)
      b(1:k-1) = b(1:k-1) - a(1:k-1,k) * b(k)
    end do
!
!  Solve transpose ( A ) * X = B.
!
  else
!
!  Solve transpose ( U ) * Y = B.
!
    do k = 1, n
      b(k) = ( b(k) - dot_product ( b(1:k-1), a(1:k-1,k) ) ) / a(k,k)
    end do
!
!  Solve transpose ( PL ) * X = Y.
!
    do k = n - 1, 1, -1

      b(k) = b(k) + dot_product ( b(k+1:n), a(k+1:n,k) )

      l = pivot(k)

      if ( l /= k ) then
        call r8_swap ( b(l), b(k) )
      end if

    end do

  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  character ( len = 8 ) date
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  character ( len = 10 )  time
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end