00001
00002
00003
00008
00009
00010
00011
00012 #define WANT_MATH
00013
00014
00015 #include "include.h"
00016
00017 #include "newmat.h"
00018 #include "newmatrm.h"
00019
00020 #ifdef use_namespace
00021 namespace NEWMAT {
00022 #endif
00023
00024 #ifdef DO_REPORT
00025 #define REPORT { static ExeCounter ExeCount(__LINE__,14); ++ExeCount; }
00026 #else
00027 #define REPORT {}
00028 #endif
00029
00030
00031
00032
00033
00034
00035
00036 ReturnMatrix Cholesky(const SymmetricMatrix& S)
00037 {
00038 REPORT
00039 Tracer trace("Cholesky");
00040 int nr = S.Nrows();
00041 LowerTriangularMatrix T(nr);
00042 Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
00043 for (int i=0; i<nr; i++)
00044 {
00045 Real* tj = t; Real sum; int k;
00046 for (int j=0; j<i; j++)
00047 {
00048 Real* tk = ti; sum = 0.0; k = j;
00049 while (k--) { sum += *tj++ * *tk++; }
00050 *tk = (*s++ - sum) / *tj++;
00051 }
00052 sum = 0.0; k = i;
00053 while (k--) { sum += square(*ti++); }
00054 Real d = *s++ - sum;
00055 if (d<=0.0) Throw(NPDException(S));
00056 *ti++ = sqrt(d);
00057 }
00058 T.release(); return T.for_return();
00059 }
00060
00061 ReturnMatrix Cholesky(const SymmetricBandMatrix& S)
00062 {
00063 REPORT
00064 Tracer trace("Band-Cholesky");
00065 int nr = S.Nrows(); int m = S.lower_val;
00066 LowerBandMatrix T(nr,m);
00067 Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
00068
00069 for (int i=0; i<nr; i++)
00070 {
00071 Real* tj = t; Real sum; int l;
00072 if (i<m) { REPORT l = m-i; s += l; ti += l; l = i; }
00073 else { REPORT t += (m+1); l = m; }
00074
00075 for (int j=0; j<l; j++)
00076 {
00077 Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
00078 while (k--) { sum += *tj++ * *tk++; }
00079 *tk = (*s++ - sum) / *tj++;
00080 }
00081 sum = 0.0;
00082 while (l--) { sum += square(*ti++); }
00083 Real d = *s++ - sum;
00084 if (d<=0.0) Throw(NPDException(S));
00085 *ti++ = sqrt(d);
00086 }
00087
00088 T.release(); return T.for_return();
00089 }
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 void update_Cholesky(UpperTriangularMatrix &chol, RowVector x)
00103 {
00104 int nc = chol.Nrows();
00105 ColumnVector cGivens(nc); cGivens = 0.0;
00106 ColumnVector sGivens(nc); sGivens = 0.0;
00107
00108 for(int j = 1; j <= nc; ++j)
00109 {
00110
00111 for(int k = 1; k < j; ++k)
00112 GivensRotation(cGivens(k), sGivens(k), chol(k,j), x(j));
00113
00114
00115 pythag(chol(j,j), x(j), cGivens(j), sGivens(j));
00116
00117
00118 {
00119 Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * x(j);
00120 chol(j,j) = tmp0; x(j) = 0.0;
00121 }
00122
00123 }
00124
00125 }
00126
00127
00128
00129 void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x)
00130 {
00131 int nRC = chol.Nrows();
00132
00133
00134 LowerTriangularMatrix L = chol.t();
00135 ColumnVector a(nRC); a = 0.0;
00136 int i, j;
00137
00138 for (i = 1; i <= nRC; ++i)
00139 {
00140
00141 Real subtrsum = 0.0;
00142 for(int k = 1; k < i; ++k) subtrsum += a(k) * L(i,k);
00143
00144 a(i) = (x(i) - subtrsum) / L(i,i);
00145 }
00146
00147
00148 Real squareNormA = a.SumSquare();
00149 if (squareNormA >= 1.0)
00150 Throw(ProgramException("downdate_Cholesky() fails", chol));
00151
00152 Real alpha = sqrt(1.0 - squareNormA);
00153
00154
00155 ColumnVector cGivens(nRC); cGivens = 0.0;
00156 ColumnVector sGivens(nRC); sGivens = 0.0;
00157 for(i = nRC; i >= 1; i--)
00158 alpha = pythag(alpha, a(i), cGivens(i), sGivens(i));
00159
00160
00161 ColumnVector xtilde(nRC); xtilde = 0.0;
00162 for(j = nRC; j >= 1; j--)
00163 {
00164
00165 for(int k = j; k >= 1; k--)
00166 GivensRotation(cGivens(k), -sGivens(k), chol(k,j), xtilde(j));
00167 }
00168 }
00169
00170
00171
00172
00173
00174
00175
00176 void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
00177 {
00178 int nRC = chol.Nrows();
00179 int i, j;
00180
00181
00182 Matrix cholCopy = chol;
00183
00184 ColumnVector columnL = cholCopy.Column(l);
00185
00186 for(j = l-1; j >= k; --j)
00187 cholCopy.Column(j+1) = cholCopy.Column(j);
00188
00189 cholCopy.Column(k) = 0.0;
00190 for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);
00191
00192
00193 int nGivens = l-k;
00194 ColumnVector cGivens(nGivens); cGivens = 0.0;
00195 ColumnVector sGivens(nGivens); sGivens = 0.0;
00196 for(i = l; i > k; i--)
00197 {
00198 int givensIndex = l-i+1;
00199 columnL(i-1) = pythag(columnL(i-1), columnL(i),
00200 cGivens(givensIndex), sGivens(givensIndex));
00201 columnL(i) = 0.0;
00202 }
00203
00204 cholCopy(k,k) = columnL(k);
00205
00206
00207
00208 for(j = k+1; j <= nRC; ++j)
00209 {
00210 ColumnVector columnJ = cholCopy.Column(j);
00211 int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
00212 for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
00213 {
00214
00215 int topRowIndex = k + nGivens - gIndex;
00216 GivensRotationR(cGivens(gIndex), sGivens(gIndex),
00217 columnJ(topRowIndex), columnJ(topRowIndex+1));
00218 }
00219 cholCopy.Column(j) = columnJ;
00220 }
00221
00222 chol << cholCopy;
00223 }
00224
00225
00226
00227
00228
00229
00230
00231 void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
00232 {
00233 int nRC = chol.Nrows();
00234 int i, j;
00235
00236
00237 Matrix cholCopy = chol;
00238
00239 ColumnVector columnK = cholCopy.Column(k);
00240
00241 for(j = k+1; j <= l; ++j)
00242 cholCopy.Column(j-1) = cholCopy.Column(j);
00243
00244 cholCopy.Column(l) = 0.0;
00245 for(i = 1; i <= k; ++i)
00246 cholCopy(i,l) = columnK(i);
00247
00248
00249 int nGivens = l-k;
00250 ColumnVector cGivens(nGivens); cGivens = 0.0;
00251 ColumnVector sGivens(nGivens); sGivens = 0.0;
00252 for(j = k; j <= nRC; ++j)
00253 {
00254 ColumnVector columnJ = cholCopy.Column(j);
00255
00256
00257 int imax = j - k; if (imax > nGivens) imax = nGivens;
00258 for(int i = 1; i <= imax; ++i)
00259 {
00260 int gIndex = i;
00261 int topRowIndex = k + i - 1;
00262 GivensRotationR(cGivens(gIndex), sGivens(gIndex),
00263 columnJ(topRowIndex), columnJ(topRowIndex+1));
00264 }
00265
00266
00267 if(j < l)
00268 {
00269 int gIndex = j-k+1;
00270 columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
00271 sGivens(gIndex));
00272 columnJ(j+1) = 0.0;
00273 }
00274
00275 cholCopy.Column(j) = columnJ;
00276 }
00277
00278 chol << cholCopy;
00279
00280 }
00281
00282
00283
00284
00285 #ifdef use_namespace
00286 }
00287 #endif
00288