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