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

#*****************************************************************************80
#
## CONDITION_HAGER estimates the L1 condition number of a matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    William Hager,
#    Condition Estimates,
#    SIAM Journal on Scientific and Statistical Computing,
#    Volume 5, Number 2, June 1984, pages 311-316.
#
#  Parameters:
#
#    Input, integer N, the dimension of the matrix.
#
#    Input, real A(N,N), the matrix.
#
#    Output, real VALUE, an estimate of the L1 condition number.
#
  import numpy as np
  from r8_sign import r8_sign
  from r8mat_norm_l1 import r8mat_norm_l1
  from r8vec_max_abs_index import r8vec_max_abs_index

  i1 = -1
  c1 = 0.0
#
#  Factor the matrix.
#
  b = np.zeros ( n )
  for i in range ( 0, n ):
    b[i] = 1.0 / float ( n )

  while ( True ):

    b2 = np.linalg.solve ( a, b )

    for i in range ( 0, n ):
      b[i] = b2[i]

    c2 = 0.0
    for i in range ( 0, n ):
      c2 = c2 + abs ( b[i] )

    for i in range ( 0, n ):
      b[i] = r8_sign ( b[i] )

    b2 = np.linalg.solve ( np.transpose ( a ), b )

    for i in range ( 0, n ):
      b[i] = b2[i]

    i2 = r8vec_max_abs_index ( n, b )

    if ( 0 <= i1 ):
      if ( i1 == i2 or c2 <= c1 ):
        break

    i1 = i2
    c1 = c2

    for i in range ( 0, n ):
      b[i] = 0.0
    b[i1] = 1.0

  value = c2 * r8mat_norm_l1 ( n, n, a )

  return value

def condition_hager_test ( ):

#*****************************************************************************80
#
## CONDITION_HAGER_TEST tests CONDITION_HAGER.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from combin import combin
  from combin import combin_inverse
  from conex1 import conex1
  from conex1 import conex1_inverse
  from conex2 import conex2
  from conex2 import conex2_inverse
  from conex3 import conex3
  from conex3 import conex3_inverse
  from conex4 import conex4
  from conex4 import conex4_inverse
  from kahan import kahan
  from kahan import kahan_inverse
  from r8mat_norm_l1 import r8mat_norm_l1
  from r8mat_uniform_01 import r8mat_uniform_01

  print ( '' )
  print ( 'CONDITION_HAGER_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  CONDITION_HAGER estimates the L1 condition number' )
  print ( '  for a matrix in general storage.' )
  print ( '' )
  print ( '  Matrix               Order   Condition         Hager' )
  print ( '' )
#
#  Combinatorial matrix.
#
  name = 'Combinatorial'
  n = 4
  alpha = 2.0
  beta = 3.0
  a = combin ( alpha, beta, n )
  a_inverse = combin_inverse ( alpha, beta, n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  CONEX1
#
  name = 'CONEX1'
  n = 4
  alpha = 100.0
  a = conex1 ( alpha )
  a_inverse = conex1_inverse ( alpha )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  CONEX2
#
  name = 'CONEX2'
  n = 3
  alpha = 100.0
  a = conex2 ( alpha )
  a_inverse = conex2_inverse ( alpha )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  CONEX3
#
  name = 'CONEX3'
  n = 5
  a = conex3 ( n )
  a_inverse = conex3_inverse ( n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  CONEX4
#
  name = 'CONEX4'
  n = 4
  a = conex4 ( )
  a_inverse = conex4_inverse ( )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  KAHAN
#
  name = 'KAHAN'
  n = 4
  alpha = 0.25
  a = kahan ( alpha, n, n )
  a_inverse = kahan_inverse ( alpha, n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1  = a_norm_l1 * a_inverse_norm_l1
  cond = condition_hager ( n, a )
  print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  Random
#
  seed = 123456789

  for j in range ( 0, 5 ):
    name = 'RANDOM'
    n = 4
    a, seed = r8mat_uniform_01 ( n, n, seed )
    a_inverse = np.linalg.inv ( a )
    a_norm_l1 = r8mat_norm_l1 ( n, n, a )
    a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
    cond_l1  = a_norm_l1 * a_inverse_norm_l1
    cond = condition_hager ( n, a )
    print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'CONDITION_HAGER_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  condition_hager_test ( )
  timestamp ( )