#! /usr/bin/env python
#
def dqrlss ( a, lda, m, n, kr, b, jpvt, qraux ):

#*****************************************************************************80
#
## DQRLSS solves a linear system in a least squares sense.
#
#  Discussion:
#
#    DQRLSS must be preceded by a call to DQRANK.
#
#    The system is to be solved is
#      A * X = B
#    where
#      A is an M by N matrix with rank KR, as determined by DQRANK,
#      B is a given M-vector,
#      X is the N-vector to be computed.
#
#    A solution X, with at most KR nonzero components, is found which
#    minimizes the 2-norm of the residual (A*X-B).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    Python version by John Burkardt.
#
#  Parameters:
#
#    Input, real A(LDA,N), the QR factorization information
#    from DQRANK.  The triangular matrix R of the QR factorization is
#    contained in the upper triangle and information needed to recover
#    the orthogonal matrix Q is stored below the diagonal in A and in
#    the vector QRAUX.
#
#    Input, integer LDA, the leading dimension of A, which must
#    be at least M.
#
#    Input, integer M, the number of rows of A.
#
#    Input, integer N, the number of columns of A.
#
#    Input, integer KR, the rank of the matrix, as estimated
#    by DQRANK.
#
#    Input, real B(M), the right hand side of the linear system.
#
#    Output, real X(N), a least squares solution to the
#    linear system.
#
#    Output, real R(M), the residual, B - A*X.  RSD may
#    overwite B.
#
#    Input, integer JPVT(N), the pivot information from DQRANK.
#    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
#    independent to within the tolerance TOL and the remaining columns
#    are linearly dependent.
#
#    Input, real QRAUX(N), auxiliary information from DQRANK
#    defining the QR factorization.
#
  import numpy as np
  from dqrsl import dqrsl
#
#  Solve the reduced system of rank KR.
#
  if ( 0 < kr ):
    job = 110
    qy, qty, x, r, ab, info = dqrsl ( a, lda, m, kr, qraux, b, job );
#
#  Reverse the pivoting to recover the original variable ordering.
#
  jpvt[0:n] = - jpvt[0:n]

  for j in range ( kr, n ):
    x = np.append ( x, 0.0 )

  for j in range ( 0, n ):

    if ( jpvt[j] <= 0 ):

      k = - jpvt[j]
      jpvt[j] = k

      while ( k != j + 1 ):

        t = x[j]
        x[j] = x[k-1]
        x[k-1] = t

        jpvt[k-1] = - jpvt[k-1]

        k = jpvt[k-1]

  return x, r

def dqrlss_test ( ):

#*****************************************************************************80
#
## DQRLSS_TEST tests DQRLSS.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from dqrank import dqrank
  from r8vec_print import r8vec_print

  print ( '' )
  print ( 'DQRLSS_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  DQRLSS solves a rectangular linear system A*x = b' )
  print ( '  with possibly non-full rank, in the least squares sense.' )
  print ( '  Compare a tabulated solution X1 to the computed result X2.' )
#
#  Set the matrix A.
#
  m = 6;
  n = 6;
  a = np.zeros ( [ m, n ] )
  b = np.zeros ( m )
  x1 = np.zeros ( n )
  
  for j in range ( 0, n ):
    x1[j] = float ( j + 1 )
    for i in range ( 0, m ):
      a[i,j] = 1.0 / float ( i + j + 1 )
      b[i] = b[i] + a[i,j] * x1[j]

  print ( '' )
  print ( '  The matrix A:' )
  print ( '' )

  for i in range ( 0, m ):
    for j in range ( 0, n ):
      print ( '  %14.6g' % ( a[i,j] ) ),
    print ( '' )

  r8vec_print ( n, x1, '  Solution x1:' )
  r1 = b - np.dot ( a, x1 )
  r8vec_print ( m, r1, '  Residual b-A*x1:' )

  lda = m
  tol = 0.000001
#
#  Factor the matrix.
#
  a_qr = a.copy ( )
  kr, jpvt, qraux, a_qr = dqrank ( a_qr, lda, m, n, tol )

  print ( '' )
  print ( '  Matrix is of order m=%d by n=%d' % ( m, n ) )
  print ( '  Condition number tolerance is %g' % ( tol ) )
  print ( '  DQRANK estimates rank as %d' % ( kr ) )
#
#  Solve the linear least squares problem.
#
  x2, r = dqrlss ( a_qr, lda, m, n, kr, b, jpvt, qraux )

  r8vec_print ( n, x2, '  Computed solution x2:' )
  r2 = b - np.dot ( a, x2 )
  r8vec_print ( m, r2, '  Residual b-A*x2:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'DQRLSS_TEST' )
  print ( '  Normal end of execution.' )
  return

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