source: ntrip/trunk/BNC/newmat/cholesky.cpp@ 6381

Last change on this file since 6381 was 2013, checked in by mervart, 15 years ago

* empty log message *

File size: 8.0 KB
RevLine 
[2013]1/// \ingroup newmat
2///@{
3
4/// \file cholesky.cpp
5/// Cholesky decomposition.
6/// Cholesky decomposition of symmetric and band symmetric matrices,
7/// update, downdate, manipulate a Cholesky decomposition
8
9
10// Copyright (C) 1991,2,3,4: R B Davies
11
12#define WANT_MATH
13//#define WANT_STREAM
14
15#include "include.h"
16
17#include "newmat.h"
18#include "newmatrm.h"
19
20#ifdef use_namespace
21namespace NEWMAT {
22#endif
23
24#ifdef DO_REPORT
25#define REPORT { static ExeCounter ExeCount(__LINE__,14); ++ExeCount; }
26#else
27#define REPORT {}
28#endif
29
30/********* Cholesky decomposition of a positive definite matrix *************/
31
32// Suppose S is symmetrix and positive definite. Then there exists a unique
33// lower triangular matrix L such that L L.t() = S;
34
35
36ReturnMatrix Cholesky(const SymmetricMatrix& S)
37{
38 REPORT
39 Tracer trace("Cholesky");
40 int nr = S.Nrows();
41 LowerTriangularMatrix T(nr);
42 Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
43 for (int i=0; i<nr; i++)
44 {
45 Real* tj = t; Real sum; int k;
46 for (int j=0; j<i; j++)
47 {
48 Real* tk = ti; sum = 0.0; k = j;
49 while (k--) { sum += *tj++ * *tk++; }
50 *tk = (*s++ - sum) / *tj++;
51 }
52 sum = 0.0; k = i;
53 while (k--) { sum += square(*ti++); }
54 Real d = *s++ - sum;
55 if (d<=0.0) Throw(NPDException(S));
56 *ti++ = sqrt(d);
57 }
58 T.release(); return T.for_return();
59}
60
61ReturnMatrix Cholesky(const SymmetricBandMatrix& S)
62{
63 REPORT
64 Tracer trace("Band-Cholesky");
65 int nr = S.Nrows(); int m = S.lower_val;
66 LowerBandMatrix T(nr,m);
67 Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
68
69 for (int i=0; i<nr; i++)
70 {
71 Real* tj = t; Real sum; int l;
72 if (i<m) { REPORT l = m-i; s += l; ti += l; l = i; }
73 else { REPORT t += (m+1); l = m; }
74
75 for (int j=0; j<l; j++)
76 {
77 Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
78 while (k--) { sum += *tj++ * *tk++; }
79 *tk = (*s++ - sum) / *tj++;
80 }
81 sum = 0.0;
82 while (l--) { sum += square(*ti++); }
83 Real d = *s++ - sum;
84 if (d<=0.0) Throw(NPDException(S));
85 *ti++ = sqrt(d);
86 }
87
88 T.release(); return T.for_return();
89}
90
91
92
93
94// Contributed by Nick Bennett of Schlumberger-Doll Research; modified by RBD
95
96// The enclosed routines can be used to update the Cholesky decomposition of
97// a positive definite symmetric matrix. A good reference for this routines
98// can be found in
99// LINPACK User's Guide, Chapter 10, Dongarra et. al., SIAM, Philadelphia, 1979
100
101// produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
102void update_Cholesky(UpperTriangularMatrix &chol, RowVector x)
103{
104 int nc = chol.Nrows();
105 ColumnVector cGivens(nc); cGivens = 0.0;
106 ColumnVector sGivens(nc); sGivens = 0.0;
107
108 for(int j = 1; j <= nc; ++j) // process the jth column of chol
109 {
110 // apply the previous Givens rotations k = 1,...,j-1 to column j
111 for(int k = 1; k < j; ++k)
112 GivensRotation(cGivens(k), sGivens(k), chol(k,j), x(j));
113
114 // determine the jth Given's rotation
115 pythag(chol(j,j), x(j), cGivens(j), sGivens(j));
116
117 // apply the jth Given's rotation
118 {
119 Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * x(j);
120 chol(j,j) = tmp0; x(j) = 0.0;
121 }
122
123 }
124
125}
126
127
128// produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
129void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x)
130{
131 int nRC = chol.Nrows();
132
133 // solve R^T a = x
134 LowerTriangularMatrix L = chol.t();
135 ColumnVector a(nRC); a = 0.0;
136 int i, j;
137
138 for (i = 1; i <= nRC; ++i)
139 {
140 // accumulate subtr sum
141 Real subtrsum = 0.0;
142 for(int k = 1; k < i; ++k) subtrsum += a(k) * L(i,k);
143
144 a(i) = (x(i) - subtrsum) / L(i,i);
145 }
146
147 // test that l2 norm of a is < 1
148 Real squareNormA = a.SumSquare();
149 if (squareNormA >= 1.0)
150 Throw(ProgramException("downdate_Cholesky() fails", chol));
151
152 Real alpha = sqrt(1.0 - squareNormA);
153
154 // compute and apply Givens rotations to the vector a
155 ColumnVector cGivens(nRC); cGivens = 0.0;
156 ColumnVector sGivens(nRC); sGivens = 0.0;
157 for(i = nRC; i >= 1; i--)
158 alpha = pythag(alpha, a(i), cGivens(i), sGivens(i));
159
160 // apply Givens rotations to the jth column of chol
161 ColumnVector xtilde(nRC); xtilde = 0.0;
162 for(j = nRC; j >= 1; j--)
163 {
164 // only the first j rotations have an affect on chol,0
165 for(int k = j; k >= 1; k--)
166 GivensRotation(cGivens(k), -sGivens(k), chol(k,j), xtilde(j));
167 }
168}
169
170
171
172// produces the Cholesky decomposition of EAE where A = chol.t() * chol
173// and E produces a RIGHT circular shift of the rows and columns from
174// 1,...,k-1,k,k+1,...l,l+1,...,p to
175// 1,...,k-1,l,k,k+1,...l-1,l+1,...p
176void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
177{
178 int nRC = chol.Nrows();
179 int i, j;
180
181 // I. compute shift of column l to the kth position
182 Matrix cholCopy = chol;
183 // a. grab column l
184 ColumnVector columnL = cholCopy.Column(l);
185 // b. shift columns k,...l-1 to the RIGHT
186 for(j = l-1; j >= k; --j)
187 cholCopy.Column(j+1) = cholCopy.Column(j);
188 // c. copy the top k-1 elements of columnL into the kth column of cholCopy
189 cholCopy.Column(k) = 0.0;
190 for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);
191
192 // II. determine the l-k Given's rotations
193 int nGivens = l-k;
194 ColumnVector cGivens(nGivens); cGivens = 0.0;
195 ColumnVector sGivens(nGivens); sGivens = 0.0;
196 for(i = l; i > k; i--)
197 {
198 int givensIndex = l-i+1;
199 columnL(i-1) = pythag(columnL(i-1), columnL(i),
200 cGivens(givensIndex), sGivens(givensIndex));
201 columnL(i) = 0.0;
202 }
203 // the kth entry of columnL is the new diagonal element in column k of cholCopy
204 cholCopy(k,k) = columnL(k);
205
206 // III. apply these Given's rotations to subsequent columns
207 // for columns k+1,...,l-1 we only need to apply the last nGivens-(j-k) rotations
208 for(j = k+1; j <= nRC; ++j)
209 {
210 ColumnVector columnJ = cholCopy.Column(j);
211 int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
212 for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
213 {
214 // apply gIndex Given's rotation
215 int topRowIndex = k + nGivens - gIndex;
216 GivensRotationR(cGivens(gIndex), sGivens(gIndex),
217 columnJ(topRowIndex), columnJ(topRowIndex+1));
218 }
219 cholCopy.Column(j) = columnJ;
220 }
221
222 chol << cholCopy;
223}
224
225
226
227// produces the Cholesky decomposition of EAE where A = chol.t() * chol
228// and E produces a LEFT circular shift of the rows and columns from
229// 1,...,k-1,k,k+1,...l,l+1,...,p to
230// 1,...,k-1,k+1,...l,k,l+1,...,p to
231void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
232{
233 int nRC = chol.Nrows();
234 int i, j;
235
236 // I. compute shift of column k to the lth position
237 Matrix cholCopy = chol;
238 // a. grab column k
239 ColumnVector columnK = cholCopy.Column(k);
240 // b. shift columns k+1,...l to the LEFT
241 for(j = k+1; j <= l; ++j)
242 cholCopy.Column(j-1) = cholCopy.Column(j);
243 // c. copy the elements of columnK into the lth column of cholCopy
244 cholCopy.Column(l) = 0.0;
245 for(i = 1; i <= k; ++i)
246 cholCopy(i,l) = columnK(i);
247
248 // II. apply and compute Given's rotations
249 int nGivens = l-k;
250 ColumnVector cGivens(nGivens); cGivens = 0.0;
251 ColumnVector sGivens(nGivens); sGivens = 0.0;
252 for(j = k; j <= nRC; ++j)
253 {
254 ColumnVector columnJ = cholCopy.Column(j);
255
256 // apply the previous Givens rotations to columnJ
257 int imax = j - k; if (imax > nGivens) imax = nGivens;
258 for(int i = 1; i <= imax; ++i)
259 {
260 int gIndex = i;
261 int topRowIndex = k + i - 1;
262 GivensRotationR(cGivens(gIndex), sGivens(gIndex),
263 columnJ(topRowIndex), columnJ(topRowIndex+1));
264 }
265
266 // compute a new Given's rotation when j < l
267 if(j < l)
268 {
269 int gIndex = j-k+1;
270 columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
271 sGivens(gIndex));
272 columnJ(j+1) = 0.0;
273 }
274
275 cholCopy.Column(j) = columnJ;
276 }
277
278 chol << cholCopy;
279
280}
281
282
283
284
285#ifdef use_namespace
286}
287#endif
288
289///@}
Note: See TracBrowser for help on using the repository browser.