#! /usr/bin/env python3
#
def linpack_d_test ( ):

#*****************************************************************************80
#
## linpack_d_test() tests linpack_d().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  from numpy.random import default_rng
  import numpy as np
  import platform

  rng = default_rng ( )

  print ( '' )
  print ( 'linpack_d_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  Test linpack_d().' )

  dgeco_test ( )
  dgedet_test ( )
  dgefa_test ( )
  dgeinv_test ( )
  dgesl_test ( )
  dgeslt_test ( )

  dpofa_test ( )

  dqrdc_test ( )
  dqrsl_test ( )

  drotg_test ( rng )

  dsvdc_test ( rng )
#
#  Terminate.
#
  print ( '' )
  print ( 'linpack_d_test():' )
  print ( '  Normal end of execution.' )
  return

def dgeco ( A, n ):

#*****************************************************************************80
#
## dgeco() factors a matrix and estimates its condition number.
#
#  Discussion:
#
#    The value rcond returned is actually the reciprocal of the condition.
#    rcond = 1/(norm(a) * (estimate of norm(inverse(a)))) .
#
#  Input:
#
#    real A[n,n]: the matrix.
#
#    integer n: the order of the matrix.
#
#  Output:
#
#    real rcond: the reciprocal condition number.
#
#    real ALU[n,n]: the LU factor information.
#
#    integer ipvt[n]: the pivot vector.
#
  import numpy as np
#
#  Compute the L1-norm of A.
#
  anorm = 0.0
  for j in range ( 0, n ):
    anorm = max ( anorm, np.sum ( np.abs ( A[0:n,j] ) ) )
#
#  Compute the LU factors of A.
#
  ALU, ipvt, info = dgefa ( A, n )
#
#  Solve U'*w = e
#
  ek = 1.0
  z = np.zeros ( n )

  for k in range ( 0, n ):

    if ( 0.0 < z[k] ):
      ek = - np.abs ( ek )
    elif ( z[k] < 0.0 ):
      ek = np.abs ( ek )

    if ( np.abs ( ALU[k,k] ) < np.abs ( ek - z[k] ) ):
      s = np.abs ( ALU[k,k] ) / np.abs ( ek - z[k] )
      z[0:n] = z[0:n] * s
      ek = ek * s

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

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

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

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

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

    z[k] = wk
#
#  Normalize z.
#
  s = np.sum ( np.abs ( z[0:n] ) )
  z[0:n] = z[0:n] / s
#
#  Solve L' * y = w
#
  for k in range ( n - 1, -1, -1 ):
    if ( k < n - 1 ):
      z[k] = z[k] + np.dot ( ALU[k+1:n,k], z[k+1:n] )
    if ( 1.0 < np.abs ( z[k] ) ):
      z[0:n] = z[0:n] / np.abs ( z[k] )

    l = ipvt[k]
    t = z[l]
    z[l] = z[k]
    z[k] = t
#
#  Normalize z.
#
  s = np.sum ( np.abs ( z[0:n] ) )
  z[0:n] = z[0:n] / s

  ynorm = 1.0
#
#  Solve L * v = y
#
  for k in range ( 0, n ):
    l = ipvt[k]
    t =  z[l]
    z[l] = z[k]
    z[k] = t

    if ( k < n - 1 ):
      z[k+1:n] = z[k+1:n] + t * ALU[k+1:n,k]
    if ( 1.0 < np.abs ( z[k] ) ):
      s = np.abs( z[k] )
      z[0:n] = s[0:n] / s
      ynorm = ynorm / s

  s = np.sum ( np.abs ( z[0:n] ) )
  z[0:n] = z[0:n] / s
  ynorm = ynorm / s
#
#  Solve  u * z = v
#
  for k in range ( n - 1, -1, -1 ):
     if ( np.abs ( ALU[k,k] ) < np.abs ( z[k] ) ):
       s = np.abs ( ALU[k,k] ) / np.abs ( z[k] )
       z[0:n] = z[0:n] * s
       ynorm = s * ynorm
     if ( ALU[k,k] != 0.0 ):
       z[k] = z[k] / ALU[k,k]
     else:
       z[k] = 1.0
     t = - z[k]
     z[0:k] = z[0:k] + t * ALU[0:k,k]
#
#  Normalize.
#
  s = np.sum ( np.abs ( z[0:n] ) )
  z[0:n] = z[0:n] / s 
  ynorm = ynorm / s

  if ( anorm != 0.0 ):
    rcond = ynorm / anorm
  else:
    rcond = 0.0

  return rcond, ALU, ipvt

def dgeco_test ( ):

#*****************************************************************************80
#
## dgeco_test() tests dgeco().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgeco_test()' )
  print ( '  dgeco() computes the condition number of a matrix stored in' )
  print ( '  general format.' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  dgeco() factors the matrix.
#
  print ( '' )
  print ( '  dgeco() factors and analyzes the matrix' )

  rcond, ALU, ipvt = dgeco ( A, n )

  print ( '' )
  print ( '  The reciprocal condition number RCOND:' )
  print ( rcond )

  return

def dgedet ( ALU, n, ipvt ):

#*****************************************************************************80
#
## dgedet() computes the determinant of a matrix factored by dgeco() or dgefa().
#
#  Discussion:
#
#    dgedet() is a variation of the linpack() function dgedi().
#    Unlike dgedi(), it only computes the determinant of the matrix,
#    but not the inverse.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 June 2024
#
#  Author:
#
#    This version by John Burkardt
#
#  Reference:
#
#    Dongarra, Moler, Bunch and Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Input:
#
#    real ALU[n,n]: an LU factorization computed by dgeco() or dgefa().
#
#    integer n: the order of the matrix.
#
#    integer ipvt[n]: the pivot vector from dgeco() or dgefa().
#
#  Output:
#
#    real det: the determinant of the matrix.
#
  det = 1.0

  for i in range ( 0, n ):

    if ( ipvt[i] != i ):
      det = - det

    det = det * ALU[i,i]

    if ( det == 0.0 ):
      break

  return det

def dgedet_test ( ):

#*****************************************************************************80
#
## dgedet_test() tests dgedet().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgedet_test()' )
  print ( '  dgedet() computes the determinant of a matrix stored in general' )
  print ( '  format, after it has been factored by dgefa().' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  dgefa() factors the matrix.
#
  print ( '' )
  print ( '  dgefa() factors the matrix' )

  ALU, ipvt, info = dgefa ( A, n )

  if ( info != 0 ):
    print ( '  dgefa() returned an error flag INFO = ', info )
    return
#
#  dgedet() computes the determinant.
#
  print ( '' )
  print ( '  dgedet() computes the determinant:' )

  det = dgedet ( ALU, n, ipvt )

  print ( '' )
  print ( '  The determinant DET:' )
  print ( det )

  return

def dgefa ( A, n ):

#*****************************************************************************80
#
## dgefa() factors a real matrix in general storage format.
#
#  Discussion:
#
#    Gaussian elimination with partial pivoting.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 June 2024
#
#  Author:
#
#    This version by John Burkardt
#
#  Reference:
#
#    Dongarra, Moler, Bunch and Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Input:
#
#    real A(N,N), the matrix to be factored.
#
#    integer N, the order of the matrix A.
#
#  Output:
#
#    real ALU(N,N), an upper triangular matrix and the multipliers 
#    used to obtain it.  The factorization can be written A=L*U, where L is 
#    a product of permutation and unit lower triangular matrices, and U is 
#    upper triangular.
#
#    integer IPVT(N), the pivot indices.
#
#    integer INFO, singularity indicator.
#    0, normal value.
#    K, if U(K-1,K-1) == 0.  This is not an error condition for this subroutine,
#    but it does indicate that DGESL or DGEDI will divide by zero if called.
#    Use RCOND in DGECO for a reliable indication of singularity.
#
  import numpy as np

  ALU = A.copy ( )

  ipvt = np.zeros ( n, dtype = int )
  info = 0

  for k in range ( 0, n - 1 ):
#
#  Find L = pivot index,
#  K <= L < N,
#  index of largest entry in column K.
#
    l = -1
    amax = - np.inf
   
    for i in range ( k, n ):
      if ( amax < np.abs ( ALU[i,k] ) ):
        l = i
        amax = np.abs ( ALU[i,k] )

    ipvt[k] = l
#
#  Zero pivot implies this column already triangularized.
#
    if ( ALU[l,k] == 0.0 ):
      info = k + 1
      continue
#
#  Interchange row K and pivot row L if necessary.
#
    if ( l != k ):
      t        = ALU[l,k]
      ALU[l,k] = ALU[k,k]
      ALU[k,k] = t
#
#  Compute multipliers.
#
    ALU[k+1:n,k] = - ALU[k+1:n,k] / ALU[k,k]
#
#  Row elimination with column indexing.
#
    for j in range ( k + 1, n ):

      t = ALU[l,j]
      if ( l != k ):
        ALU[l,j] = ALU[k,j]
        ALU[k,j] = t

      ALU[k+1:n,j] = ALU[k+1:n,j] + t * ALU[k+1:n,k]

  ipvt[n-1] = n - 1

  if ( ALU[n-1,n-1] == 0.0 ):
    info = n

  return ALU, ipvt, info

def dgefa_test ( ):

#*****************************************************************************80
#
## dgefa_test() tests dgefa().
#
#  Discussion:
#
#    Compute the pivot vector and LU factors of a matrix A.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgefa_test()' )
  print ( '  dgefa() computes pivot vector and LU factors of a matrix' )
  print ( '  stored in general format.' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  Factor the matrix.
#
  print ( '' )
  print ( '  dgefa() factors the matrix' )

  ALU, ipvt, info = dgefa ( A, n )

  if ( info != 0 ):
    print ( '  dgefa() returned an error flag INFO = ', info )
    return

  print ( '' )
  print ( '  The matrix ALU:' )
  print ( ALU )

  print ( '' )
  print ( '  The pivot vector IPVT:' )
  print ( ipvt )

  return

def dgeinv ( ALU, n, ipvt ):

#*****************************************************************************80
#
## dgeinv() computes the inverse of a matrix that was factored by dgefa() or dgeco().
#
#  Input:
#
#    real ALU[n,n]: the LU factors from dgeco() or dgefa().
#
#    integer n: the order of the matrix.
#
#    integer ipvt[n]: the pivot vector from dgeco() or dgefa().
#
#  Output:
#
#    real AINV[n,n]: the inverse matrix.
#
  import numpy as np

  AINV = ALU.copy()

  for k in range ( 0, n ):
    AINV[k,k] = 1.0 / AINV[k,k]
    t = - AINV[k,k]
    AINV[0:k,k] = AINV[0:k,k] * t
    for j in range ( k + 1, n ):
      t = AINV[k,j]
      AINV[k,j] = 0.0
      AINV[0:k+1,j] = AINV[0:k+1,j] + t * AINV[0:k+1,k]
#
#  Form inverse(u)*inverse(l)
#
  work = np.zeros ( n )

  for k in range ( n - 1, -1, -1 ):

    for i in range ( k + 1, n ):
      work[i] = AINV[i,k]
      AINV[i,k] = 0.0

    for j in range ( k + 1, n ):
      t = work[j]
      AINV[0:n,k] = AINV[0:n,k] + t * AINV[0:n,j]

    l = ipvt[k]
    if ( l != k ):
      work[0:n] = AINV[0:n,k]
      AINV[0:n,k] = AINV[0:n,l]
      AINV[0:n,l] = work[0:n]

  return AINV

def dgeinv_test ( ):

#*****************************************************************************80
#
## dgeinv_test() tests dgeinv().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgeinv_test()' )
  print ( '  dgeinv() computes inverse of a matrix stored in general' )
  print ( '  format, which has been factored by dgefa().' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  Factor the matrix.
#
  print ( '' )
  print ( '  dgefa() factors the matrix' )

  ALU, ipvt, info = dgefa ( A, n )

  if ( info != 0 ):
    print ( '  dgefa() returned an error flag INFO = ', info )
    return
#
#  Compute the inverse.
#
  print ( '' )
  print ( '  dgeinv() computes the inverse matrix' )

  AINV = dgeinv ( ALU, n, ipvt )

  print ( '' )
  print ( '  The inverse matrix AINV:' )
  print ( AINV )
#
#  Compute A*AINV = I?
#
  B = np.matmul ( A, AINV )

  print ( '' )
  print ( '  The product A * AINV:' )
  print ( B )

  return

def dgesl ( ALU, n, ipvt, b ):

#*****************************************************************************80
#
## dgesl() solves a real general linear system A * X = B.
#
#  Discussion:
#
#    The system matrix must have been factored by DGECO or DGEFA.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 June 2024
#
#  Author:
#
#    This version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Moler, Bunch and Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Input:
#
#    real ALU(N,N), the LU factorization from DGECO or DGEFA.
#
#    integer N, the order of the matrix A.
#
#    integer IPVT(N), the pivot vector from DGECO or DGEFA.
#
#    real B(N), the right hand side vector.
#
#  Output:
#
#    real X(N), the solution vector.
#
  import numpy as np

  x = b.copy()

  for k in range ( 0, n - 1 ):

    l = ipvt[k]
    t = x[l]
    if ( l != k ):
      x[l] = x[k]
      x[k] = t

    x[k+1:n] = x[k+1:n] + t * ALU[k+1:n,k]

  for k in range ( n - 1, -1, -1 ):
    x[k] = x[k] / ALU[k,k]
    t = - x[k]
    x[0:k] = x[0:k] + t * ALU[0:k,k]

  return x

def dgesl_test ( ):

#*****************************************************************************80
#
## dgesl_test() tests dgesl(), after dgefa() has factored a small matrix.
#
#  Discussion:
#
#    Solve A*x = b where A is a given matrix, and B a right hand side.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgesl_test()' )
  print ( '  dgesl() solves a linear system involving a matrix' )
  print ( '  stored in general format,' )
  print ( '  after dgefa() has computed the LU factorization;' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  Factor the matrix.
#
  print ( '' )
  print ( '  dgefa() factors the matrix' )

  ALU, ipvt, info = dgefa ( A, n )

  if ( info != 0 ):
    print ( '  dgefa() returned an error flag INFO = ', info )
    return
#
#  Set the values of the right hand side vector B.
#
  b = np.array ( [ 14.0, 32.0, 23.0 ] )

  print ( '' )
  print ( '  The right hand side B' )
  print ( b )
#
#  Use DGESL to solve the system.
#
  print ( '' )
  print ( '  dgesl() solves the linear system.' )

  job = 0
  x = dgesl ( ALU, n, ipvt, b )

  print ( '' )
  print ( '  Computed solution X (should be (1,2,3))' )
  print ( x )

  return

def dgeslt ( ALU, n, ipvt, b ):

#*****************************************************************************80
#
## dgeslt() solves a transposed linear system A' * X = B.
#
#  Discussion:
#
#    The system matrix must have been factored by DGECO or DGEFA.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    This version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Moler, Bunch and Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Input:
#
#    real ALU(N,N), the LU factorization from DGECO or DGEFA.
#
#    integer N, the order of the matrix A.
#
#    integer IPVT(N), the pivot vector from DGECO or DGEFA.
#
#    real B(N), the right hand side vector.
#
#  Output:
#
#    real X(N), the solution vector.
#
  import numpy as np

  x = b.copy()

  for k in range ( 0, n ):
    t = np.dot ( x[0:k], ALU[0:k,k] )
    x[k] = ( x[k] - t ) / ALU[k,k]

  for k in range ( n - 2, -1, -1 ):

    x[k] = x[k] + np.dot ( x[k+1:n], ALU[k+1:n,k] )
    l = ipvt[k]

    if ( l != k ):
      t    = x[l]
      x[l] = x[k]
      x[k] = t

  return x

def dgeslt_test ( ):

#*****************************************************************************80
#
## dgeslt_test() tests dgeslt(), after dgefa() has factored a small matrix.
#
#  Discussion:
#
#    Solve A'*x = b where A is a given matrix, and B a right hand side.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 June 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3

  print ( '' )
  print ( 'dgeslt_test()' )
  print ( '  dgeslt() solves a transposed linear system involving a matrix' )
  print ( '  stored in general format,' )
  print ( '  after dgefa() has computed the LU factorization;' )
  print ( '' )
  print ( '  The number of equations is N = ', n )
#
#  Set the values of the matrix A.
#
  A = np.array ( [ \
    [ 1.0, 2.0, 3.0 ], \
    [ 4.0, 5.0, 6.0 ], \
    [ 7.0, 8.0, 0.0 ] ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  Factor the matrix.
#
  print ( '' )
  print ( '  dgefa() factors the matrix' )

  ALU, ipvt, info = dgefa ( A, n )

  if ( info != 0 ):
    print ( '  dgefa() returned an error flag INFO = ', info )
    return
#
#  Set the values of the right hand side vector B.
#
  b = np.array ( [ 30.0, 36.0, 15.0 ] )

  print ( '' )
  print ( '  The right hand side B' )
  print ( b )
#
#  Use DGESLT to solve the system.
#
  print ( '' )
  print ( '  dgeslt() solves the linear system.' )

  job = 0
  x = dgeslt ( ALU, n, ipvt, b )

  print ( '' )
  print ( '  Computed solution X (should be (1,2,3))' )
  print ( x )

  return

def dpofa ( a, lda, n ):

#*****************************************************************************80
#
## dpofa() factors a real symmetric positive definite matrix.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    FORTRAN77 version by Cleve Moler.
#    This version by John Burkardt.
#
#  Parameters:
#
#    Input, real A(LDA,N), the symmetric matrix to be factored.  Only the 
#    diagonal and upper triangle are accessed.
#
#    Input, integer LDA, the leading dimension of the array A.
#    N <= LDA.
#
#    Input, integer N, the order of the matrix.
#
#    Output, real A_FAC(LDA,N), an upper triangular matrix R such that
#    A = R' * R.  If INFO is nonzero, the factorization was not completed.
#
#    Output, integer INFO, error flag.
#    0, no error was detected.
#    K, the leading minor of order K is not positive definite.
#
  import numpy as np

  a_fac = a.copy ( )

  for i in range ( 1, n ):
    for j in range ( 0, i ):
      a_fac[i,j] = 0.0

  info = 0

  for j in range ( 0, n ):
    s = 0.0
    for k in range ( 0, j ):
      t = a_fac[k,j]
      for i in range ( 0, k ):
        t = t - a_fac[i,k] * a_fac[i,j]
      t = t / a_fac[k,k]
      a_fac[k,j] = t
      s = s + t * t

    s = a_fac[j,j] - s
    if ( s <= 0.0 ):
      info = j + 1
      return a_fac, info

    a_fac[j,j] = np.sqrt ( s )

  return a_fac, info

def dpofa_test ( ):

#*****************************************************************************80
#
## dpofa_test() tests dpofa().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 5
  lda = n

  print ( '' )
  print ( 'dpofa_test():' )
  print ( '  dpofa() computes the LU factors of a positive definite symmetric matrix,' )
#
#  Set the matrix A.
#
  a = np.zeros ( [ n, n ] )

  for i in range ( 0, n ):
    a[i,i] = 2.0
    if ( 0 < i ):
      a[i,i-1] = -1.0
    if ( i < n - 1 ):
      a[i,i+1] = -1.0

  print ( '' )
  print ( '  Matrix A:' )
  print ( a )
#
#  Factor the matrix.
#
  print ( '' )
  print ( '  Call DPOFA to factor the matrix.' )

  a_lu, info = dpofa ( a, lda, n )
 
  if ( info != 0 ):
    print ( '' )
    print ( '  Error, DPOFA returns INFO = %d' % ( info ) )
    return

  print ( '' )
  print ( '  Upper triangular factor U:' )
  print ( a_lu )

  uut = np.dot ( np.transpose ( a_lu ), a_lu )

  print ( '' )
  print ( '  Product Ut * U:' )
  print ( uut )

  return

def dqrdc ( a, lda, n, p, jpvt, job ):

#*****************************************************************************80
#
## dqrdc() computes the QR factorization of a real rectangular matrix.
#
#  Discussion:
#
#    The code uses Householder transformations.
#
#    Column pivoting based on the 2-norms of the reduced columns may be
#    performed at the user's option.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    This version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Moler, Bunch, Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Parameters:
#
#    Input, real A(LDA,P), the N by P matrix whose decomposition is to 
#    be computed.
#
#    Input, integer LDA, the leading dimension of the array A.  LDA must
#    be at least N.
#
#    Input, integer N, the number of rows of the matrix A.
#
#    Input, integer P, the number of columns of the matrix A.
#
#    Input, integer JPVT(P), integers that control the selection of the pivot 
#    columns.  The K-th column A(*,K) of A is placed in one of three classes 
#    according to the value of JPVT(K).
#      > 0, then A(K) is an initial column.
#      = 0, then A(K) is a free column.
#      < 0, then A(K) is a final column.
#    Before the decomposition is computed, initial columns are moved to
#    the beginning of the array A and final columns to the end.  Both
#    initial and final columns are frozen in place during the computation
#    and only free columns are moved.  At the K-th stage of the
#    reduction, if A(*,K) is occupied by a free column it is interchanged
#    with the free column of largest reduced norm.  JPVT is not referenced
#    if JOB == 0.  
#
#    Input, integer JOB, initiates column pivoting.
#    0, no pivoting is done.
#    nonzero, pivoting is done.
#
#    Output, real A(LDA,P), contains in its upper triangle the upper 
#    triangular matrix R of the QR factorization.  Below its diagonal 
#    A contains information from which the orthogonal part of the 
#    decomposition can be recovered.  Note that if pivoting has been 
#    requested, the decomposition is not that of the original matrix A 
#    but that of A with its columns permuted as described by JPVT.
#
#    Output, real QRAUX(P), contains further information required
#    to recover the orthogonal part of the decomposition.
#
#    Output, integer JPVT(P), JPVT(K) contains the index of the column of the
#    original matrix that has been interchanged into the K-th column, if
#    pivoting was requested.
#
  import numpy as np

  pl = 1
  pu = 0
  qraux = np.zeros ( p )
  work = np.zeros ( p )
#
#  If pivoting is requested, rearrange the columns.
#
  if ( job != 0 ):

    for j in range ( 0, p ):

      swapj = ( 0 < jpvt[j] )

      if ( jpvt[j] < 0 ):
        jpvt[j] = - ( j + 1 )
      else:
        jpvt[j] = ( j + 1 )

      if ( swapj ):

        if ( j + 1 != pl ):
          for i in range ( 0, n ):
            t = a[i,pl-1]
            a[i,pl-1] = a[i,j]
            a[i,j] = t

        jpvt[j] = jpvt[pl-1]
        jpvt[pl-1] = j + 1
        pl = pl + 1

    pu = p

    for j in range ( p - 1, -1, -1 ):

      if ( jpvt[j] < 0 ):

        jpvt[j] = - jpvt[j]

        if ( j + 1 != pu ):

          for i in range ( 0, n ):
            t = a[i,pu-1]
            a[i,pu-1] = a[i,j]
            a[i,j] = t

          jp = jpvt[pu-1]
          jpvt[pu-1] = jpvt[j]
          jpvt[j] = jp

        pu = pu - 1
#
#  Compute the norms of the free columns.
#
  for j in range ( pl - 1, pu ):
    t = 0.0
    for i in range ( 0, n ):
      t = t + a[i,j] ** 2
    qraux[j] = np.sqrt ( t )
    work[j] = qraux[j]
#
#  Perform the Householder reduction of A.
#
  lup = min ( n, p )

  for l in range ( 0, lup ):
#
#  Bring the column of largest norm into the pivot position.
#
    if ( pl <= l + 1 and l + 1 < pu ):

      maxnrm = 0.0
      maxj = l

      for j in range ( l, pu ):

        if ( maxnrm < qraux[j] ):
          maxnrm = qraux[j]
          maxj = j

      if ( maxj != l ):

        for i in range ( 0, n ):
          t         = a[i,l]
          a[i,l]    = a[i,maxj]
          a[i,maxj] = t

        qraux[maxj] = qraux[l]
        work[maxj] = work[l]

        jp = jpvt[maxj]
        jpvt[maxj] = jpvt[l]
        jpvt[l] = jp
#
#  Compute the Householder transformation for column L.
#
    qraux[l] = 0.0

    if ( l + 1 != n ):

      t = 0.0
      for i in range ( l, n ):
        t = t + a[i,l] ** 2
      nrmxl = np.sqrt ( t )

      if ( nrmxl != 0.0 ):

        if ( a[l,l] < 0.0 ):
          nrmxl = - abs ( nrmxl )
        elif ( 0.0 < a[l,l] ):
          nrmxl = abs ( nrmxl )

        for i in range ( l, n ):
          a[i,l] = a[i,l] / nrmxl

        a[l,l] = 1.0 + a[l,l]
#
#  Apply the transformation to the remaining columns, updating the norms.
#
        for j in range ( l + 1, p ):

          t = 0.0
          for i in range ( l, n ):
            t = t - a[i,l] * a[i,j]
          t = t / a[l,l]

          for i in range ( l, n ):
            a[i,j] = a[i,j] + t * a[i,l]

          if ( pl <= j and j <= pu ):

            if ( qraux[j] != 0.0 ):

              tt = 1.0 - ( abs ( a[l,j] ) / qraux[j] ) ** 2
              tt = max ( tt, 0.0 )
              t = tt
              tt = 1.0 + 0.05 * tt * ( qraux[j] / work[j] ) ** 2

              if ( tt != 1.0 ):
                qraux[j] = qraux[j] * np.sqrt ( t )
              else:
                t = 0.0
                for i in range ( l + 1, n ):
                  t = t + a[i,j] ** 2
                qraux[j] = np.sqrt ( t )
                work[j] = qraux[j]
#
#  Save the transformation.
#
        qraux[l] = a[l,l]
        a[l,l] = - nrmxl

  return a, qraux, jpvt

def dqrdc_test ( ):

#*****************************************************************************80
#
## dqrdc_test() tests dqrdc().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 3
  p = 3
  lda = n

  print ( '' )
  print ( 'dqrdc_test():' )
  print ( '  dqrdc() computes the QR decomposition of a rectangular' )
  print ( '  matrix, but does not return Q and R explicitly.' )
  print ( '' )
  print ( '  Show how Q and R can be recovered using DQRSL.' )
#
#  Set the matrix A.
#
  a = np.array \
  ( [ \
    [ 1.0, 1.0, 0.0 ], \
    [ 1.0, 0.0, 1.0 ], \
    [ 0.0, 1.0, 1.0 ]  \
  ] )

  print ( '' )
  print ( '  The original matrix A:' )
  print ( a )
#
#  Decompose the matrix.
#
  print ( '' )
  print ( '  Decompose the matrix.' )

  job = 0
  ipvt = np.zeros ( p, dtype = np.int32 )

  a, qraux, ipvt = dqrdc ( a, lda, n, p, ipvt, job )
#
#  Print out what DQRDC has stored in A...
#
  print ( '' )
  print ( '  The packed matrix A which describes Q and R:' )
  print ( a )
#
#  ...and in QRAUX.
#
  print ( '' )
  print ( '  The QRAUX vector, containing some additional' )
  print ( '  information defining Q:' )
  print ( qraux )
#
#  Print out the resulting R factor.
#
  r = np.zeros ( [ n, p ] )

  for i in range ( 0, n ):
    for j in range ( 0, p ):
      if ( i <= j ):
        r[i,j] = a[i,j]

  print ( '' )
  print ( '  The R factor:' )
  print ( r )
#
#  Call DQRSL to extract the information about the Q matrix.
#  We do this, essentially, by asking DQRSL to tell us the
#  value of Q*Y, where Y is a column of the identity matrix.
#
  q = np.zeros ( [ n, n ] )

  job = 10000

  for j in range ( 0, n ):
#
#  Set the vector Y.
#
    y = np.zeros ( n )

    y[j] = 1.0
#
#  Ask DQRSL to tell us what Q*Y is.
#
    qy, qty, b, rsd, xb, info = dqrsl ( a, lda, n, p, qraux, y, job )

    if ( info != 0 ):
      print ( '  Error!  DQRSL returns INFO = %d' % ( info ) )
      return
#
#  Copy QY into the appropriate column of Q.
#
    for i in range ( 0, n ):
      q[i,j] = qy[i]
#
#  Now print out the Q matrix we have extracted.
#
  print ( '' )
  print ( '  The Q factor:' )
  print ( q )
#
#  Compute Q*R to verify that it equals A.
#
  b = np.dot ( q, r )
#
#  Print the result.
#
  print ( '' )
  print ( '  The product Q * R:' )
  print ( b )

  return

def dqrsl ( a, lda, n, k, qraux, y, job ):

#*****************************************************************************80
#
## dqrsl() computes transformations, projections, and least squares solutions.
#
#  Discussion:
#
#    The code requires the output of dqrdc().
#
#    For K <= min(N,P), let AK be the matrix
#
#      AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
#
#    formed from columns JPVT(1), ..., JPVT(K) of the original
#    N by P matrix A that was input to DQRDC.  If no pivoting was
#    done, AK consists of the first K columns of A in their
#    original order.  DQRDC produces a factored orthogonal matrix Q
#    and an upper triangular matrix R such that
#
#      AK = Q * (R)
#               (0)
#
#    This information is contained in coded form in the arrays
#    A and QRAUX.
#
#    The parameters QY, QTY, B, RSD, and AB are not referenced
#    if their computation is not requested and in this case
#    can be replaced by dummy variables in the calling program.
#    To save storage, the user may in some cases use the same
#    array for different parameters in the calling sequence.  A
#    frequently occuring example is when one wishes to compute
#    any of B, RSD, or AB and does not need Y or QTY.  In this
#    case one may identify Y, QTY, and one of B, RSD, or AB, while
#    providing separate arrays for anything else that is to be
#    computed.
#
#    Thus the calling sequence
#
#      call dqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
#
#    will result in the computation of B and RSD, with RSD
#    overwriting Y.  More generally, each item in the following
#    list contains groups of permissible identifications for
#    a single calling sequence.
#
#      1. (Y,QTY,B) (RSD) (AB) (QY)
#
#      2. (Y,QTY,RSD) (B) (AB) (QY)
#
#      3. (Y,QTY,AB) (B) (RSD) (QY)
#
#      4. (Y,QY) (QTY,B) (RSD) (AB)
#
#      5. (Y,QY) (QTY,RSD) (B) (AB)
#
#      6. (Y,QY) (QTY,AB) (B) (RSD)
#
#    In any group the value returned in the array allocated to
#    the group corresponds to the last member of the group.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    This version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Moler, Bunch, Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Parameters:
#
#    Input, real A(LDA,P), contains the output of DQRDC.
#
#    Input, integer LDA, the leading dimension of the array A.
#
#    Input, integer N, the number of rows of the matrix AK.  It must
#    have the same value as N in DQRDC.
#
#    Input, integer K, the number of columns of the matrix AK.  K
#    must not be greater than min(N,P), where P is the same as in the
#    calling sequence to DQRDC.
#
#    Input, real QRAUX(P), the auxiliary output from DQRDC.
#
#    Input, real Y(N), a vector to be manipulated by DQRSL.
#
#    Input, integer JOB, specifies what is to be computed.  JOB has
#    the decimal expansion ABCDE, with the following meaning:
#      if A /= 0, compute QY.
#      if B /= 0, compute QTY.
#      if C /= 0, compute QTY and B.
#      if D /= 0, compute QTY and RSD.
#      if E /= 0, compute QTY and AB.
#    Note that a request to compute B, RSD, or AB automatically triggers
#    the computation of QTY, for which an array must be provided in the
#    calling sequence.
#
#    Output, real QY(N), contains Q * Y, if requested.
#
#    Output, real QTY(N), contains Q' * Y, if requested.
#
#    Output, real B(K), the solution of the least squares problem
#      minimize norm2 ( Y - AK * B),
#    if its computation has been requested.  Note that if pivoting was
#    requested in DQRDC, the J-th component of B will be associated with
#    column JPVT(J) of the original matrix A that was input into DQRDC.
#
#    Output, real RSD(N), the least squares residual Y - AK * B,
#    if its computation has been requested.  RSD is also the orthogonal
#    projection of Y onto the orthogonal complement of the column space
#    of AK.
#
#    Output, real AB(N), the least squares approximation Ak * B,
#    if its computation has been requested.  AB is also the orthogonal
#    projection of Y onto the column space of A.
#
#    Output, integer INFO, is zero unless the computation of B has
#    been requested and R is exactly singular.  In this case, INFO is the
#    index of the first zero diagonal element of R, and B is left unaltered.
#
  import numpy as np

  qy = np.zeros ( n )
  qty = np.zeros ( n )
  b = np.zeros ( k )
  rsd = np.zeros ( n )
  ab = np.zeros ( n )
#
#  Set info flag.
#
  info = 0
#
#  Determine what is to be computed.
#
  cqy =  int   ( job / 10000         ) != 0
  cqty =       ( job % 10000 )         != 0
  cb =   int ( ( job %  1000 ) / 100 ) != 0
  cr =   int ( ( job %   100 ) /  10 ) != 0
  cab =        ( job %    10 )         != 0

  ju = min ( k, n - 1 )
#
#  Special action when N = 1.
#
  if ( ju == 0 ):

    qy[0] = y[0]
    qty[0] = y[0]
    ab[0] = y[0]

    if ( a[0,0] == 0.0 ):
      info = 1
    else:
      b[0] = y[0] / a[0]

    rsd[0] = 0.0

    return qy, qty, b, rsd, ab, info
#
#  Set up to compute QY or QTY.
#
  for i in range ( 0, n ):
    qy[i] = y[i]
    qty[i] = y[i]
#
#  Compute QY.
#
  if ( cqy ):

    for j in range ( ju - 1, -1, -1 ):

      if ( qraux[j] != 0.0 ):
        ajj = a[j,j]
        a[j,j] = qraux[j]
        t = 0.0
        for i in range ( j, n ):
          t = t - a[i,j] * qy[i]
        t = t / a[j,j]
        for i in range ( j, n ):
          qy[i] = qy[i] + t * a[i,j]
        a[j,j] = ajj
#
#  Compute Q'*Y.
#
  if ( cqty ):

    for j in range ( 0, ju ):
      if ( qraux[j] != 0.0 ):
        ajj = a[j,j]
        a[j,j] = qraux[j]
        t = 0.0
        for i in range ( j, n ):
          t = t - a[i,j] * qty[i]
        t = t / a[j,j]
        for i in range ( j, n ):
          qty[i] = qty[i] + t * a[i,j]
        a[j,j] = ajj
#
#  Set up to compute B, RSD, or AB.
#
  for i in range ( 0, k ):
    b[i] = qty[i]
    ab[i] = qty[i]

  for i in range ( k, n ):
    rsd[i] = qty[i]
#
#  Compute B.
#
  if ( cb ):

    for j in range ( k - 1, -1, -1 ):

      if ( a[j,j] == 0.0 ):
        info = j + 1
        break

      b[j] = b[j] / a[j,j]

      if ( 0 < j ):
        t = - b[j]
        for i in range ( 0, j ):
          b[i] = b[i] + t * a[i,j]
#
#  Compute RSD or AB as required.
#
  if ( cr or cab ):

    for j in range ( ju - 1, -1, -1 ):

      if ( qraux[j] != 0.0 ):

        ajj = a[j,j]
        a[j,j] = qraux[j]

        if ( cr ):
          t = 0.0
          for i in range ( j, n ):
            t = t - a[i,j] * rsd[i]
          t = t / a[j,j]
          for i in range ( j, n ):
            rsd[i] = rsd[i] + t * a[i,j]

        if ( cab ):
          t = 0.0
          for i in range ( j, n ):
            t = t - a[i,j] * ab[i]
          t = t / a[j,j]
          for i in range ( j, n ):
            ab[i] = ab[i] + t * a[i,j]

        a[j,j] = ajj

  return qy, qty, b, rsd, ab, info

def dqrsl_test ( ):

#*****************************************************************************80
#
## dqrsl_test() tests dqrsl().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 5
  p = 3
  lda = n

  print ( '' )
  print ( 'dqrsl_test():' )
  print ( '  dqrsl() solves a rectangular linear system A*x=b in the' )
  print ( '  least squares sense after A has been factored by DQRDC.' )
#
#  Set the matrix A.
#
  a = np.array \
  ( [ \
    [ 1.0,  1.0,  1.0 ], \
    [ 1.0,  2.0,  4.0 ], \
    [ 1.0,  3.0,  9.0 ], \
    [ 1.0,  4.0, 16.0 ], \
    [ 1.0,  5.0, 25.0 ]  \
  ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( a )
#
#  Decompose the matrix.
#
  print ( '' )
  print ( '  Decompose the matrix.' )

  job = 0
  ipvt = np.zeros ( n, dtype = np.int32 )

  a, qraux, ipvt = dqrdc ( a, lda, n, p, ipvt, job )
#
#  Call DQRSL to compute the least squares solution A*x=b.
#
  job = 110

  y = np.array \
  ( [ \
    1.0, \
    2.3, \
    4.6, \
    3.1, \
    1.2  \
  ] )  

  b2 = np.array \
  ( [ \
    -3.02, \
     4.4914286, \
    -0.72857143 
  ] )

  qy, qty, b, r, xb, info = dqrsl ( a, lda, n, p, qraux, y, job )

  if ( info != 0 ):
    print ( '' )
    print ( 'DQRSL_TEST - Warning!' )
    print ( '  DQRSL returns INFO = %d' % ( info ) )
    return

  print ( '' )
  print ( '      X          X(expected):' )
  print ( '' )

  for i in range ( 0, p ):
    print ( '  %14.6g  %14.6g' % ( b[i], b2[i] ) )

  return

def drotg ( a, b ):

#*****************************************************************************80
#
## drotg() constructs a Givens plane rotation.
#
#  Discussion:
#
#    Given values A and B, this routine computes
#
#    SIGMA = sign ( A ) if abs ( A ) >  abs ( B )
#          = sign ( B ) if abs ( A ) <= abs ( B )
#
#    R     = SIGMA * ( A * A + B * B )
#
#    C = A / R if R is not 0
#      = 1     if R is 0
#
#    S = B / R if R is not 0,
#        0     if R is 0.
#
#    The computed numbers then satisfy the equation
#
#    (  C  S ) ( A ) = ( R )
#    ( -S  C ) ( B ) = ( 0 )
#
#    The routine also computes
#
#    Z = S     if abs ( A ) > abs ( B, end = '' )
#      = 1 / C if abs ( A ) <= abs ( B ) and C is not 0,
#      = 1     if C is 0.
#
#    The single value Z encodes C and S, and hence the rotation:
#
#    If Z = 1, set C = 0 and S = 1
#    If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z
#    if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C )
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 September 2016
#
#  Author:
#
#    This version by John Burkardt.
#
#  Reference:
#
#    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
#    LINPACK User's Guide,
#    SIAM, 1979.
#
#    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
#    Basic Linear Algebra Subprograms for Fortran Usage,
#    Algorithm 539,
#    ACM Transactions on Mathematical Software,
#    Volume 5, Number 3, September 1979, pages 308-323.
#
#  Parameters:
#
#    Input, real A, B, the values A and B.
#
#    Output, real C, S, the cosine and sine of the Givens rotation.
#
#    Output, real R, Z, the values R and Z.
#
  import numpy as np

  if ( abs ( b ) < abs ( a ) ):
    roe = a
  else:
    roe = b

  scale = abs ( a ) + abs ( b )

  if ( scale == 0.0 ):
    c = 1.0
    s = 0.0
    r = 0.0
  else:
    r = scale * np.sqrt ( ( a / scale ) ** 2 + ( b / scale ) ** 2 )
    if ( roe < 0.0 ):
      r = - r
    c = a / r
    s = b / r

  if ( 0.0 < abs ( c ) and abs ( c ) <= s ):
    z = 1.0 / c
  else:
    z = s

  return c, s, r, z

def drotg_test ( rng ):

#*****************************************************************************80
#
## drotg_test() tests drotg().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 September 2016
#
#  Author:
#
#    John Burkardt
#
  test_num = 5

  print ( '' )
  print ( 'drotg_test():' )
  print ( '  drotg() generates a real Givens rotation' )
  print ( '    (  C  S ) * ( A ) = ( R )' )
  print ( '    ( -S  C )   ( B )   ( 0 )' )
  print ( '' )

  for test in range ( 0, test_num ):

    a = rng.random ( )
    b = rng.random ( )

    c, s, r, z = drotg ( a, b )

    print ( '' )
    print ( '  A =  %g  B =  %g' % ( a, b ) )
    print ( '  C =  %g  S =  %g' % ( c, s ) )
    print ( '  R =  %g  Z =  %g' % ( r, z ) )
    print ( '   C*A+S*B = %g' % (  c * a + s * b ) )
    print ( '  -S*A+C*B = %g' % ( -s * a + c * b ) )

  return

def drot ( n, x, incx, y, incy, c, s ):

#*****************************************************************************80
#
## drot() applies a plane rotation.
#
#  Discussion:
#
#    Note that the FORTRAN version of this function allowed users to pass in
#    X and Y data that was noncontiguous, (such as rows of a FORTRAN matrix).
#    The rotated data overwrote the input data, and so it might therefore
#    also be noncontiguous.
#
#    This function does not assume that the output overwrites the input,
#    and treats the output vectors as new items of length exactly N.
#
#    Note, moreover, that Python does NOT allow a matrix to be treated as a 
#    vector in quite the simple way that FORTRAN does.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 September 2016
#
#  Author:
#
#    This version by John Burkardt.
#
#  Reference:
#
#    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
#    LINPACK User's Guide,
#    SIAM, 1979.
#
#    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
#    Basic Linear Algebra Subprograms for Fortran Usage,
#    Algorithm 539,
#    ACM Transactions on Mathematical Software,
#    Volume 5, Number 3, September 1979, pages 308-323.
#
#  Parameters:
#
#    Input, integer N, the number of entries in the vectors.
#
#    Input, real X(INCX*N), one of the vectors to be rotated.
#
#    Input, integer INCX, the increment between successive entries of X.
#
#    Input, real Y(INCX*N), one of the vectors to be rotated.
#
#    Input, integer INCY, the increment between successive elements of Y.
#
#    Input, real C, S, parameters (presumably the cosine and
#    sine of some angle) that define a plane rotation.
#
#    Output, real XR(N), the rotated vector.
#
#    Output, real YR(N), the rotated vector.
#
  import numpy as np

  xr = np.zeros ( n )
  yr = np.zeros ( n )

  if ( n <= 0 ):

    pass

  elif ( incx == 1 and incy == 1 ):

    xr[0:n] = c * x[0:n] + s * y[0:n]
    yr[0:n] = c * y[0:n] - s * x[0:n]

  else:

    if ( 0 <= incx ):
      ix = 0
    else:
      ix = ( - n + 1 ) * incx

    if ( 0 <= incy ):
      iy = 0
    else:
      iy = ( - n + 1 ) * incy

    for i in range ( 0, n ):
      xr[ix] = c * x[ix] + s * y[iy]
      yr[iy] = c * y[iy] - s * x[ix]
      ix = ix + incx
      iy = iy + incy

  return xr, yr

def drot_test ( ):

#*****************************************************************************80
#
## drot_test() tests drot().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 September 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
 
  n = 6
  x = np.zeros ( n )
  y = np.zeros ( n )

  for i in range ( 0, n ):
    x[i] = float ( i + 1 )
    y[i] = x[i] * x[i] - 12.0
 
  print ( '' )
  print ( 'drot_test():' )
  print ( '  drot() carries out a Givens rotation.' )
  print ( '' )
  print ( '  Vectors X and Y' )
  print ( '' )
  for i in range ( 0, n ):
    print ( '  %6d  %14.6g  %14.6g' % ( i, x[i], y[i] ) )

  c = 0.5
  s = np.sqrt ( 1.0 - c * c )
  xr, yr = drot ( n, x, 1, y, 1, c, s )
  print ( '' )
  print ( '  xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) )
  print ( '' )
  print ( '  Rotated vectors XR and YR' )
  print ( '' )
  for i in range ( 0, n ):
    print ( '  %6d  %14.6g  %14.6g' % ( i, xr[i], yr[i] ) )

  t = np.arctan2 ( y[0], x[0] )
  c = np.cos ( t )
  s = np.sin ( t )
  xr, yr = drot ( n, x, 1, y, 1, c, s )
  print ( '' )
  print ( '  xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) )
  print ( '' )
  print ( '  Rotated vectors XR and YR' )
  print ( '' )
  for i in range ( 0, n ):
    print ( '  %6d  %14.6g  %14.6g' % ( i, xr[i], yr[i] ) )

  return

def dsvdc ( a, lda, m, n, ldu, ldv, job ):

#*****************************************************************************80
#
## dsvdc() computes the singular value decomposition of a real rectangular matrix.
#
#  Discussion:
#
#    This routine reduces an M by N matrix A to diagonal form by orthogonal
#    transformations U and V.  The diagonal elements S(I) are the singular
#    values of A.  The columns of U are the corresponding left singular
#    vectors, and the columns of V the right singular vectors.
#
#    The form of the singular value decomposition is then
#
#      A(MxN) = U(MxM) * S(MxN) * V(NxN)'
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    This version by John Burkardt.
#
#  Reference:
#
#    Dongarra, Moler, Bunch and Stewart,
#    LINPACK User's Guide,
#    SIAM, (Society for Industrial and Applied Mathematics),
#    3600 University City Science Center,
#    Philadelphia, PA, 19104-2688.
#    ISBN 0-89871-172-X
#
#  Parameters:
#
#    Input, real A(LDA,N), the M by N matrix whose singular 
#    value decomposition is to be computed.
#
#    Input, integer LDA, the leading dimension of the array A.
#    LDA must be at least N.
#
#    Input, integer M, the number of rows of the matrix.
#
#    Input, integer N, the number of columns of the matrix A.
#
#    Input, integer LDU, the leading dimension of the array U.
#    LDU must be at least M.
#
#    Input, integer LDV, the leading dimension of the array V.
#    LDV must be at least N.
#
#    Input, integer JOB, controls the computation of the singular
#    vectors.  It has the decimal expansion AB with the following meaning:
#      A =  0, for not compute the left singular vectors.
#      A =  1, return the M left singular vectors in U.
#      A >= 2, return the first min(M,N) singular vectors in U.
#      B =  0, for not compute the right singular vectors.
#      B =  1, return the right singular vectors in V.
#
#    Output, real A(LDA,N), the matrix has been destroyed.  Depending on the 
#    user's requests, the matrix may contain other useful information.
#
#    Output, real S(MM), where MM = max(M+1,N).  The first
#    min(M,N) entries of S contain the singular values of A arranged in
#    descending order of magnitude.
#
#    Output, real E(MM), where MM = max(M+1,N), ordinarily contains zeros.
#    However see the discussion of INFO for exceptions.
#
#    Output, real U(LDU,K).  If JOBA = 1 then K = M
#    if 2 <= JOBA, then K = min(M,N).  U contains the M by M matrix of
#    left singular vectors.  U is not referenced if JOBA = 0.  If M <= N
#    or if JOBA = 2, then U may be identified with A in the subroutine call.
#
#    Output, real V(LDV,N), the N by N matrix of right singular
#    vectors.  V is not referenced if JOB is 0.  If N <= M, then V may be
#    identified with A in the subroutine call.
#
#    Output, integer INFO, status indicator.
#    The singular values (and their corresponding singular vectors)
#    S(INFO+1), S(INFO+2),...,S(MN) are correct.  Here MN = min ( M, N ).
#    Thus if INFO is 0, all the singular values and their vectors are
#    correct.  In any event, the matrix B = U' * A * V is the bidiagonal
#    matrix with the elements of S on its diagonal and the elements of E on
#    its superdiagonal.  Thus the singular values of A and B are the same.
#
  import numpy as np

  maxit = 30
  mm = max ( m + 1, n )

  s = np.zeros ( mm )
  e = np.zeros ( mm )
  joba = ( job // 10 )
  if ( joba == 1 ):
    u = np.zeros ( [ ldu, m ] )
  elif ( 1 <= joba ):
    u = np.zeros ( [ ldu, min ( m, n ) ] )
  v = np.zeros ( [ ldv, n ] )
  work = np.zeros ( m )

  info = 0
#
#  Determine what is to be computed.
#
  wantu = 0
  wantv = 0
  jobu = int ( ( job % 100 ) / 10 )

  if ( 1 < jobu ):
    ncu = min ( m, n )
  else:
    ncu = m

  if ( jobu != 0 ):
    wantu = 1

  if ( ( job % 10 ) != 0 ):
    wantv = 1
#
#  Reduce A to bidiagonal form, storing the diagonal elements
#  in S and the super-diagonal elements in E.
#
  info = 0
  nct = min ( m - 1, n )
  nrt = max ( 0, min ( m, n - 2 ) )
  lu = max ( nct, nrt )

  for l in range ( 0, lu ):
#
#  Compute the transformation for the L-th column and
#  place the L-th diagonal in S(L).
#
    if ( l <= nct - 1 ):

      t = 0.0
      for i in range ( l, m ):
        t = t + a[i,l] ** 2
      t = np.sqrt ( t )
      s[l] = t

      if ( s[l] != 0.0 ):
        if ( a[l,l] < 0.0 ):
          s[l] = - s[l]
        for i in range ( l, m ):
          a[i,l] = a[i,l] / s[l]
        a[l,l] = 1.0 + a[l,l]

      s[l] = - s[l]

    for j in range ( l + 1, n ):
#
#  Apply the transformation.
#
      if ( l <= nct - 1 and s[l] != 0.0 ):
        t = 0.0
        for i in range ( l, m ):
          t = t - a[i,l] * a[i,j]
        t = t / a[l,l]
        for i in range ( l, m ):
          a[i,j] = a[i,j] + t * a[i,l]
#
#  Place the L-th row of A into E for the
#  subsequent calculation of the row transformation.
#
      e[j] = a[l,j]
#
#  Place the transformation in U for subsequent back multiplication.
#
    if ( wantu and l <= nct - 1 ):
      for i in range ( l, m ):
        u[i,l] = a[i,l]
#
#  Compute the L-th row transformation and place the
#  L-th superdiagonal in E(L).
#
    if ( l <= nrt - 1 ):

      t = 0.0
      for i in range ( l + 1, n ):
        t = t + e[i] ** 2
      e[l] = np.sqrt ( t )
      if ( e[l] != 0.0 ):
        if ( e[l+1] < 0.0 ):
          e[l] = - e[l]
        e[l+1:n] = e[l+1:n] / e[l]
        e[l+1] = 1.0 + e[l+1]
      e[l] = - e[l]
#
#  Apply the transformation.
#
      if ( l + 1 <= m - 1 and e[l] != 0.0 ):

        work[l+1:m] = 0.0

        for j in range ( l + 1, n ):
          for i in range ( l + 1, m ):
            work[i] = work[i] + e[j] * a[i,j]

        for j in range ( l + 1, n ):
          for i in range ( l + 1, m ):
            a[i,j] = a[i,j] - e[j] / e[l+1] * work[i]
#
#  Place the transformation in V for subsequent back multiplication.
#
      if ( wantv ):
        v[l+1:n,l] = e[l+1:n]
#
#  Set up the final bidiagonal matrix of order MN.
#
  mn = min ( m + 1, n )
  nctp1 = nct + 1
  nrtp1 = nrt + 1

  if ( nct < n ):
    s[nctp1-1] = a[nctp1-1,nctp1-1]

  if ( m < mn ):
    s[mn-1] = 0.0

  if ( nrtp1 < mn ):
    e[nrtp1-1] = a[nrtp1-1,mn-1]

  e[mn-1] = 0.0
#
#  If required, generate U.
#
  if ( wantu ):

    for j in range ( nctp1 - 1, ncu ):
      for i in range ( 0, m ):
        u[i,j] = 0.0

    for j in range ( nctp1 - 1, ncu ):
      u[j,j] = 1.0

    for l in range ( nct - 1, -1, -1 ):
 
      if ( s[l] != 0.0 ):

        for j in range ( l + 1, ncu ):
          t = 0.0
          for i in range ( l, m ):
            t = t - u[i,l] * u[i,j]
          t = t / u[l,l]
          for i in range ( l, m ):
            u[i,j] = u[i,j] + t * u[i,l]

        for i in range ( 0, l ):
          u[i,l] = 0.0
        for i in range ( l, m ):
          u[i,l] = - u[i,l]
        u[l,l] = 1.0 + u[l,l]

      else:

        for i in range ( 0, m ):
          u[i,l] = 0.0
        u[l,l] = 1.0
#
#  If it is required, generate V.
#
  if ( wantv ):

    for l in range ( n - 1, -1, -1 ):

      if ( l <= nrt - 1 and e[l] != 0.0 ):

        for j in range ( l + 1, n ):
          t = 0.0
          for i in range ( l + 1, n ):
            t = t - v[i,l] * v[i,j]
          t = t / v[l+1,l]
          for i in range ( l + 1, n ):
            v[i,j] = v[i,j] + t * v[i,l]

      for i in range ( 0, n ):
        v[i,l] = 0.0
      v[l,l] = 1.0
#
#  Main iteration loop for the singular values.
#
  mm = mn
  iter = 0

  while ( 0 < mn ):
#
#  If too many iterations have been performed, set flag and return.
#
    if ( maxit <= iter ):
      print ( '' )
      print ( 'DSVDC - Warning!' )
      print ( '  MAXIT = %d <= ITER = %d' % ( maxit, iter ) )
      info = mn
      return a, s, e, u, v, info
#
#  This section of the program inspects for
#  negligible elements in the S and E arrays.
#
#  On completion the variables KASE and L are set as follows:
#
#  KASE = 1     if S(MN) and E(L-1) are negligible and L < MN
#  KASE = 2     if S(L) is negligible and L < MN
#  KASE = 3     if E(L-1) is negligible, L < MN, and
#               S(L), ..., S(MN) are not negligible (QR step).
#  KASE = 4     if E(MN-1) is negligible (convergence).
#
    for l in range ( mn - 2, -2, -1 ):

      if ( l == -1 ):
        break

      test = abs ( s[l] ) + abs ( s[l+1] )
      ztest = test + abs ( e[l] )

      if ( ztest == test ):
        e[l] = 0.0
        break

    if ( l == mn - 2 ):

      kase = 4

    else:

      for ls in range ( mn - 1, l - 1, -1 ):

        if ( ls == l ):
          break

        test = 0.0
        if ( ls != mn - 1 ):
          test = test + abs ( e[ls] )

        if ( ls != l + 1 ):
          test = test + abs ( e[ls-1] )

        ztest = test + abs ( s[ls] )

        if ( ztest == test ):
          s[ls] = 0.0
          break

      if ( ls == l ):
        kase = 3
      elif ( ls == mn - 1 ):
        kase = 1
      else:
        kase = 2
        l = ls

    l = l + 1
#
#  Deflate negligible S(MN).
#
    if ( kase == 1 ):

      mm1 = mn - 1
      f = e[mn-2]
      e[mn-2] = 0.0

      for k in range ( mm1 - 1, l - 1, -1 ):
 
        t1 = s[k]
        cs, sn, t1, f = drotg ( t1, f )
        s[k] = t1

        if ( k != l ):
          f = - sn * e[k-1]
          e[k-1] = cs * e[k-1]

        if ( wantv ):
          v1, v2 = drot ( n, v[0:n,k], 1, v[0:n,mn-1], 1, cs, sn )
          for i in range ( 0, n ):
            v[i,k]    = v1[i]
            v[i,mn-1] = v2[i]
#
#  Split at negligible S(L).
#
    elif ( kase == 2 ):

      f = e[l-1]
      e[l-1] = 0.0

      for k in range ( l, mn ):

        t1 = s[k]
        cs, sn, t1, f = drotg ( t1, f )
        s[k] = t1
        f = - sn * e[k]
        e[k] = cs * e[k]
        if ( wantu ):
          u1, u2 = drot ( m, u[0:m,k], 1, u[0:m,l-1], 1, cs, sn )
          for i in range ( 0, m ):
            u[i,k] = u1[i]
            u[i,l-1] = u2[i]
#
#  Perform one QR step.
#
    elif ( kase == 3 ):
#
#  Calculate the shift.
#
      scale = max ( abs ( s[mn-1] ), \
              max ( abs ( s[mn-2] ), \
              max ( abs ( e[mn-2] ), \
              max ( abs ( s[l]    ), \
                    abs ( e[l] ) ) ) ) )

      sm = s[mn-1] / scale
      smm1 = s[mn-2] / scale
      emm1 = e[mn-2] / scale
      sl = s[l] / scale
      el = e[l] / scale
      b = ( ( smm1 + sm ) * ( smm1 - sm ) + emm1 * emm1 ) / 2.0
      c = sm  * sm * emm1 * emm1
      shift = 0.0

      if ( b != 0.0 or c != 0.0 ):
        shift = np.sqrt ( b * b + c )
        if ( b < 0.0 ):
          shift = - shift
        shift = c / ( b + shift )

      f = ( sl + sm ) * ( sl - sm ) - shift
      g = sl * el
#
#  Chase zeros.
#
      mm1 = mn - 1

      for k in range ( l, mm1 ):

        cs, sn, f, g = drotg ( f, g )

        if ( k != l ):
          e[k-1] = f

        f = cs * s[k] + sn * e[k]
        e[k] = cs * e[k] - sn * s[k]
        g = sn * s[k+1]
        s[k+1] = cs * s[k+1]

        if ( wantv ):
          v1, v2 = drot ( n, v[0:n,k], 1, v[0:n,k+1], 1, cs, sn )
          for i in range ( 0, n ):
            v[i,k] = v1[i]
            v[i,k+1] = v2[i]

        cs, sn, f, g = drotg ( f, g )
        s[k] = f
        f = cs * e[k] + sn * s[k+1]
        s[k+1] = - sn * e[k] + cs * s[k+1]
        g = sn * e[k+1]
        e[k+1] = cs * e[k+1]

        if ( wantu and k < m ):
          u1, u2 = drot ( m, u[0:m,k], 1, u[0:m,k+1], 1, cs, sn )
          for i in range ( 0, m ):
            u[i,k] = u1[i]
            u[i,k+1] = u2[i]

      e[mn-2] = f

      iter = iter + 1
#
#  Convergence.
#
    elif ( kase == 4 ):
#
#  Make the singular value nonnegative.
#
      if ( s[l] < 0.0 ):
        s[l] = - s[l]
        if ( wantv ):
          for i in range ( 0, n ):
            v[i,l] = - v[i,l]
#
#  Order the singular values.
#
      while ( True ):

        if ( l + 1 == mm ):
          break

        if ( s[l+1] <= s[l] ):
          break

        t = s[l]
        s[l] = s[l+1]
        s[l+1] = t

        if ( wantv and l < n - 1 ):
          for i in range ( 0, n ):
            t = v[i,l]
            v[i,l] = v[i,l+1]
            v[i,l+1] = t
 
        if ( wantu and l < m - 1 ):
          for i in range ( 0, m ):
            t = u[i,l]
            u[i,l] = u[i,l+1]
            u[i,l+1] = t
 
        l = l + 1

      iter = 0
      mn = mn - 1

  return a, s, e, u, v, info

def dsvdc_test ( rng ):

#*****************************************************************************80
#
## dsvdc_test() tests dsvdc().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 September 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  m = 6
  n = 4

  print ( '' )
  print ( 'dsvdc_test():' )
  print ( '  dsvdc() computes the singular value decomposition' )
  print ( '  for an MxN matrix A in general storage.' )
  print ( '    A = U * S * V\'' )
  print ( '' )
  print ( '  Matrix rows M =    %d' % ( m ) )
  print ( '  Matrix columns N = %d' % ( n ) )
#
#  Set A.
#
  A = rng.random ( size = [ m, n ] )

  print ( '' )
  print ( '  The matrix A:' )
  print ( A )
#
#  Decompose the matrix.
#
  print ( '' )
  print ( '  Decompose the matrix.' )

  job = 11
  lda = m
  ldu = m
  ldv = n

  A, s, e, u, v, info = dsvdc ( A, lda, m, n, ldu, ldv, job )

  if ( info != 0 ):
    print ( '' )
    print ( 'DSVDC_TEST - Warning:' )
    print ( '  DSVDC returned nonzero INFO = ', info )
#   return

  print ( '' )
  print ( '  Singular values:' )
  print ( s )

  print ( '' )
  print ( '  Left singular vectors U:' )
  print ( u )

  print ( '' )
  print ( '  Right singular vectors V:' )
  print ( v )

  sigma = np.zeros ( [ m, n ] )

  for i in range ( 0, min ( m, n ) ):
    sigma[i,i] = s[i]

  usv = np.dot ( u, np.dot ( sigma, v.transpose ( ) ) )

  print ( '' )
  print ( '  Product U * S * V'' should equal A:' )
  print ( usv )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

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

  return

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