svd.cpp
Go to the documentation of this file.
1 
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
18 namespace 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 
30 void 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");
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 
200 void 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 
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
void SVD(const Matrix &A, DiagonalMatrix &Q, Matrix &U, Matrix &V, bool withU, bool withV)
Definition: svd.cpp:30
static Real Maximum()
Definition: precisio.h:235
#define REPORT
Definition: svd.cpp:24
Real & element(int, int)
Definition: newmat6.cpp:748
virtual void ReSize(int m, int n)
Definition: newmat.h:662
double Real
Definition: include.h:307
Real SumSquare() const
Definition: newmatrm.cpp:55
Real pythag(Real f, Real g, Real &c, Real &s)
Definition: newmatrm.cpp:189
void Up()
Definition: newmatrm.h:80
Real square(Real x)
Definition: bandmat.cpp:672
void ReSize(int m)
Definition: newmat.h:935
int Nrows() const
Definition: newmat.h:494
void Reset(const Matrix &, int, int, int)
Definition: newmatrm.cpp:28
static Real Minimum()
Definition: precisio.h:241
void ComplexScale(RectMatrixCol &U, RectMatrixCol &V, Real x, Real y)
Definition: newmatrm.cpp:135
void sort_descending(GeneralMatrix &)
Definition: sort.cpp:49
void Left()
Definition: newmatrm.h:81
void Divide(const RectMatrixRowCol &, Real)
Definition: newmatrm.cpp:97
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
void UpDiag()
Definition: newmatrm.h:46
Real & First()
Definition: newmatrm.h:44
int Ncols() const
Definition: newmat.h:495
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
Diagonal matrix.
Definition: newmat.h:896
void DownDiag()
Definition: newmatrm.h:45
static Real Epsilon()
Definition: precisio.h:231
Covergence failure exception.
Definition: newmat.h:1922
void AddScaled(const RectMatrixRowCol &, Real)
Definition: newmatrm.cpp:82
Row vector.
Definition: newmat.h:953
void Right()
Definition: newmatrm.h:63
void Down()
Definition: newmatrm.h:62
void SortSV(DiagonalMatrix &D, Matrix &U, bool ascending=false)
Definition: sort.cpp:194
void Right()
Definition: newmatrm.h:79
Real & element(int)
Definition: newmat6.cpp:792
void Down()
Definition: newmatrm.h:78


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:17