#! /usr/bin/env python
#
def c8mat_expm1 ( n, a ):

#*****************************************************************************80
#
## C8MAT_EXPM1 is essentially MATLAB's built-in matrix exponential algorithm.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    05 March 2013
#
#  Author:
#
#    Python version by John Burkardt.
#
#  Reference:
#
#    Cleve Moler, Charles VanLoan,
#    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
#    Twenty-Five Years Later,
#    SIAM Review,
#    Volume 45, Number 1, March 2003, pages 3-49.
#
#  Parameters:
#
#    Input, int N, the dimension of the matrix.
#
#    Input, double complex A[N*N], the matrix.
#
#    Output, double complex C8MAT_EXPM1[N*N], the estimate for exp(A).
#
  import numpy as np

  q = 6

  a2 = a.copy ( )

  a_norm = np.linalg.norm ( a2, np.inf )

  ee = ( int ) ( np.log2 ( a_norm ) ) + 1
  
  s = max ( 0, ee + 1 )

  a2 = a2 / ( 2.0 ** s )

  x = a2.copy ( )

  c = 0.5

  e = np.eye ( n, dtype = np.complex64 ) + c * a2

  d = np.eye ( n, dtype = np.complex64 ) - c * a2

  p = True

  for k in range ( 2, q + 1 ):
  
    c = c * float ( q - k + 1 ) / float ( k * ( 2 * q - k + 1 ) )

    x = np.dot ( a2, x )

    e = e + c * x

    if ( p ):
      d = d + c * x
    else:
      d = d - c * x

    p = not p
#
#  E -> inverse(D) * E
#
  e = np.linalg.solve ( d, e )
#
#  E -> E^(2*S)
#
  for k in range ( 0, s ):
    e = np.dot ( e, e )

  return e

def r8mat_expm1 ( n, a ):

#******************************************************************************/
#
## R8MAT_EXPM1 is essentially MATLAB's built-in matrix exponential algorithm.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    05 February 2017
#
#  Author:
#
#    Cleve Moler, Charles Van Loan
#
#  Reference:
#
#    Cleve Moler, Charles VanLoan,
#    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
#    Twenty-Five Years Later,
#    SIAM Review,
#    Volume 45, Number 1, March 2003, pages 3-49.
#
#  Parameters:
#
#    Input, int N, the dimension of the matrix.
#
#    Input, double A[N,N], the matrix.
#
#    Output, double E[N,N], the estimate for exp(A).
#
  import numpy as np

  q = 6
 
  a2 = a.copy ( )

  a_norm = np.linalg.norm ( a2, np.inf )

  ee = ( int ) ( np.log2 ( a_norm ) ) + 1
  
  s = max ( 0, ee + 1 )

  a2 = a2 / ( 2.0 ** 2 )

  x = a2.copy ( )

  c = 0.5

  e = np.eye ( n ) + c * a2

  d = np.eye ( n ) - c * a2

  p = True

  for k in range ( 2, q + 1 ):

    c = c * float ( q - k + 1 ) / float ( k * ( 2 * q - k + 1 ) )

    x = np.dot ( a2, x )

    e = e + c * x

    if ( p ):
      d = d + c * x
    else:
      d = d - c * x

    p = not p
#
#  E -> inverse(D) * E
#
  e = np.linalg.solve ( d, e )
#
#  E -> E^(2*S)
#
  for k in range ( 0, 2 ):
    e = np.dot ( e, e )

  return e

def r8mat_expm2 ( n, a ):

#******************************************************************************/
#
#  Purpose:
#
#    R8MAT_EXPM2 uses the Taylor series for the matrix exponential.
#
#  Discussion:
#
#    Formally,
#
#      exp ( A ) = I + A + 1/2 A^2 + 1/3! A^3 + ...
#
#    This function sums the series until a tolerance is satisfied.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    05 February 2017
#
#  Author:
#
#    Cleve Moler, Charles Van Loan
#
#  Reference:
#
#    Cleve Moler, Charles VanLoan,
#    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
#    Twenty-Five Years Later,
#    SIAM Review,
#    Volume 45, Number 1, March 2003, pages 3-49.
#
#  Parameters:
#
#    Input, int N, the dimension of the matrix.
#
#    Input, double A[N,N], the matrix.
#
#    Output, double E[N,N], the estimate for exp(A).
#
  import numpy as np

  e = np.zeros ( [ n, n ] )

  f = np.eye ( n )

  k = 1

  while ( r8mat_is_significant ( n, n, e, f ) ):
  
    e = e + f

    f = np.dot ( a, f )

    f = f / float ( k )

    k = k + 1

  return e

def r8mat_expm3 ( n, a ):

#******************************************************************************/
#
## R8MAT_EXPM3 approximates the matrix exponential using an eigenvalue approach.
#
#  Discussion:
#
#    exp(A) = V * D * inv(V)
#
#    where V is the matrix of eigenvectors of A, and D is the diagonal matrix
#    whose i-th diagonal entry is exp(lambda(i)), for lambda(i) an eigenvalue
#    of A.
#
#    This function is accurate for matrices which are symmetric, orthogonal,
#    or normal.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    05 February 2017
#
#  Author:
#
#    Cleve Moler, Charles Van Loan
#
#  Reference:
#
#    Cleve Moler, Charles VanLoan,
#    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
#    Twenty-Five Years Later,
#    SIAM Review,
#    Volume 45, Number 1, March 2003, pages 3-49.
#
#  Parameters:
#
#    Input, int N, the dimension of the matrix.
#
#    Input, double A[N,N], the matrix.
#
#    Output, double E[N,N], the estimate for exp(A).
#
  import numpy as np

  cevals, cevecs = np.linalg.eig ( a )
#
#  Need to take the real part of these quantities!
#
  evals = cevals.real
  evecs = cevecs.real

  exp_evals = np.exp ( evals )
  d2 = np.diag ( exp_evals )
#
#  Pardon this godawful circumlocution.
#
  b = np.dot ( evecs, d2 )
  bt = b.transpose ( )
  a = evecs
  at = a.transpose ( )
  et, residuals, rank, s = np.linalg.lstsq ( at, bt, rcond = None )
  e = et.transpose ( )

  return e

def r8mat_is_significant ( m, n, r, s ):

#*****************************************************************************80
#
## R8MAT_IS_SIGNIFICANT determines if an R8MAT is relatively significant.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 February 217
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer M, N, the dimension of the matrices.
#
#    Input, real R(M,N), the vector to be compared against.
#
#    Input, real S(M,N), the vector to be compared.
#
#    Output, logical R8MAT_IS_SIGNIFICANT, is true if S is significant
#    relative to R.
#
  eps = 2.220446049250313E-016

  value = False

  for j in range ( 0, n ):
    for i in range ( 0, m ):

      t = r[i,j] + s[i,j]
      tol = eps * abs ( r[i,j] )

      if ( tol < abs ( r[i,j] - t ) ):
        value = True
        break
  
  return value

def matrix_exponential_test ( ):

#******************************************************************************/
#
## MATRIX_EXPONENTIAL_TEST tests MATRIX_EXPONENTIAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    05 February 2017
#
#  Author:
#
#    John Burkardt
#
  import platform
  import sys

  sys.path.append ( '/home/jburkardt/public_html/py_src/test_matrix_exponential' ) 

  print ( '' )
  print ( 'MATRIX_EXPONENTIAL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test the MATRIX_EXPONENTIAL library.' )
  print ( '  This test needs the TEST_MATRIX_EXPONENTIAL library.' )

  matrix_exponential_test01 ( )
  matrix_exponential_test02 ( )

#
#  Terminate.
#
  print ( '' )
  print ( 'MATRIX_EXPONENTIAL_TEST:' )
  print ( '  Normal end of execution.' )
  return 

def matrix_exponential_test01 ( ):

#******************************************************************************/
#
## MATRIX_EXPONENTIAL_TEST01 compares real matrix exponential algorithms.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 February 2017
#
#  Author:
#
#    John Burkardt
#
  from test_matrix_exponential import r8mat_exp_a
  from test_matrix_exponential import r8mat_exp_expa
  from test_matrix_exponential import r8mat_exp_n
  from test_matrix_exponential import r8mat_exp_story
  from test_matrix_exponential import r8mat_exp_test_num
  from test_matrix_exponential import r8mat_print

  print ( '' )
  print ( 'MATRIX_EXPONENTIAL_TEST01:' )
  print ( '  R8MAT_EXPM1 is an equivalent to EXPM' )
  print ( '  R8MAT_EXPM2 uses a Taylor series approach' )
  print ( '  R8MAT_EXPM3 relies on an eigenvalue calculation.' )

  test_num = r8mat_exp_test_num ( )

  for test in range ( 1, test_num + 1 ):

    print ( '' )
    print ( '  Test #%d' % (test ) )

    r8mat_exp_story ( test )

    n = r8mat_exp_n ( test )

    print ( '  Matrix order N = %d' % ( n ) )

    a = r8mat_exp_a ( test, n )

    r8mat_print ( n, n, a, '  Matrix:' )

    a_exp = r8mat_expm1 ( n, a )
    r8mat_print ( n, n, a_exp, '  R8MAT_EXPM1(A):' )

    a_exp = r8mat_expm2 ( n, a )
    r8mat_print ( n, n, a_exp, "  R8MAT_EXPM2(A):" )

    a_exp = r8mat_expm3 ( n, a )
    r8mat_print ( n, n, a_exp, "  R8MAT_EXPM3(A):" )

    a_exp = r8mat_exp_expa ( test, n )
    r8mat_print ( n, n, a_exp, "  Exact Exponential:" )

  return

def matrix_exponential_test02 ( ):

#*****************************************************************************80
#
## MATRIX_EXPONENTIAL_TEST02 compares complex matrix exponential algorithms.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 February 2017
#
#  Author:
#
#    John Burkardt
#
  from test_matrix_exponential import c8mat_exp_a
  from test_matrix_exponential import c8mat_exp_expa
  from test_matrix_exponential import c8mat_exp_n
  from test_matrix_exponential import c8mat_exp_story
  from test_matrix_exponential import c8mat_exp_test_num
  from test_matrix_exponential import c8mat_print

  print ( '' )
  print ( 'MATRIX_EXPONENTIAL_TEST02:' )
  print ( '  EXPM is MATLAB\'s matrix exponential function' )
  print ( '  C8MAT_EXPM1 is an equivalent to EXPM' )
  print ( '  C8MAT_EXPM2 uses a Taylor series approach' )
  print ( '  C8MAT_EXPM3 relies on an eigenvalue calculation.' )

  test_num = c8mat_exp_test_num ( )

  for test in range ( 1, test_num + 1 ): 
  
    print ( '' )
    print ( '  Test #%d' % ( test ) )

    c8mat_exp_story ( test )

    n = c8mat_exp_n ( test )

    print ( '  Matrix order N = %d' % ( n ) )

    a = c8mat_exp_a ( test, n )
    c8mat_print ( n, n, a, '  Matrix:' )

    a_exp = c8mat_expm1 ( n, a )
    c8mat_print ( n, n, a_exp, '  C8MAT_EXPM1(A):' )
#
#   a_exp = c8mat_expm2 ( n, a )
#   c8mat_print ( n, n, a_exp, '  C8MAT_EXPM2(A):' )
#
#   a_exp = c8mat_expm3 ( n, a )
#   c8mat_print ( n, n, a_exp, '  C8MAT_EXPM3(A):' )
#
    a_exp = c8mat_exp_expa ( test, n )
    c8mat_print ( n, n, a_exp, '  Exact Exponential:' )

  return

def timestamp ( ):

#*****************************************************************************80
#
## TIMESTAMP prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

def timestamp_test ( ):

#*****************************************************************************80
#
## TIMESTAMP_TEST tests TIMESTAMP.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    03 December 2014
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import platform

  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  TIMESTAMP prints a timestamp of the current date and time.' )
  print ( '' )

  timestamp ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  timestamp ( )
  matrix_exponential_test ( )
  timestamp ( )