hholder.cpp
Go to the documentation of this file.
1 
8 
9 // Copyright (C) 1991,2,3,4: R B Davies
10 
11 #define WANT_MATH
12 //#define WANT_STREAM
13 
14 #include "include.h"
15 
16 #include "newmatap.h"
17 
18 #ifdef use_namespace
19 namespace NEWMAT {
20 #endif
21 
22 #ifdef DO_REPORT
23 #define REPORT { static ExeCounter ExeCount(__LINE__,16); ++ExeCount; }
24 #else
25 #define REPORT {}
26 #endif
27 
28 
29 /*************************** QR decompositions ***************************/
30 
31 inline Real square(Real x) { return x*x; }
32 
34 {
35  REPORT
36  Tracer et("QRZT(1)");
37  int n = X.Ncols(); int s = X.Nrows(); L.resize(s);
38  if (n == 0 || s == 0) { L = 0.0; return; }
39  Real* xi = X.Store(); int k;
40  for (int i=0; i<s; i++)
41  {
42  Real sum = 0.0;
43  Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
44  sum = sqrt(sum);
45  if (sum == 0.0)
46  {
47  REPORT
48  k=n; while(k--) { *xi0++ = 0.0; }
49  for (int j=i; j<s; j++) L.element(j,i) = 0.0;
50  }
51  else
52  {
53  L.element(i,i) = sum;
54  Real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
55  for (int j=i+1; j<s; j++)
56  {
57  sum=0.0;
58  xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
59  xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
60  L.element(j,i) = sum;
61  }
62  }
63  }
64 }
65 
66 void QRZT(const Matrix& X, Matrix& Y, Matrix& M)
67 {
68  REPORT
69  Tracer et("QRZT(2)");
70  int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
71  if (Y.Ncols() != n)
72  { Throw(ProgramException("Unequal row lengths",X,Y)); }
73  M.resize(t,s);
74  Real* xi = X.Store(); int k;
75  for (int i=0; i<s; i++)
76  {
77  Real* xj0 = Y.Store(); Real* xi0 = xi;
78  for (int j=0; j<t; j++)
79  {
80  Real sum=0.0;
81  xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
82  xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
83  M.element(j,i) = sum;
84  }
85  }
86 }
87 
88 /*
89 void QRZ(Matrix& X, UpperTriangularMatrix& U)
90 {
91  Tracer et("QRZ(1)");
92  int n = X.Nrows(); int s = X.Ncols(); U.resize(s);
93  Real* xi0 = X.Store(); int k;
94  for (int i=0; i<s; i++)
95  {
96  Real sum = 0.0;
97  Real* xi = xi0; k=n; while(k--) { sum += square(*xi); xi+=s; }
98  sum = sqrt(sum);
99  U.element(i,i) = sum;
100  if (sum==0.0) Throw(SingularException(U));
101  Real* xj0=xi0; k=n; while(k--) { *xj0 /= sum; xj0+=s; }
102  xj0 = xi0;
103  for (int j=i+1; j<s; j++)
104  {
105  sum=0.0;
106  xi=xi0; k=n; xj0++; Real* xj=xj0;
107  while(k--) { sum += *xi * *xj; xi+=s; xj+=s; }
108  xi=xi0; k=n; xj=xj0;
109  while(k--) { *xj -= sum * *xi; xj+=s; xi+=s; }
110  U.element(i,j) = sum;
111  }
112  xi0++;
113  }
114 }
115 */
116 
118 {
119  REPORT
120  Tracer et("QRZ(1)");
121  int n = X.Nrows(); int s = X.Ncols(); U.resize(s); U = 0.0;
122  if (n == 0 || s == 0) return;
123  Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
124  int j, k; int J = s; int i = s;
125  while (i--)
126  {
127  Real* xj0 = xi0; Real* xi = xi0; k = n;
128  if (k) for (;;)
129  {
130  u = u0; Real Xi = *xi; Real* xj = xj0;
131  j = J; while(j--) *u++ += Xi * *xj++;
132  if (!(--k)) break;
133  xi += s; xj0 += s;
134  }
135 
136  Real sum = sqrt(*u0); *u0 = sum; u = u0+1;
137  if (sum == 0.0)
138  {
139  REPORT
140  j = J - 1; while(j--) *u++ = 0.0;
141 
142  xj0 = xi0++; k = n;
143  if (k) for (;;)
144  {
145  *xj0 = 0.0;
146  if (!(--k)) break;
147  xj0 += s;
148  }
149  u0 += J--;
150  }
151  else
152  {
153  int J1 = J-1; j = J1; while(j--) *u++ /= sum;
154 
155  xj0 = xi0; xi = xi0++; k = n;
156  if (k) for (;;)
157  {
158  u = u0+1; Real Xi = *xi; Real* xj = xj0;
159  Xi /= sum; *xj++ = Xi;
160  j = J1; while(j--) *xj++ -= *u++ * Xi;
161  if (!(--k)) break;
162  xi += s; xj0 += s;
163  }
164  u0 += J--;
165  }
166  }
167 }
168 
169 void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
170 {
171  REPORT
172  Tracer et("QRZ(2)");
173  int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
174  if (Y.Nrows() != n)
175  { Throw(ProgramException("Unequal column lengths",X,Y)); }
176  M.resize(s,t); M = 0;Real* m0 = M.Store(); Real* m;
177  Real* xi0 = X.Store();
178  int j, k; int i = s;
179  while (i--)
180  {
181  Real* xj0 = Y.Store(); Real* xi = xi0; k = n;
182  if (k) for (;;)
183  {
184  m = m0; Real Xi = *xi; Real* xj = xj0;
185  j = t; while(j--) *m++ += Xi * *xj++;
186  if (!(--k)) break;
187  xi += s; xj0 += t;
188  }
189 
190  xj0 = Y.Store(); xi = xi0++; k = n;
191  if (k) for (;;)
192  {
193  m = m0; Real Xi = *xi; Real* xj = xj0;
194  j = t; while(j--) *xj++ -= *m++ * Xi;
195  if (!(--k)) break;
196  xi += s; xj0 += t;
197  }
198  m0 += t;
199  }
200 }
201 
202 /*
203 
204 void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
205 {
206  Tracer et("QRZ(2)");
207  int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
208  if (Y.Nrows() != n)
209  { Throw(ProgramException("Unequal column lengths",X,Y)); }
210  M.resize(s,t);
211  Real* xi0 = X.Store(); int k;
212  for (int i=0; i<s; i++)
213  {
214  Real* xj0 = Y.Store();
215  for (int j=0; j<t; j++)
216  {
217  Real sum=0.0;
218  Real* xi=xi0; Real* xj=xj0; k=n;
219  while(k--) { sum += *xi * *xj; xi+=s; xj+=t; }
220  xi=xi0; k=n; xj=xj0++;
221  while(k--) { *xj -= sum * *xi; xj+=t; xi+=s; }
222  M.element(i,j) = sum;
223  }
224  xi0++;
225  }
226 }
227 */
228 
230 {
231  REPORT
232  Tracer et("updateQRZT");
233  int n = X.Ncols(); int s = X.Nrows();
234  if (s != L.Nrows())
235  Throw(ProgramException("Incompatible dimensions",X,L));
236  if (n == 0 || s == 0) return;
237  Real* xi = X.Store(); int k;
238  for (int i=0; i<s; i++)
239  {
240  Real r = L.element(i,i);
241  Real sum = 0.0;
242  Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
243  sum = sqrt(sum + square(r));
244  if (sum == 0.0)
245  {
246  REPORT
247  k=n; while(k--) { *xi0++ = 0.0; }
248  for (int j=i; j<s; j++) L.element(j,i) = 0.0;
249  }
250  else
251  {
252  Real frs = fabs(r) + sum;
253  Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
254  if (r <= 0) { REPORT L.element(i,i) = sum; alpha = -alpha; }
255  else { REPORT L.element(i,i) = -sum; }
256  Real* xj0=xi0; k=n; while(k--) { *xj0++ *= alpha; }
257  for (int j=i+1; j<s; j++)
258  {
259  sum = 0.0;
260  xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
261  sum += a0 * L.element(j,i);
262  xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
263  L.element(j,i) -= sum * a0;
264  }
265  }
266  }
267 }
268 
270 {
271  REPORT
272  Tracer et("updateQRZ");
273  int n = X.Nrows(); int s = X.Ncols();
274  if (s != U.Ncols())
275  Throw(ProgramException("Incompatible dimensions",X,U));
276  if (n == 0 || s == 0) return;
277  Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
278  RowVector V(s); Real* v0 = V.Store(); Real* v; V = 0.0;
279  int j, k; int J = s; int i = s;
280  while (i--)
281  {
282  Real* xj0 = xi0; Real* xi = xi0; k = n;
283  if (k) for (;;)
284  {
285  v = v0; Real Xi = *xi; Real* xj = xj0;
286  j = J; while(j--) *v++ += Xi * *xj++;
287  if (!(--k)) break;
288  xi += s; xj0 += s;
289  }
290 
291  Real r = *u0;
292  Real sum = sqrt(*v0 + square(r));
293 
294  if (sum == 0.0)
295  {
296  REPORT
297  u = u0; v = v0;
298  j = J; while(j--) { *u++ = 0.0; *v++ = 0.0; }
299  xj0 = xi0++; k = n;
300  if (k) for (;;)
301  {
302  *xj0 = 0.0;
303  if (!(--k)) break;
304  xj0 += s;
305  }
306  u0 += J--;
307  }
308  else
309  {
310  Real frs = fabs(r) + sum;
311  Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
312  if (r <= 0) { REPORT alpha = -alpha; *u0 = sum; }
313  else { REPORT *u0 = -sum; }
314 
315  j = J - 1; v = v0 + 1; u = u0 + 1;
316  while (j--)
317  { *v = a0 * *u + alpha * *v; *u -= a0 * *v; ++v; ++u; }
318 
319  xj0 = xi0; xi = xi0++; k = n;
320  if (k) for (;;)
321  {
322  v = v0 + 1; Real Xi = *xi; Real* xj = xj0;
323  Xi *= alpha; *xj++ = Xi;
324  j = J - 1; while(j--) *xj++ -= *v++ * Xi;
325  if (!(--k)) break;
326  xi += s; xj0 += s;
327  }
328 
329  j = J; v = v0;
330  while (j--) *v++ = 0.0;
331 
332  u0 += J--;
333  }
334  }
335 }
336 
337 // Matrix A's first n columns are orthonormal
338 // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
339 // Fill out the remaining columns of A to make them orthonormal
340 // so A.t() * A is the identity matrix
341 void extend_orthonormal(Matrix& A, int n)
342 {
343  REPORT
344  Tracer et("extend_orthonormal");
345  int nr = A.nrows(); int nc = A.ncols();
346  if (nc > nr) Throw(IncompatibleDimensionsException(A));
347  if (n > nc) Throw(IncompatibleDimensionsException(A));
348  ColumnVector SSR;
349  { Matrix A1 = A.Columns(1,n); SSR = A1.sum_square_rows(); }
350  for (int i = n; i < nc; ++i)
351  {
352  // pick row with smallest SSQ
353  int k; SSR.minimum1(k);
354  // orthogonalise column with 1 at element k, 0 elsewhere
355  // next line is rather inefficient
356  ColumnVector X = - A.Columns(1, i) * A.SubMatrix(k, k, 1, i).t();
357  X(k) += 1.0;
358  // normalise
359  X /= sqrt(X.SumSquare());
360  // update row sums of squares
361  for (k = 1; k <= nr; ++k) SSR(k) += square(X(k));
362  // load new column into matrix
363  A.Column(i+1) = X;
364  }
365 }
366 
367 
368 
369 
370 
371 #ifdef use_namespace
372 }
373 #endif
374 
375 
void QRZT(Matrix &X, LowerTriangularMatrix &L)
Definition: hholder.cpp:33
void updateQRZT(Matrix &X, LowerTriangularMatrix &L)
Definition: hholder.cpp:229
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
GetSubMatrix Columns(int f, int l) const
Definition: newmat.h:2153
Real & element(int, int)
Definition: newmat6.cpp:732
double Real
Definition: include.h:307
Real & element(int, int)
Definition: newmat6.cpp:682
int Nrows() const
Definition: newmat.h:494
Upper triangular matrix.
Definition: newmat.h:799
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
Real square(Real x)
Definition: hholder.cpp:31
TransposedMatrix t() const
Definition: newmat6.cpp:320
Real * Store() const
Definition: newmat.h:497
void extend_orthonormal(Matrix &A, int n)
Definition: hholder.cpp:341
#define REPORT
Definition: hholder.cpp:25
#define Throw(E)
Definition: myexcept.h:191
ReturnMatrix sum_square_rows() const
Definition: newmat8.cpp:736
The usual rectangular matrix.
Definition: newmat.h:625
Real SumSquare() const
Definition: newmat.h:346
Incompatible dimensions exception.
Definition: newmat.h:1998
int Ncols() const
Definition: newmat.h:495
void QRZ(Matrix &X, UpperTriangularMatrix &U)
Definition: hholder.cpp:117
void updateQRZ(Matrix &X, UpperTriangularMatrix &U)
Definition: hholder.cpp:269
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const
Definition: newmat.h:2146
int ncols() const
Definition: newmat.h:500
virtual void resize(int, int)
Definition: newmat4.cpp:289
FloatVector FloatVector FloatVector * alpha
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
int nrows() const
Definition: newmat.h:499
Row vector.
Definition: newmat.h:953
Column vector.
Definition: newmat.h:1008
Real minimum1(int &i) const
Definition: newmat8.cpp:286


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