#! /usr/bin/env python
#
def dqrdc ( a, lda, n, p, jpvt, job ):

#*****************************************************************************80
#
## DQRDC computes the QR factorization of a real rectangular matrix.
#
#  Discussion:
#
#    DQRDC 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 GNU LGPL license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    Python 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.
#
#  Discussion:
#
#    DQRDC computes the QR factorization.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from dqrsl import dqrsl

  n = 3
  p = 3
  lda = n

  print ( '' )
  print ( 'DQRDC_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  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 ( '' )

  for i in range ( 0, n ):
    for j in range ( 0, p ):
      print ( '  %14.6g' % ( a[i,j] ), end = '' )
    print ( '' )
#
#  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 ( '' )

  for i in range ( 0, n ):
    for j in range ( 0, p ):
      print ( '  %14.6g' % ( a[i,j] ), end = '' )
    print ( '' )
#
#  ...and in QRAUX.
#
  print ( '' )
  print ( '  The QRAUX vector, containing some additional' )
  print ( '  information defining Q:' )
  print ( '' )

  for i in range ( 0, n ):
    print ( '  %14.6g' % ( qraux[i] ) )
  print ( '' )
#
#  Print out the resulting R factor.
#
  r = np.zeros ( [ n, p ] )

  print ( '' )
  print ( '  The R factor:' )
  print ( '' )

  for i in range ( 0, n ):
    for j in range ( 0, p ):
      if ( i <= j ):
        r[i,j] = a[i,j]
      print ( '  %14.6g' % ( r[i,j] ), end = '' )
    print ( '' )
#
#  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 ( '' )

  for i in range ( 0, n ):
    for j in range ( 0, n ):
      print ( '  %14.6g' % ( q[i,j] ), end = '' )
    print ( '' )
#
#  Compute Q*R to verify that it equals A.
#
  b = np.dot ( q, r )
#
#  Print the result.
#
  print ( '' )
  print ( '  The product Q * R:' )
  print ( '' )

  for i in range ( 0, n ):
    for j in range ( 0, p ):
      print ( '  %14.6g' % ( b[i,j] ), end = '' )
    print ( '' )
#
#  Terminate.
#
  print ( '' )
  print ( 'DQRDC_TEST' )
  print ( '  Normal end of execution.' )
  return

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