Bug in Singular Value Decomposition

["Brian Shier" <abcshier@hotmail.com>]

In an earlier message about the JAMA SVD code I mentioned that the need to 
initialize the elements of U to zero in the constructor is suspicious. I 
believe that the need for this is removed by the change in the following 
code snippet.

  // If required, generate U.

  if (wantu)
    {
	for (j = nct; j < nu; j++)
      {
      for (i = 0; i < m; i++)
        U[i][j] = 0.0;

      U[j][j] = 1.0;
      }

    for (k = nct-1; k >= 0; k--)
      {
      if (s[k] != 0.0)
        {
	    for (j = k+1; j < nu; j++)
          {
          t = 0.0;

          for (i = k; i < m; i++)
            t += U[i][k]*U[i][j];

          t = -t/U[k][k];

          for (i = k; i < m; i++)
            U[i][j] += t*U[i][k];
          }

        for (i = k; i < m; i++ )
          U[i][k] = -U[i][k];

        U[k][k] = 1.0 + U[k][k];

//      for (i = 0; i < k-1; i++)       // Error in JAMA version ?
        for (i = 0; i < k; i++)         // Correction
          U[i][k] = 0.0;
        }
      else
        {
        for (i = 0; i < m; i++)
          U[i][k] = 0.0;

        U[k][k] = 1.0;
        }
      }
    }

Please investigate.

-- Brian Shier

_________________________________________________________________
The new MSN 8: smart spam protection and 2 months FREE*  
http://join.msn.com/?page=features/junkmail