function [ a, q ] = minfit ( n, tol, a )

%*****************************************************************************80
%
%% minfit() computes the singular value decomposition of an N by N array.
%
%  Discussion:
%
%    This is an improved version of the EISPACK routine MINFIT
%    restricted to the case M = N and P = 0.
%
%    The singular values of the array A are returned in Q.  A is
%    overwritten with the orthogonal matrix V such that U * diag(Q) = A * V,
%    where U is another orthogonal matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 July 2016
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization with Derivatives,
%    Prentice Hall, 1973,
%    Reprinted by Dover, 2002.
%
%    James Wilkinson, Christian Reinsch,
%    Handbook for Automatic Computation,
%    Volume II, Linear Algebra, Part 2,
%    Springer Verlag, 1971.
%
%    Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, Yasuhiko Ikebe,
%    Virginia Klema, Cleve Moler,
%    Matrix Eigensystem Routines, EISPACK Guide,
%    Lecture Notes in Computer Science, Volume 6,
%    Springer Verlag, 1976,
%    ISBN13: 978-3540075462,
%    LC: QA193.M37.
%
%  Input:
%
%    integer N, the order of the matrix A.
%
%    real TOL, a tolerance which determines when a vector
%    (a column or part of a column of the matrix) may be considered
%    "essentially" equal to zero.
%
%    real A(N,N), an array whose singular value decomposition is desired.  
%
%  Output:
%
%    real A(N,N), the SVD orthogonal matrix factor V.
%
%    real Q(N), the singular values.
%
  kt_max = 30;

  e = zeros ( n, 1 );
  q = zeros ( n, 1 );
%
%  Householder's reduction to bidiagonal form.
%
  if ( n == 1 )
    q(1) = a(1,1);
    a(1,1) = 1.0;
    return
  end

  g = 0.0;
  x = 0.0;

  for i = 1 : n

    e(i) = g;
    l = i + 1;

    s = sum ( a(i:n,i) .^ 2 );

    g = 0.0;

    if ( tol <= s )

      f = a(i,i);

      g = sqrt ( s );
      if ( 0.0 <= f )
        g = - g;
      end

      h = f * g - s;
      a(i,i) = f - g;

      for j = l : n

        f = ( a(i:n,i)' * a(i:n,j) ) / h;

        a(i:n,j) = a(i:n,j) + f * a(i:n,i);

      end

    end

    q(i) = g;

    s = sum ( a(i,l:n) .^ 2 );

    g = 0.0;

    if ( tol <= s )

      if ( i ~= n )
        f = a(i,i+1);
      end

      g = sqrt ( s );
      if ( 0.0 <= f )
        g = - g;
      end

      h = f * g - s;

      if ( i ~= n )

        a(i,i+1) = f - g;
        e(l:n) = a(i,l:n)' / h;

        for j = l : n

          s = a(j,l:n) * a(i,l:n)';

          a(j,l:n) = a(j,l:n) + s * e(l:n)';

        end

      end

    end

    y = abs ( q(i) ) + abs ( e(i) );

    x = max ( x, y );

  end
%
%  Accumulation of right-hand transformations.
%
  a(n,n) = 1.0;
  g = e(n);
  l = n;

  for i = n - 1 : -1 : 1

    if ( g ~= 0.0 )

      h = a(i,i+1) * g;

      a(l:n,i) = a(i,l:n) / h;

      for j = l : n

        s = a(i,l:n) * a(l:n,j);

        a(l:n,j) = a(l:n,j) + s * a(l:n,i);

      end

    end

    a(i,l:n) = 0.0;
    a(l:n,i) = 0.0;
    a(i,i) = 1.0;

    g = e(i);

    l = i;

  end
%
%  Diagonalization of the bidiagonal form.
%
  epsx = eps * x;

  for k = n : -1 : 1

    kt = 0;

    while ( true )

      kt = kt + 1;

      if ( kt_max < kt )
        e(k) = 0.0;
        fprintf ( 1, '\n' );
        fprintf ( 1, 'MINFIT - Fatal error!\n' );
        fprintf ( 1, '  The QR algorithm failed to converge.\n' );
        error ( 'MINFIT - Fatal error!' );
      end

      skip = false;

      for l2 = k : -1 : 1

        l = l2;

        if ( abs ( e(l) ) <= epsx )
          skip = true;
          break
        end

        if ( l ~= 1 )
          if ( abs ( q(l-1) ) <= epsx )
            break
          end
        end

      end
%
%  Cancellation of E(L) if 1 < L.
%
      if ( ~ skip )

        c = 0.0;
        s = 1.0;

        for i = l : k

          f = s * e(i);
          e(i) = c * e(i);
          if ( abs ( f ) <= epsx )
            break
          end
          g = q(i);
%
%  q(i) = h = sqrt(g*g + f*f).
%
          h = r8_hypot ( f, g );

          q(i) = h;

          if ( h == 0.0 )
            g = 1.0;
            h = 1.0;
          end

          c =   g / h;
          s = - f / h;

        end

      end
%
%  Test for convergence for this index K.
%
      z = q(k);

      if ( l == k )
        if ( z < 0.0 )
          q(k) = - z;
          a(1:n,k) = - a(1:n,k);
        end
        break
      end
%
%  Shift from bottom 2*2 minor.
%
      x = q(l);
      y = q(k-1);
      g = e(k-1);
      h = e(k);
      f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) / ( 2.0 * h * y );

      g = r8_hypot ( f, 1.0 );

      if ( f < 0.0 )
        temp = f - g;
      else
        temp = f + g;
      end

      f = ( ( x - z ) * ( x + z ) + h * ( y / temp - h ) ) / x;
%
%  Next QR transformation.
%
      c = 1.0;
      s = 1.0;

      for i = l + 1 : k

        g = e(i);
        y = q(i);
        h = s * g;
        g = g * c;

        z = r8_hypot ( f, h );

        e(i-1) = z;

        if ( z == 0.0 )
          f = 1.0;
          z = 1.0;
        end

        c = f / z;
        s = h / z;
        f =   x * c + g * s;
        g = - x * s + g * c;
        h = y * s;
        y = y * c;

        for j = 1 : n
          x = a(j,i-1);
          z = a(j,i);
          a(j,i-1) = x * c + z * s;
          a(j,i) = - x * s + z * c;
        end

        z = r8_hypot ( f, h );

        q(i-1) = z;

        if ( z == 0.0 )
          f = 1.0;
          z = 1.0;
        end

        c = f / z;
        s = h / z;
        f =   c * g + s * y;
        x = - s * g + c * y;

      end

      e(l) = 0.0;
      e(k) = f;
      q(k) = x;

    end

  end

  return
end
