Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

svd.cpp

Go to the documentation of this file.
00001 //$$svd.cpp                           singular value decomposition
00002 
00003 // Copyright (C) 1991,2,3,4,5: R B Davies
00004 // Updated 17 July, 1995
00005 
00006 #ifndef WANT_MATH
00007 #define WANT_MATH
00008 #endif
00009 
00010 #include "include.h"
00011 #include "newmatap.h"
00012 #include "newmatrm.h"
00013 #include "precisio.h"
00014 
00015 #ifdef use_namespace
00016 namespace NEWMAT {
00017 #endif
00018 
00019 #ifdef DO_REPORT
00020 #define REPORT { static ExeCounter ExeCount(__LINE__,15); ++ExeCount; }
00021 #else
00022 #define REPORT {}
00023 #endif
00024 
00025 
00026 
00027 
00028 void SVD(const Matrix& A, DiagonalMatrix& Q, Matrix& U, Matrix& V,
00029    bool withU, bool withV)
00030 // from Wilkinson and Reinsch: "Handbook of Automatic Computation"
00031 {
00032    REPORT
00033    Tracer trace("SVD");
00034    Real eps = FloatingPointPrecision::Epsilon();
00035    Real tol = FloatingPointPrecision::Minimum()/eps;
00036 
00037    int m = A.Nrows(); int n = A.Ncols();
00038    if (m<n)
00039       Throw(ProgramException("Want no. Rows >= no. Cols", A));
00040    if (withV && &U == &V)
00041       Throw(ProgramException("Need different matrices for U and V", U, V));
00042    U = A; Real g = 0.0; Real f,h; Real x = 0.0; int i;
00043    RowVector E(n); RectMatrixRow EI(E,0); Q.ReSize(n);
00044    RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1);
00045 
00046    if (n) for (i=0;;)
00047    {
00048       EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare();
00049       if (s<tol) { REPORT Q.element(i) = 0.0; }
00050       else
00051       {
00052          REPORT
00053          f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g;
00054          Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i;
00055          while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); }
00056       }
00057 
00058       s = URI.SumSquare();
00059       if (s<tol) { REPORT g = 0.0; }
00060       else
00061       {
00062          REPORT
00063          f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g;
00064          EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i;
00065          while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); }
00066       }
00067 
00068       Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) { REPORT x = y; }
00069       if (++i == n) { REPORT break; }
00070       UCI.DownDiag(); URI.DownDiag();
00071    }
00072 
00073    if (withV)
00074    {
00075       REPORT
00076       V.ReSize(n,n); V = 0.0; RectMatrixCol VCI(V,n-1,n-1,1);
00077       if (n) { VCI.First() = 1.0; g=E.element(n-1); if (n!=1) URI.UpDiag(); }
00078       for (i=n-2; i>=0; i--)
00079       {
00080          VCI.Left();
00081          if (g!=0.0)
00082          {
00083             VCI.Divide(URI, URI.First()*g); int j = n-i;
00084             RectMatrixCol VCJ = VCI;
00085             while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); }
00086          }
00087          VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i);
00088          if (i==0) break;
00089          URI.UpDiag();
00090       }
00091    }
00092 
00093    if (withU)
00094    {
00095       REPORT
00096       for (i=n-1; i>=0; i--)
00097       {
00098          g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero();
00099          if (g!=0.0)
00100          {
00101             h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI;
00102             while (--j)
00103             {
00104                UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ;
00105                UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h);
00106             }
00107             UCI.Divide(g);
00108          }
00109          else UCI.Zero();
00110          UCI.First() += 1.0;
00111          if (i==0) break;
00112          UCI.UpDiag();
00113       }
00114    }
00115 
00116    eps *= x;
00117    for (int k=n-1; k>=0; k--)
00118    {
00119       Real z = -FloatingPointPrecision::Maximum(); // to keep Gnu happy
00120       Real y; int limit = 50; int l = 0;
00121       while (limit--)
00122       {
00123          Real c, s; int i; int l1=k; bool tfc=false;
00124          for (l=k; l>=0; l--)
00125          {
00126 //          if (fabs(E.element(l))<=eps) goto test_f_convergence;
00127             if (fabs(E.element(l))<=eps) { REPORT tfc=true; break; }
00128             if (fabs(Q.element(l-1))<=eps) { REPORT l1=l; break; }
00129             REPORT
00130          }
00131          if (!tfc)
00132          {
00133             REPORT
00134             l=l1; l1=l-1; s = -1.0; c = 0.0;
00135             for (i=l; i<=k; i++)
00136             {
00137                f = - s * E.element(i); E.element(i) *= c;
00138 //             if (fabs(f)<=eps) goto test_f_convergence;
00139                if (fabs(f)<=eps) { REPORT break; }
00140                g = Q.element(i); h = pythag(g,f,c,s); Q.element(i) = h;
00141                if (withU)
00142                {
00143                   REPORT
00144                   RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1);
00145                   ComplexScale(UCJ, UCI, c, s);
00146                }
00147             }
00148          }
00149 //       test_f_convergence: z = Q.element(k); if (l==k) goto convergence;
00150          z = Q.element(k);  if (l==k) { REPORT break; }
00151 
00152          x = Q.element(l); y = Q.element(k-1);
00153          g = E.element(k-1); h = E.element(k);
00154          f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y);
00155          if (f>1)         { REPORT g = f * sqrt(1 + square(1/f)); }
00156          else if (f<-1)   { REPORT g = -f * sqrt(1 + square(1/f)); }
00157          else             { REPORT g = sqrt(f*f + 1); }
00158             { REPORT f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x; }
00159 
00160          c = 1.0; s = 1.0;
00161          for (i=l+1; i<=k; i++)
00162          {
00163             g = E.element(i); y = Q.element(i); h = s*g; g *= c;
00164             z = pythag(f,h,c,s); E.element(i-1) = z;
00165             f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c;
00166             if (withV)
00167             {
00168                REPORT
00169                RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1);
00170                ComplexScale(VCI, VCJ, c, s);
00171             }
00172             z = pythag(f,h,c,s); Q.element(i-1) = z;
00173             f = c*g + s*y; x = -s*g + c*y;
00174             if (withU)
00175             {
00176                REPORT
00177                RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1);
00178                ComplexScale(UCI, UCJ, c, s);
00179             }
00180          }
00181          E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x;
00182       }
00183       if (l!=k) { Throw(ConvergenceException(A)); }
00184 // convergence:
00185       if (z < 0.0)
00186       {
00187          REPORT
00188          Q.element(k) = -z;
00189          if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); }
00190       }
00191    }
00192    if (withU & withV) SortSV(Q, U, V);
00193    else if (withU) SortSV(Q, U);
00194    else if (withV) SortSV(Q, V);
00195    else SortDescending(Q);
00196 }
00197 
00198 void SVD(const Matrix& A, DiagonalMatrix& D)
00199 { REPORT Matrix U; SVD(A, D, U, U, false, false); }
00200 
00201 
00202 
00203 #ifdef use_namespace
00204 }
00205 #endif
00206 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3