25 #define REPORT { static ExeCounter ExeCount(__LINE__,14); ++ExeCount; } 43 for (
int i=0; i<nr; i++)
46 for (
int j=0; j<i; j++)
48 Real* tk = ti; sum = 0.0; k = j;
49 while (k--) { sum += *tj++ * *tk++; }
50 *tk = (*s++ -
sum) / *tj++;
53 while (k--) { sum +=
square(*ti++); }
69 for (
int i=0; i<nr; i++)
72 if (i<m) {
REPORT l = m-i; s += l; ti += l; l = i; }
73 else {
REPORT t += (m+1); l =
m; }
75 for (
int j=0; j<l; j++)
77 Real* tk = ti; sum = 0.0;
int k = j; tj += (m-j);
78 while (k--) { sum += *tj++ * *tk++; }
79 *tk = (*s++ -
sum) / *tj++;
82 while (l--) { sum +=
square(*ti++); }
104 int nc = chol.
Nrows();
108 for(
int j = 1; j <= nc; ++j)
111 for(
int k = 1; k < j; ++k)
115 pythag(chol(j,j), x(j), cGivens(j), sGivens(j));
119 Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * x(j);
120 chol(j,j) = tmp0; x(j) = 0.0;
131 int nRC = chol.
Nrows();
138 for (i = 1; i <= nRC; ++i)
142 for(
int k = 1; k < i; ++k) subtrsum +=
a(k) * L(i,k);
144 a(i) = (x(i) - subtrsum) / L(i,i);
149 if (squareNormA >= 1.0)
157 for(i = nRC; i >= 1; i--)
158 alpha =
pythag(alpha,
a(i), cGivens(i), sGivens(i));
162 for(j = nRC; j >= 1; j--)
165 for(
int k = j; k >= 1; k--)
178 int nRC = chol.
Nrows();
186 for(j = l-1; j >= k; --j)
190 for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);
196 for(i = l; i > k; i--)
198 int givensIndex = l-i+1;
199 columnL(i-1) =
pythag(columnL(i-1), columnL(i),
200 cGivens(givensIndex), sGivens(givensIndex));
204 cholCopy(k,k) = columnL(k);
208 for(j = k+1; j <= nRC; ++j)
211 int imin = nGivens - (j-k) + 1;
if (imin < 1) imin = 1;
212 for(
int gIndex = imin; gIndex <= nGivens; ++gIndex)
215 int topRowIndex = k + nGivens - gIndex;
217 columnJ(topRowIndex), columnJ(topRowIndex+1));
219 cholCopy.
Column(j) = columnJ;
233 int nRC = chol.
Nrows();
241 for(j = k+1; j <= l; ++j)
245 for(i = 1; i <= k; ++i)
246 cholCopy(i,l) = columnK(i);
252 for(j = k; j <= nRC; ++j)
257 int imax = j - k;
if (imax > nGivens) imax = nGivens;
258 for(
int i = 1; i <= imax; ++i)
261 int topRowIndex = k + i - 1;
263 columnJ(topRowIndex), columnJ(topRowIndex+1));
270 columnJ(j) =
pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
275 cholCopy.
Column(j) = columnJ;
void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x)
void GivensRotation(Real cGivens, Real sGivens, Real &x, Real &y)
void GivensRotationR(Real cGivens, Real sGivens, Real &x, Real &y)
Miscellaneous exception (details in character string).
ReturnMatrix for_return() const
Real pythag(Real f, Real g, Real &c, Real &s)
void update_Cholesky(UpperTriangularMatrix &chol, RowVector x)
GetSubMatrix Column(int f) const
TransposedMatrix t() const
The usual rectangular matrix.
void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
ReturnMatrix Cholesky(const SymmetricMatrix &S)
FloatVector FloatVector * a
Real trace(const BaseMatrix &B)
void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
Lower triangular band matrix.
FloatVector FloatVector FloatVector * alpha
Real sum(const BaseMatrix &B)
Not positive definite exception.