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

#*****************************************************************************80
#
## CONDITION_LINPACK estimates the L1 condition number of a matrix.
#
#  Discussion:
#
#    The R8GE storage format is used for a general M by N matrix.  A storage 
#    space is made for each logical entry.  The two dimensional logical
#    array is mapped to a vector, in which storage is by columns.
#
#    For the system A * X = B, relative perturbations in A and B
#    of size EPSILON may cause relative perturbations in X of size
#    EPSILON*RCOND.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 July 2015
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Bunch, Moler, Stewart,
#    LINPACK User's Guide,
#    SIAM, 1979
#
#  Parameters:
#
#    Input, integer N, the order of the matrix A.
#
#    Input, real A(N,N), a matrix to be factored.
#
#    Output, real VALUE, an estimate of the condition number of A.
#
  import numpy as np
  from r8_sign import r8_sign
  from r8ge import r8ge_fa
  from r8mat_norm_l1 import r8mat_norm_l1
#
#  Compute the L1 norm of A.
#
  anorm = r8mat_norm_l1 ( n, n, a )
#
#  Compute the LU factorization.
#
  a_lu, pivot, info = r8ge_fa ( n, a )
#
#  COND = norm(A) * (estimate of norm(inverse(A))) 
#
#  estimate of norm(inverse(A)) = norm(Z) / norm(Y)
#
#  where
#    A * Z = Y
#  and
#    A' * Y = E
#
#  The components of E are chosen to cause maximum local growth in the
#  elements of W, where U'*W = E.  The vectors are frequently rescaled
#  to avoid overflow.
#
#  Solve U' * W = E.
#
  ek = 1.0
  z = np.zeros ( n )

  for k in range ( 0, n ):

    if ( z[k] != 0.0 ):
      ek = - r8_sign ( z[k] ) * abs ( ek )

    if ( abs ( a_lu[k,k] ) < abs ( ek - z[k] ) ):
      s = abs ( a_lu[k,k] ) / abs ( ek - z[k] )
      for i in range ( 0, n ):
        z[i] = s * z[i]
      ek = s * ek

    wk = ek - z[k]
    wkm = - ek - z[k]
    s = abs ( wk )
    sm = abs ( wkm )

    if ( a_lu[k,k] != 0.0 ):
      wk = wk / a_lu[k,k]
      wkm = wkm / a_lu[k,k]
    else:
      wk = 1.0
      wkm = 1.0

    if ( k + 1 <= n - 1 ):

      for j in range ( k + 1, n ):
        sm = sm + abs ( z[j] + wkm * a_lu[k,j] )
        z[j] = z[j] + wk * a_lu[k,j]
        s = s + abs ( z[j] )

      if ( s < sm ):
        t = wkm - wk
        wk = wkm
        for j in range ( k + 1, n ):
          z[j] = z[j] + t * a_lu[k,j]

    z[k] = wk

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

  for i in range ( 0, n ):
    z[i] = z[i] / t
#
#  Solve L' * Y = W
#
  for k in range ( n - 1, -1, -1 ):

    for i in range ( k + 1, n ):
      z[k] = z[k] + z[i] * a_lu[i,k]

    t = abs ( z[k] )

    if ( 1.0 < t ):
      for i in range ( 0, n ):
        z[i] = z[i] / t

    l = pivot[k]

    t = z[l]
    z[l] = z[k]
    z[k] = t

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

  for i in range ( 0, n ):
    z[i] = z[i] / t

  ynorm = 1.0
#
#  Solve L * V = Y.
#
  for k in range ( 0, n ):

    l = pivot[k]

    t = z[l]
    z[l] = z[k]
    z[k] = t

    for i in range ( k + 1, n ):
      z[i] = z[i] + t * a_lu[i,k]

    if ( 1.0 < abs ( z[k] ) ):
      t = abs ( z[k] )
      ynorm = ynorm / t
      for i in range ( 0, n ):
        z[i] = z[i] / t

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

  for i in range ( 0, n ):
    z[i] = z[i] / t

  ynorm = ynorm / t
#
#  Solve U * Z = V.
#
  for k in range ( n - 1, -1, -1 ):

    if ( abs ( a_lu[k,k] ) < abs ( z[k] ) ):
      s = abs ( a_lu[k,k] ) / abs ( z[k] )
      for i in range ( 0, n ):
        z[i] = s * z[i]
      ynorm = s * ynorm

    if ( a_lu[k,k] != 0.0 ):
      z[k] = z[k] / a_lu[k,k]
    else:
      z[k] = 1.0

    for i in range ( 0, k ):
      z[i] = z[i] - a_lu[i,k] * z[k]
#
#  Normalize Z in the L1 norm.
#
  t = 0.0
  for i in range ( 0, n ):
    t = t + abs ( z[i] )

  for i in range ( 0, n ):
    z[i] = z[i] / t

  ynorm = ynorm / t

  value = anorm / ynorm

  return value

def condition_linpack_test ( ):

#*****************************************************************************80
#
## CONDITION_LINPACK_TEST tests CONDITION_LINPACK.
#
#  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_LINPACK_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  CONDITION_LINPACK estimates the L1 condition number' )
  print ( '  for a matrix in general storage.' )
  print ( '' )
  print ( '  Matrix               Order   Condition       Linpack' )
  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_linpack ( 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_linpack ( 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_linpack ( 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_linpack ( 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_linpack ( 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_linpack ( 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_linpack ( n, a )
    print ( '  %20s  %4d  %14.6g  %14.6g' % ( name, n, cond_l1, cond ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'CONDITION_LINPACK_TEST' )
  print ( '  Normal end of execution.' )
  return

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