home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / SVD.CXX < prev    next >
C/C++ Source or Header  |  1992-11-30  |  5KB  |  161 lines
  1. //$$svd.cxx                           singular value decomposition
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8. #include "newmat.h"
  9. #include "newmatrm.h"
  10. #include "precisio.h"
  11.  
  12.  
  13. void SVD(const Matrix& A, DiagonalMatrix& Q, Matrix& U, Matrix& V,
  14.    Boolean withU, Boolean withV)
  15. // from Wilkinson and Reinsch: "Handbook of Automatic Computation"
  16. {
  17.    Tracer trace("SVD");
  18.    Real eps = FloatingPointPrecision::Epsilon();
  19.    Real tol = FloatingPointPrecision::Minimum()/eps;
  20.  
  21.    int m = A.Nrows(); int n = A.Ncols();
  22.    if (m<n) 
  23.       Throw(ProgramException("Want no. Rows >= no. Cols", A));
  24.  
  25.    U = A; Real g = 0.0; Real f,h; Real x = 0.0;
  26.    RowVector E(n); RectMatrixRow EI(E,0); Q.ReDimension(n);
  27.    RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1);
  28.  
  29.    for (int i=0; i<n; i++)
  30.    {
  31.       EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare();
  32.       if (s<tol) Q.element(i) = 0.0;
  33.       else
  34.       {
  35.      f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g;
  36.      Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i;
  37.      while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); }
  38.       }
  39.  
  40.       s = URI.SumSquare();
  41.       if (s<tol) g = 0.0;
  42.       else
  43.       {
  44.      f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g;
  45.      EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i;
  46.      while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); }
  47.       }
  48.  
  49.       Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) x = y;
  50.       UCI.DownDiag(); URI.DownDiag();
  51.    }
  52.  
  53.    if (withV)
  54.    {
  55.       V.ReDimension(n,n); V = 0.0; RectMatrixCol VCI(V,n,n,0);
  56.       for (i=n-1; i>=0; i--)
  57.       {
  58.      URI.UpDiag(); VCI.Left();
  59.      if (g!=0.0)
  60.      {
  61.         VCI.Divide(URI, URI.First()*g); int j = n-i;
  62.         RectMatrixCol VCJ = VCI;
  63.         while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); }
  64.      }
  65.      VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i);
  66.       }
  67.    }
  68.  
  69.    if (withU)
  70.    {
  71.       for (i=n-1; i>=0; i--)
  72.       {
  73.      UCI.UpDiag(); g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero();
  74.          if (g!=0.0)
  75.          {
  76.         h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI;
  77.         while (--j)
  78.             {
  79.            UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ;
  80.                UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h);
  81.             }
  82.         UCI.Divide(g);
  83.          }
  84.          else UCI.Zero();
  85.          UCI.First() += 1.0;
  86.       }
  87.    }
  88.  
  89.    eps *= x;
  90.    for (int k=n-1; k>=0; k--)
  91.    {
  92.       Real z; Real y; int limit = 50; int l;
  93.       while (limit--)
  94.       {
  95.      Real c=0.0; Real s=1.0; int i; int l1=k; Boolean tfc=FALSE;
  96.      for (l=k; l>=0; l--)
  97.      {
  98. //          if (fabs(E.element(l))<=eps) goto test_f_convergence;
  99.         if (fabs(E.element(l))<=eps) { tfc=TRUE; break; }
  100.         if (fabs(Q.element(l-1))<=eps) { l1=l; break; }
  101.      }
  102.          if (!tfc)
  103.          {
  104.         l=l1; l1=l-1;
  105.         for (i=l; i<=k; i++)
  106.         {
  107.            f = s * E.element(i); E.element(i) *= c;
  108. //             if (fabs(f)<=eps) goto test_f_convergence;
  109.            if (fabs(f)<=eps) break;
  110.            g = Q.element(i); h = sqrt(f*f + g*g); Q.element(i) = h;
  111.            if (withU)
  112.            {
  113.               RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1);
  114.               ComplexScale(UCI, UCJ, g/h, -f/h);
  115.                }
  116.             }
  117.      }
  118. //    test_f_convergence: z = Q.element(k); if (l==k) goto convergence;
  119.      z = Q.element(k);  if (l==k) break;
  120.  
  121.      x = Q.element(l); y = Q.element(k-1);
  122.      g = E.element(k-1); h = E.element(k);
  123.      f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y); g = sqrt(f*f + 1);
  124.      f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x;
  125.  
  126.      c = 1.0; s = 1.0;
  127.      for (i=l+1; i<=k; i++)
  128.      {
  129.         g = E.element(i); y = Q.element(i); h = s*g; g *= c;
  130.         z = sqrt(f*f + h*h); E.element(i-1) = z; c = f/z; s = h/z;
  131.         f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c;
  132.         if (withV)
  133.         {
  134.            RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1);
  135.            ComplexScale(VCI, VCJ, c, s);
  136.         }
  137.         z = sqrt(f*f + h*h); Q.element(i-1) = z; c = f/z; s = h/z;
  138.         f = c*g + s*y; x = -s*g + c*y;
  139.         if (withU)
  140.         {
  141.            RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1);
  142.            ComplexScale(UCI, UCJ, c, s);
  143.         }
  144.      }
  145.      E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x;
  146.       }
  147.       if (l!=k) { Throw(ConvergenceException(A)); }
  148. // convergence:
  149.       if (z < 0.0)
  150.       {
  151.      Q.element(k) = -z;
  152.      if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); }
  153.       }
  154.    }
  155. }
  156.  
  157. void SVD(const Matrix& A, DiagonalMatrix& D)
  158. { Matrix U; SVD(A, D, U, U, FALSE, FALSE); }
  159.  
  160.  
  161.