source: ntrip/branches/BNC_2.12/newmat/svd.cpp@ 8142

Last change on this file since 8142 was 2013, checked in by mervart, 13 years ago

* empty log message *

File size: 5.9 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file svd.cpp
5/// Singular value decomposition.
6
7// Copyright (C) 1991,2,3,4,5: R B Davies
8// Updated 17 July, 1995
9
10#define WANT_MATH
11
12#include "include.h"
13#include "newmatap.h"
14#include "newmatrm.h"
15#include "precisio.h"
16
17#ifdef use_namespace
18namespace NEWMAT {
19#endif
20
21#ifdef DO_REPORT
22#define REPORT { static ExeCounter ExeCount(__LINE__,15); ++ExeCount; }
23#else
24#define REPORT {}
25#endif
26
27
28
29
30void SVD(const Matrix& A, DiagonalMatrix& Q, Matrix& U, Matrix& V,
31 bool withU, bool withV)
32// from Wilkinson and Reinsch: "Handbook of Automatic Computation"
33{
34 REPORT
35 Tracer trace("SVD");
36 Real eps = FloatingPointPrecision::Epsilon();
37 Real tol = FloatingPointPrecision::Minimum()/eps;
38
39 int m = A.Nrows(); int n = A.Ncols();
40 if (m<n)
41 Throw(ProgramException("Want no. Rows >= no. Cols", A));
42 if (withV && &U == &V)
43 Throw(ProgramException("Need different matrices for U and V", U, V));
44 U = A; Real g = 0.0; Real f,h; Real x = 0.0; int i;
45 RowVector E(n); RectMatrixRow EI(E,0); Q.ReSize(n);
46 RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1);
47
48 if (n) for (i=0;;)
49 {
50 EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare();
51 if (s<tol) { REPORT Q.element(i) = 0.0; }
52 else
53 {
54 REPORT
55 f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g;
56 Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i;
57 while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); }
58 }
59
60 s = URI.SumSquare();
61 if (s<tol) { REPORT g = 0.0; }
62 else
63 {
64 REPORT
65 f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g;
66 EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i;
67 while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); }
68 }
69
70 Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) { REPORT x = y; }
71 if (++i == n) { REPORT break; }
72 UCI.DownDiag(); URI.DownDiag();
73 }
74
75 if (withV)
76 {
77 REPORT
78 V.ReSize(n,n); V = 0.0; RectMatrixCol VCI(V,n-1,n-1,1);
79 if (n) { VCI.First() = 1.0; g=E.element(n-1); if (n!=1) URI.UpDiag(); }
80 for (i=n-2; i>=0; i--)
81 {
82 VCI.Left();
83 if (g!=0.0)
84 {
85 VCI.Divide(URI, URI.First()*g); int j = n-i;
86 RectMatrixCol VCJ = VCI;
87 while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); }
88 }
89 VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i);
90 if (i==0) break;
91 URI.UpDiag();
92 }
93 }
94
95 if (withU)
96 {
97 REPORT
98 for (i=n-1; i>=0; i--)
99 {
100 g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero();
101 if (g!=0.0)
102 {
103 h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI;
104 while (--j)
105 {
106 UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ;
107 UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h);
108 }
109 UCI.Divide(g);
110 }
111 else UCI.Zero();
112 UCI.First() += 1.0;
113 if (i==0) break;
114 UCI.UpDiag();
115 }
116 }
117
118 eps *= x;
119 for (int k=n-1; k>=0; k--)
120 {
121 Real z = -FloatingPointPrecision::Maximum(); // to keep Gnu happy
122 Real y; int limit = 50; int l = 0;
123 while (limit--)
124 {
125 Real c, s; int i; int l1=k; bool tfc=false;
126 for (l=k; l>=0; l--)
127 {
128// if (fabs(E.element(l))<=eps) goto test_f_convergence;
129 if (fabs(E.element(l))<=eps) { REPORT tfc=true; break; }
130 if (fabs(Q.element(l-1))<=eps) { REPORT l1=l; break; }
131 REPORT
132 }
133 if (!tfc)
134 {
135 REPORT
136 l=l1; l1=l-1; s = -1.0; c = 0.0;
137 for (i=l; i<=k; i++)
138 {
139 f = - s * E.element(i); E.element(i) *= c;
140// if (fabs(f)<=eps) goto test_f_convergence;
141 if (fabs(f)<=eps) { REPORT break; }
142 g = Q.element(i); h = pythag(g,f,c,s); Q.element(i) = h;
143 if (withU)
144 {
145 REPORT
146 RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1);
147 ComplexScale(UCJ, UCI, c, s);
148 }
149 }
150 }
151// test_f_convergence: z = Q.element(k); if (l==k) goto convergence;
152 z = Q.element(k); if (l==k) { REPORT break; }
153
154 x = Q.element(l); y = Q.element(k-1);
155 g = E.element(k-1); h = E.element(k);
156 f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y);
157 if (f>1) { REPORT g = f * sqrt(1 + square(1/f)); }
158 else if (f<-1) { REPORT g = -f * sqrt(1 + square(1/f)); }
159 else { REPORT g = sqrt(f*f + 1); }
160 { REPORT f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x; }
161
162 c = 1.0; s = 1.0;
163 for (i=l+1; i<=k; i++)
164 {
165 g = E.element(i); y = Q.element(i); h = s*g; g *= c;
166 z = pythag(f,h,c,s); E.element(i-1) = z;
167 f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c;
168 if (withV)
169 {
170 REPORT
171 RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1);
172 ComplexScale(VCI, VCJ, c, s);
173 }
174 z = pythag(f,h,c,s); Q.element(i-1) = z;
175 f = c*g + s*y; x = -s*g + c*y;
176 if (withU)
177 {
178 REPORT
179 RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1);
180 ComplexScale(UCI, UCJ, c, s);
181 }
182 }
183 E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x;
184 }
185 if (l!=k) { Throw(ConvergenceException(A)); }
186// convergence:
187 if (z < 0.0)
188 {
189 REPORT
190 Q.element(k) = -z;
191 if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); }
192 }
193 }
194 if (withU & withV) SortSV(Q, U, V);
195 else if (withU) SortSV(Q, U);
196 else if (withV) SortSV(Q, V);
197 else sort_descending(Q);
198}
199
200void SVD(const Matrix& A, DiagonalMatrix& D)
201{ REPORT Matrix U; SVD(A, D, U, U, false, false); }
202
203
204
205#ifdef use_namespace
206}
207#endif
208
209///@}
Note: See TracBrowser for help on using the repository browser.