subroutine dpoly_val ( n, p, x, pval )

!*****************************************************************************80
!
!! DPOLY_VAL evaluates a double precision polynomial.
!
!  Modified:
!
!    08 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the polynomial.
!
!    Input, double precision PCOF(0:N), the polynomial coefficients.
!    P(I) is the coefficient of X**I.
!
!    Input, double precision X, the evaluation point.
!
!    Output, double precision PVAL, the value of the polynomial at X.
!
  implicit none

  integer n

  integer i
  double precision p(0:n)
  double precision pval
  double precision x

  pval = p(0)
  do i = 1, n
    pval = pval + p(i) * x**i
  end do

  return
end
subroutine dpoly_val_horner ( n, p, x, pval )

!*****************************************************************************80
!
!! DPOLY_VAL_HORNER evaluates a double precision polynomial using Horner's method.
!
!  Modified:
!
!    08 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the polynomial.
!
!    Input, double precision PCOF(0:N), the polynomial coefficients.
!    P(I) is the coefficient of X**I.
!
!    Input, double precision X, the evaluation point.
!
!    Output, double precision PVAL, the value of the polynomial at X.
!
  implicit none

  integer n

  integer i
  double precision p(0:n)
  double precision pval
  double precision x

  pval = p(n)
  do i = n - 1, 0, -1
    pval = pval * x + p(i)
  end do

  return
end
subroutine dpoly2_roots ( p, r )

!*****************************************************************************80
!
!! DPOLY2_ROOTS finds the roots of a quadratic polynomial.
!
!  Discussion:
!
!    The standard quadratic formula is used.
!
!  Modified:
!
!    09 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision PCOF(0:2), the polynomial coefficients.
!    P(I) is the coefficient of X**I.
!
!    Output, complex R(2), the roots of the polynomial.
!
  implicit none

  double precision disc
  double precision p(0:2)
  complex ( kind = 16 ) r(2)

  if ( p(2) == 0.0D+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'DPOLY2_ROOTS - Fatal error!'
    write ( *, * ) '  Quadratic coefficient is zero.'
    stop
  end if

  disc = p(1)**2 - 4.0D+00 * p(2) * p(0)

  if ( disc >= 0.0D+00 ) then

    r(1) = cmplx ( 0.5D+00 * ( - p(1) + sqrt ( disc ) ) / p(2), 0.0D+00 )
    r(2) = cmplx ( 0.5D+00 * ( - p(1) - sqrt ( disc ) ) / p(2), 0.0D+00 )

  else if ( disc < 0.0D+00 ) then

    r(1) = cmplx ( - 0.5D+00 * p(1) / p(2), 0.5D+00 * sqrt ( - disc ) / p(2) )
    r(2) = cmplx ( - 0.5D+00 * p(1) / p(2), - 0.5D+00 * sqrt ( - disc ) / p(2) )
 
  end if

  return
end
subroutine dpoly2_roots2 ( p, r, ierror )

!*****************************************************************************80
!
!! DPOLY2_ROOTS2 finds the roots of a quadratic polynomial.
!
!  Discussion:
!
!    An alternate form of the quadratic formula is used.
!
!  Modified:
!
!    09 December 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision PCOF(0:2), the polynomial coefficients.
!    P(I) is the coefficient of X**I.
!
!    Output, double complex R(2), the roots of the polynomial.
!
!    Output, integer IERROR, error flag.
!    0, no error;
!    1, an error occurred.
!
  implicit none

  double precision disc
  integer ierror
  double precision p(0:2)
  complex ( kind = 16 ) r(2)

  ierror = 0

  if ( p(2) == 0.0D+00 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DPOLY2_ROOTS2 - Fatal error!'
    write ( *, '(a)' ) '  Quadratic coefficient is zero.'
    stop
  end if

  disc = p(1)**2 - 4.0D+00 * p(2) * p(0)

  if ( disc >= 0.0D+00 ) then

    if ( - p(1) + sqrt ( disc ) == 0.0D+00 ) then
      ierror = 1
      r(1) = cmplx ( 0.0D+00, 0.0D+00 )
    else
      r(1) = cmplx ( 2.0D+00 * p(0) / ( - p(1) + sqrt ( disc ) ), 0.0D+00 )
    end if

    if ( - p(1) - sqrt ( disc ) == 0.0D+00 ) then
      ierror = 1
      r(2) = cmplx ( 0.0D+00, 0.0D+00 )
    else
      r(2) = cmplx ( 2.0D+00 * p(0) / ( - p(1) - sqrt ( disc ) ), 0.0D+00 )
    end if
!
!  Need to revise this part of the calculation.
!
  else if ( disc < 0.0D+00 ) then

    r(1) = cmplx ( - 0.5D+00 * p(1) / p(2), + 0.5D+00 * sqrt ( - disc ) / p(2) )
    r(2) = cmplx ( - 0.5D+00 * p(1) / p(2), - 0.5D+00 * sqrt ( - disc ) / p(2) )

  end if

  return
end
function fmin ( ax, bx, f, tol )

!*****************************************************************************80
!
!! FMIN seeks a minimizer of a scalar function of a scalar variable.
!
!  Discussion:
!
!    FMIN seeks an approximation to the point where F attains a minimum on
!    the interval (AX,BX).
!
!    The method used is a combination of golden section search and
!    successive parabolic interpolation.  Convergence is never much
!    slower than that for a Fibonacci search.  If F has a continuous
!    second derivative which is positive at the minimum (which is not
!    at AX or BX), then convergence is superlinear, and usually of the
!    order of about 1.324....
!
!    The function F is never evaluated at two points closer together
!    than EPS * ABS ( FMIN ) + (TOL/3), where EPS is approximately the
!    square root of the relative machine precision.  If F is a unimodal
!    function and the computed values of F are always unimodal when
!    separated by at least EPS * ABS ( XSTAR ) + (TOL/3), then FMIN
!    approximates the abcissa of the global minimum of F on the
!    interval AX, BX with an error less than 3 * EPS * ABS ( FMIN ) + TOL.
!    If F is not unimodal, then FMIN may approximate a local, but
!    perhaps non-global, minimum to the same accuracy.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization without Derivatives,
!    Prentice Hall, 1973.
!
!    David Kahaner, Clever Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1988.
!
!  Parameters
!
!    Input/output, real AX, BX.  On input, the left and right endpoints
!    of the initial interval.  On output, the lower and upper bounds for
!    the minimizer.
!
!    Input, external F, a real function of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates F(X) for any X in the interval (AX,BX).
!
!    Input, real TOL, the desired length of the interval of uncertainty of the
!    final result ( >= 0.0)
!
!    Output, real FMIN, the abcissa approximating the minimizer of f.
!
  implicit none

  real a
  real ax
  real b
  real bx
  real c
  real d
  real e
  real eps
  real, external :: f
  real fmin
  real fu
  real fv
  real fw
  real fx
  real p
  real q
  real r
  real tol
  real tol1
  real tol2
  real u
  real v
  real w
  real x
  real xm

  c = 0.5E+00 * ( 3.0E+00 - sqrt ( 5.0E+00 ) )
!
!  C is the squared inverse of the golden ratio.
!
!  EPS is the square root of the relative machine precision.
!
  eps = sqrt ( epsilon ( eps ) )
!
!  Initialization.
!
  a = ax
  b = bx
  v = a + c * ( b - a )
  w = v
  x = v
  e = 0.0E+00
  fx = f(x)
  fv = fx
  fw = fx
!
!  Main loop starts here.
!
20 continue

  xm = 0.5E+00 * ( a + b )
  tol1 = eps * abs ( x ) + tol / 3.0E+00
  tol2 = 2.0E+00 * tol1
!
!  Check the stopping criterion.
!
  if ( abs ( x - xm ) <= ( tol2 - 0.5E+00 * ( b - a ) ) ) then
    fmin = x
    return
  end if
!
!  Is golden-section necessary?
!
  if ( abs ( e ) <= tol1 ) then
    go to 40
  end if
!
!  Fit a parabola.
!
  r = ( x - w ) * ( fx - fv )
  q = ( x - v ) * ( fx - fw )
  p = ( x - v ) * q - ( x - w ) * r
  q = 2.0E+00 * ( q - r )
  if ( q > 0.0E+00 ) then
    p = -p
  end if
  q = abs ( q )
  r = e
  e = d
!
!  Is a parabola acceptable?
!
   30 continue

  if ( abs ( p ) >= abs ( 0.5E+00 * q * r ) ) then
    go to 40
  end if

  if ( p <= q * ( a - x ) ) then
    go to 40
  end if

  if ( p >= q * ( b - x ) ) then
    go to 40
  end if
!
!  A parabolic interpolation step
!
  d = p / q
  u = x + d
!
!  F must not be evaluated too close to AX or BX.
!
  if ( ( u - a ) < tol2 ) then
    d = sign ( tol1, xm - x )
  end if

  if ( ( b - u ) < tol2 ) then
    d = sign ( tol1, xm - x )
  end if

  go to 50
!
!  A golden-section step.
!
   40 continue

  if ( x >= xm ) then
    e = a - x
  else
    e = b - x
  end if

  d = c * e
!
!  F must not be evaluated too close to X.
!
   50 continue

  if ( abs ( d ) >= tol1 ) then
    u = x + d
  end if

  if ( abs ( d ) < tol1 ) then
    u = x + sign ( tol1, d )
  end if

  fu = f(u)
!
!  Update  a, b, v, w, and x
!
  if ( fu <= fx ) then

    if ( u >= x ) then
      a = x
    else
      b = x
    end if

    v = w
    fv = fw
    w = x
    fw = fx
    x = u
    fx = fu
    go to 20

  end if

60 continue

  if ( u < x ) then
    a = u
  else
    b = u
  end if

  if ( fu <= fw ) then
    go to 70
  end if

  if ( w == x ) then
    go to 70
  end if

  if ( fu <= fv ) then
    go to 80
  end if

  if ( v == x ) then
    go to 80
  end if

  if ( v == w ) then
    go to 80
  end if

  go to 20

   70 continue

  v = w
  fv = fw
  w = u
  fw = fu
  go to 20
   80 continue

  v = u
  fv = fu
  go to 20

end
function isamax ( n, x, incx )

!*****************************************************************************80
!
!! ISAMAX finds the index of the vector element of maximum absolute value.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer ISAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none

  integer i
  integer incx
  integer isamax
  integer ix
  integer n
  real samax
  real x(*)

  if ( n <= 0 ) then

    isamax = 0

  else if ( n == 1 ) then

    isamax = 1

  else if ( incx == 1 ) then

    isamax = 1
    samax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > samax ) then
        isamax = i
        samax = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    isamax = 1
    samax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        isamax = i
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function lcm_12n ( n )

!*****************************************************************************80
!
!! LCM_12N computes the least common multiple of the integers 1 through N.
!
!  Examples:
!
!    N    LCM_12N
!
!    1          1
!    2          2
!    3          3
!    4         12
!    5         60
!    6         60
!    7        420
!    8        840
!    9       2520
!   10       2520
!
!  Modified:
!
!    18 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the value of N.
!
!    Output, integer LCM_12N, the least common multiple of the integers 1 to N.
!
  implicit none

  integer i
  integer imult
  integer j
  integer lcm_12n
  integer n

  lcm_12n = 1

  do i = 2, n

    imult = i

    do j = 1, i-1

      if ( mod ( imult, (i-j) )   ==   0 ) then
        imult = imult / ( i - j )
      end if

    end do

    lcm_12n = lcm_12n * imult

  end do

  return
end
subroutine matrix_exponential_taylor ( n, a, a_exp )

!*****************************************************************************80
!
!! MAXTRIX_EXPONENTIAL_TAYLOR uses a Taylor series for the matrix exponential.
!