function [ a, s, e, u, v, info ] = 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:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB return 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 ( kind = 8 ) 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.
%
  maxit = 30;
  s = [];
  e = [];
  u = [];
  v = [];
%
%  Determine what is to be computed.
%
  wantu = 0;
  wantv = 0;
  jobu = floor ( mod ( job, 100 ) / 10 );

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

  if ( jobu ~= 0 )
    wantu = 1;
  end

  if ( mod ( job, 10 ) ~= 0 )
    wantv = 1;
  end
%
%  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 = 1 : lu
%
%  Compute the transformation for the L-th column and
%  place the L-th diagonal in S(L).
%
    if ( l <= nct )

      s(l) = dnrm2 ( m-l+1, a(l:m,l), 1 );

      if ( s(l) ~= 0.0 )
        if ( a(l,l) ~= 0.0 )
          s(l) = abs ( s(l) ) * r8_sign ( a(l,l) );
        end
        a(l:m,l) = a(l:m,l) / s(l);
        a(l,l) = 1.0 + a(l,l);
      end

      s(l) = -s(l);

    end

    for j = l+1 : n
%
%  Apply the transformation.
%
      if ( l <= nct && s(l) ~= 0.0 )
        t = -ddot ( m-l+1, a(l:m,l), 1, a(l:m,j), 1 ) / a(l,l);
        a(l:m,j) = daxpy ( m-l+1, t, a(l:m,l), 1, a(l:m,j), 1 );
      end
%
%  Place the L-th row of A into E for the
%  subsequent calculation of the row transformation.
%
      e(j) = a(l,j);

    end
%
%  Place the transformation in U for subsequent back multiplication.
%
    if ( wantu && l <= nct )
      u(l:m,l) = a(l:m,l);
    end
%
%  Compute the L-th row transformation and place the
%  L-th superdiagonal in E(L).
%
    if ( l <= nrt )

      e(l) = dnrm2 ( n-l, e(l+1:n), 1 );

      if ( e(l) ~= 0.0 )
        if ( e(l+1) ~= 0.0 )
          e(l) = abs ( e(l) ) * r8_sign ( e(l+1) );
        end
        e(l+1:n) = e(l+1:n) / e(l);
        e(l+1) = 1.0 + e(l+1);
      end

      e(l) = -e(l);
%
%  Apply the transformation.
%
      if ( l + 1 <= m && e(l) ~= 0.0 )

        work(l+1:m) = 0.0;

        for j = l+1 : n
          work(l+1:m) = daxpy ( m-l, e(j), a(l+1:m,j)', 1, work(l+1:m), 1 );
        end

        for j = l+1 : n
          a(l+1:m,j) = ...
            daxpy ( m-l, -e(j)/e(l+1), work(l+1:m)', 1, a(l+1:m,j), 1 );
        end

      end
%
%  Place the transformation in V for subsequent back multiplication.
%
      if ( wantv )
        v(l+1:n,l) = e(l+1:n);
      end

    end

  end
%
%  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) = a(nctp1,nctp1);
  end

  if ( m < mn )
    s(mn) = 0.0;
  end

  if ( nrtp1 < mn )
    e(nrtp1) = a(nrtp1,mn);
  end

  e(mn) = 0.0;
%
%  If required, generate U.
%
  if ( wantu )

    u(1:m,nctp1:ncu) = 0.0;

    for j = nctp1 : ncu
      u(j,j) = 1.0;
    end

    for ll = 1 : nct

      l = nct - ll + 1;

      if ( s(l) ~= 0.0 )

        for j = l+1 : ncu
          t = -ddot ( m-l+1, u(l:m,l), 1, u(l:m,j), 1 ) / u(l,l);
          u(l:m,j) = daxpy ( m-l+1, t, u(l:m,l), 1, u(l:m,j), 1 );
        end

        u(l:m,l) = -u(l:m,l);
        u(l,l) = 1.0 + u(l,l);
        u(1:l-1,l) = 0.0;

      else

        u(1:m,l) = 0.0;
        u(l,l) = 1.0;

      end

    end

  end
%
%  If it is required, generate V.
%
  if ( wantv )

    for ll = 1 : n

      l = n - ll + 1;

      if ( l <= nrt && e(l) ~= 0.0 )

        for j = l+1 : n
          t = -ddot ( n-l, v(l+1:n,l), 1, v(l+1:n,j), 1 ) / v(l+1,l);
          v(l+1:n,j) = daxpy ( n-l, t, v(l+1:n,l), 1, v(l+1:n,j), 1 );
        end

      end

      v(1:n,l) = 0.0;
      v(l,l) = 1.0;

    end

  end
%
%  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 )
      info = mn;
      return
    end
%
%  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 ll = 1 : mn

      l = mn - ll;

      if ( l == 0 )
        break
      end

      test = abs ( s(l) ) + abs ( s(l+1) );
      ztest = test + abs ( e(l) );

      if ( ztest == test )
        e(l) = 0.0;
        break
      end

    end

    if ( l == mn - 1 )

      kase = 4;

    else

      mp1 = mn + 1;

      for lls = l+1 : mn+1

        ls = mn - lls + l + 1;

        if ( ls == l )
          break
        end

        test = 0.0;
        if ( ls ~= mn )
          test = test + abs ( e(ls) );
        end

        if ( ls ~= l + 1 )
          test = test + abs ( e(ls-1) );
        end

        ztest = test + abs ( s(ls) );

        if ( ztest == test )
          s(ls) = 0.0;
          break
        end

      end

      if ( ls == l )
        kase = 3;
      elseif ( ls == mn )
        kase = 1;
      else
        kase = 2;
        l = ls;
      end

    end

    l = l + 1;
%
%  Deflate negligible S(MN).
%
    if ( kase == 1 )

      mm1 = mn - 1;
      f = e(mn-1);
      e(mn-1) = 0.0;

      for kk = l : mm1

        k = mm1 - kk + l;
        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);
        end

        if ( wantv )
          [ v(1:n,k), v(1:n,mn) ] = ...
            drot ( n, v(1:n,k), 1, v(1:n,mn), 1, cs, sn );
        end

      end
%
%  Split at negligible S(L).
%
    elseif ( kase == 2 )

      f = e(l-1);
      e(l-1) = 0.0;

      for k = 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 )
          [ u(1:m,k), u(1:m,l-1) ] = ...
            drot ( m, u(1:m,k), 1, u(1:m,l-1), 1, cs, sn );
        end

      end
%
%  Perform one QR step.
%
    elseif ( kase == 3 )
%
%  Calculate the shift.
%
      scale = max ( abs ( s(mn)   ), ...
             max ( abs ( s(mn-1) ), ...
             max ( abs ( e(mn-1) ), ...
             max ( abs ( s(l)    ), abs ( e(l) ) ) ) ) );

      sm = s(mn) / scale;
      smm1 = s(mn-1) / scale;
      emm1 = e(mn-1) / 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 || c ~= 0.0 )
        shift = sqrt ( b * b + c );
        if ( b < 0.0 )
          shift = -shift;
        end
        shift = c / ( b + shift );
      end

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

      for k = l : mm1

        [ cs, sn, f, g ] = drotg ( f, g );

        if ( k ~= l )
          e(k-1) = f;
        end

        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 )
          [ v(1:n,k), v(1:n,k+1) ] = ...
            drot ( n, v(1:n,k), 1, v(1:n,k+1), 1, cs, sn );
        end

        [ 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 && k < m )
          [ u(1:m,k), u(1:m,k+1) ] = ...
            drot ( m, u(1:m,k), 1, u(1:m,k+1), 1, cs, sn );
        end

      end

      e(mn-1) = f;
      iter = iter + 1;
%
%  Convergence.
%
    elseif ( kase == 4 )
%
%  Make the singular value nonnegative.
%
      if ( s(l) < 0.0 )
        s(l) = -s(l);
        if ( wantv )
          v(1:n,l) = -v(1:n,l);
        end
      end
%
%  Order the singular value.
%
      while ( 1 )

        if ( l == mm )
          break
        end

        if ( s(l+1) <= s(l) )
          break
        end

        t = s(l);
        s(l) = s(l+1);
        s(l+1) = t;

        if ( wantv && l < n )
          temp(1:n) = v(1:n,l);
          v(1:n,l) = v(1:n,l+1);
          v(1:n,l+1) = temp(1:n);
        end

        if ( wantu && l < m )
          temp(1:m) = u(1:m,l);
          u(1:m,l) = u(1:m,l+1);
          u(1:m,l+1) = temp(1:m);
        end

        l = l + 1;

      end

      iter = 0;
      mn = mn - 1;

    end

  end

  return
end
