#! /usr/bin/env python
#
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 GNU LGPL license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
#    Python 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
  from drot import drot
  from drotg import drotg

  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 ( ):

#*****************************************************************************80
#
## DSVDC_TEST tests DSVDC.
#
#  Discussion:
#
#    DSVDC computes the singular value decomposition:
#
#      A = U * S * V'
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    07 September 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8mat_print import r8mat_print
  from r8mat_uniform_01 import r8mat_uniform_01
  from r8vec_print import r8vec_print

  m = 6
  n = 4

  print ( '' )
  print ( 'DSVDC_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  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.
#
  seed = 123456789

  a, seed = r8mat_uniform_01 ( m, n, seed )

  r8mat_print ( m, n, a, '  The matrix 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 = %d' % ( info ) )
#   return

  r8vec_print ( min ( m, n ), s, '  Singular values S:' )

  r8mat_print ( m, m, u, '  Left Singular Vector Matrix U:' )

  r8mat_print ( n, n, v, '  Right Singular Vector Matrix 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 ( ) ) )

  r8mat_print ( m, n, usv, '  The product U * S * V\' (should equal A):' )
#
#  Terminate.
#
  print ( '' )
  print ( 'DSVDC_TEST' )
  print ( '  Normal end of execution.' )
  return

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