!  Discussion:
!
!    This is a simple method, but NOT recommended.   It is easy to
!    find examples for which this method fails.
!
!    Successive terms in the Taylor series are added.
!
!  Modified:
!
!    04 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Cleve Moler and Charles Van Loan,
!    19 Dubious Ways to Compute the Exponential of a Matrix, 25 Years Later,
!    SIAM Review, 
!    Volume 45, Number 1, pages 3-49, March 2003.
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input, real A(N,N), the matrix whose exponential is desired.
!
!    Output, real A_EXP(N,N), a Taylor estimate for the matrix exponential.
!
  implicit none

  integer n

  real a(n,n)
  real a_exp(n,n)
  real b_exp(n,n)
  real a_k(n,n)
  real fact_k
  integer i
  integer k
  real, parameter :: tol = 0.0E+00

  a_exp(1:n,1:n) = 0.0E+00

  do i = 1, n
    a_exp(i,i) = 1.0E+00
  end do

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

  k = 1

  do 

    b_exp(1:n,1:n) = a_exp(1:n,1:n)

    a_exp(1:n,1:n) = a_exp(1:n,1:n) + a_k(1:n,1:n)

    b_exp(1:n,1:n) = abs ( b_exp(1:n,1:n) - a_exp(1:n,1:n) )

    if ( all ( b_exp(1:n,1:n) <= tol ) ) then
      exit
    end if

    k = k + 1
 
    a_k(1:n,1:n) = matmul ( a_k(1:n,1:n), a(1:n,1:n) ) / real ( k )
    
  end do

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

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

  integer n

  integer i
  real p(0:n)
  real pval
  real x

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

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

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

  integer n

  integer i
  real p(0:n)
  real pval
  real x

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

  return
end
subroutine rpoly2_roots ( p, r )

!*****************************************************************************80
!
!! RPOLY2_ROOTS finds the roots of a quadratic polynomial.
!
!  Discussion:
!
!    The standard quadratic formula is used.
!
!  Modified:
!
!    09 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real 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

  real disc
  real p(0:2)
  complex r(2)

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

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

  if ( disc >= 0.0E+00 ) then

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

  else if ( disc < 0.0E+00 ) then

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

  return
end
subroutine rpoly2_roots2 ( p, r, ierror )

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

  real disc
  integer ierror
  real p(0:2)
  complex r(2)

  ierror = 0

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

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

  if ( disc >= 0.0E+00 ) then

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

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

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

  end if

  return
end
function samax ( n, x, incx )

!*****************************************************************************80
!
!! SAMAX returns the maximum absolute value of the entries in a vector.
!
!  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 X.
!
!    Output, real SAMAX, the maximum absolute value of an element of X.
!
  implicit none

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

  if ( n <= 0 ) then

    samax = 0.0E+00

  else if ( n == 1 ) then

    samax = abs ( x(1) )

  else if ( incx == 1 ) then

    samax = abs ( x(1) )

    do i = 2, n
      if ( abs ( x(i) ) > samax ) then
        samax = abs ( x(i) )
      end if
    end do

  else

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

    samax = abs ( x(ix) )
    ix = ix + incx

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

  end if

  return
end
subroutine saxpy ( n, sa, x, incx, y, incy )

!*****************************************************************************80
!
!! SAXPY adds a constant times one vector to another.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), the vector to which a multiple of X is to
!    be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)

  if ( n <= 0 ) then

  else if ( sa == 0.0E+00 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = y(1:n) + sa * x(1:n)

  else

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

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine scopy ( n, x, incx, y, incy )

!*****************************************************************************80
!
!! SCOPY copies one real vector into another.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be copied into Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real Y(*), the copy of X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none

  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real x(*)
  real y(*)

  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = x(1:n)

  else

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

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
function sdot ( n, x, incx, y, incy )

!*****************************************************************************80
!
!! SDOT forms the dot product of two vectors.
!
!  Modified:
!
!    02 June 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real SDOT, the dot product of X and Y.
!
  implicit none

  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sdot
  real stemp
  real x(*)
  real y(*)

  if ( n <= 0 ) then

    sdot = 0.0E+00

  else if ( incx == 1 .and. incy == 1 ) then

    sdot = dot_product ( x(1:n), y(1:n) )

  else

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

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    stemp = 0.0E+00
    do i = 1, n
      stemp = stemp + x(ix) * y(iy)
      ix = ix + incx
      iy = iy + incy
    end do

    sdot = stemp

  end if

  return
end
function sdsdot ( n, x, incx, y, incy )

!*****************************************************************************80
!
!! SDSDOT forms the dot product of two vectors using higher precision.
!
!  Modified:
!
!    02 June 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real SDSDOT, the dot product of X and Y.
!
  implicit none

  double precision dsdot
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sdsdot
  real x(*)
  real y(*)

  if ( n <= 0 ) then

    dsdot = 0.0D+00

  else if ( incx == 1 .and. incy == 1 ) then

    dsdot = dot_product ( dble ( x(1:n) ), dble ( y(1:n) ) )

  else

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

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    dsdot = 0.0D+00
    do i = 1, n