16 using namespace NEWMAT;
73 Tracer et1(
"LU1 - BandLU");
80 Tracer et1(
"LU2 - BandLU");
87 Tracer et1(
"LU3 - BandLU");
125 : m(mx), n1(n1x), n2(n2x), n3(n3x) { Reset(); }
135 Tracer et(
"Thirteenth test of Matrix package");
139 for (j=1;j<=20;j++) X(1,j) = j+1;
140 for (i=2;i<=5;i++)
for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
152 Diff = L - L1;
Clean(Diff,0.000000001);
Print(Diff);
155 Diff = L - Ut.
t();
Clean(Diff,0.000000001);
Print(Diff);
157 for (j=1;j<=20;j++) Y(1,j) = 22-j;
158 for (i=2;i<=3;i++)
for (j=1;j<=20; j++)
159 Y(i,j) = (long)Y(i-1,j) * j % 101;
162 Diff = Xt - X.
t();
Clean(Diff,0.000000001);
Print(Diff);
163 Diff = Yt - Y.
t();
Clean(Diff,0.000000001);
Print(Diff);
164 Diff = Mt - M.
t();
Clean(Diff,0.000000001);
Print(Diff);
165 Diff = Y2 * Xt2 * S.
i() - M * L.
i();
173 for (j=1;j<=5;j++) X(1,j) = j+1;
174 for (i=2;i<=5;i++)
for (j=1;j<=5; j++)
175 X(i,j) = (long)X(i-1,j) * j % 1001;
176 for (i=1;i<=5;i++) C1(i) = i*i;
185 for (j=1;j<=7;j++) X(1,j) = j+1;
186 for (i=2;i<=7;i++)
for (j=1;j<=7; j++)
187 X(i,j) = (long)X(i-1,j) * j % 1001;
189 for (i=1;i<=7;i++) C1(i) = i*i;
191 Diff = R1 * X.
i() - ( X.
t().
i() * R1.
t() ).t();
Clean(Diff,0.000000001);
198 for (j=1;j<=5;j++) X(1,j) = j+1;
199 for (i=2;i<=5;i++)
for (j=1;j<=5; j++)
200 X(i,j) = (long)X(i-1,j) * j % 1001;
202 for (i=1;i<=5;i++) C1(i) = i*i;
216 if (i<=n-1) B(i,i+1) = -4;
217 if (i<=n-2) B(i,i+2) = 1;
219 B(1,1) = 5; B(n,n) = 5;
223 for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
234 Y2 = CU * (CL * X) - Y1;
236 Y2 = B.
i() * X - Y1;
Clean(Y2, 0.000000001);
Print(Y2);
244 Y3 = CU * (CL * X) - Y1;
247 Y3 = B.
i() * X - Y1;
Clean(Y3, 0.000000001);
Print(Y3);
250 Y2 = AI*X - Y1;
Clean(Y2, 0.000000001);
Print(Y2);
258 M = A; AI << M; M = AI-A;
Clean(M, 0.000000001);
Print(M);
259 C = B; BI << C; M = BI-B;
Clean(M, 0.000000001);
Print(M);
276 M = A *
Cholesky(B); M = M * M.
t() - A * B * A;
286 for (i=1;i<=N;i++) { S(i,i)=1; B(i,i)=1; B(i,i+1)=-1; }
287 for (i=1;i<N; i++) S(i,i+1)=-.5;
289 B = B.
t()*S.
i()*B - (D-1.0/(N+1))*2.0;
297 A.
Row(1) << 3 << 2 << -1 << 4 << -3 << 5 << 9;
298 A.
Row(2) << -8 << 7 << 2 << 0 << 7 << 0 << -1;
299 A.
Row(3) << 2 << -2 << 3 << 1 << 9 << 0 << 3;
300 A.
Row(4) << -1 << 5 << 2 << 2 << 5 << -1 << 2;
301 A.
Row(5) << 4 << -4 << 1 << 9 << -8 << 7 << 5;
302 A.
Row(6) << 1 << -2 << 5 << -1 << -2 << 5 << 1;
303 A.
Row(7) << -6 << 3 << -1 << 8 << -1 << 2 << 2;
313 C = A *
Inverter2(B1) - IdentityMatrix(7);
320 C = A *
Inverter2(B1) - IdentityMatrix(7);
326 C = -A *
Inverter2(B1) - IdentityMatrix(7);
330 C = A * B.
i() - IdentityMatrix(7);
Clean(C, 0.000000001);
Print(C);
331 D(10) = (B == B1 ? 0 : 1) + (B == B2 ? 0 : 1);
336 B2 = B1; D(15) = B == B2 ? 0 : 1;
340 B1 = B; B1 = B1.
i(); C = A - B1.
i();
Clean(C, 0.000000001);
Print(C);
341 B1 = B; B1.
release(); B1 = B1; B2 = B1;
342 D(19) = B == B1 ? 0 : 1; D(20) = B == B2 ? 0 : 1;
345 B1 = GM; D(23) = B == B1 ? 0 : 1;
346 B1 = A * 0; B2 = B1; D(24) = B2.
is_singular() ? 0 : 1;
356 D(25) = (d != d1 ? 0 : 1) + (d1 == d2 ? 0 : 1)
357 + (i != i1 ? 0 : 1) + (i1 == i2 ? 0 : 1);
366 A.
Row(1) << 3 << 2 << -1;
367 A.
Row(2) << -8 << 7 << 2 << 0;
368 A.
Row(3) << 2 << -2 << 3 << 1 << 9;
369 A.
Row(4) << -1 << 5 << 2 << 2 << 5 << -1;
370 A.
Row(5) << -4 << 1 << 9 << -8 << 7 << 5;
371 A.
Row(6) << 5 << -1 << -2 << 5 << 1;
372 A.
Row(7) << 8 << -1 << 2 << 2;
382 C = A *
Inverter2(B1) - IdentityMatrix(7);
389 C = A *
Inverter2(B1) - IdentityMatrix(7);
395 C = -A *
Inverter2(B1) - IdentityMatrix(7);
399 C = A * B.
i() - IdentityMatrix(7);
Clean(C, 0.000000001);
Print(C);
400 D(10) = (B == B1 ? 0 : 1) + (B == B2 ? 0 : 1);
405 D(14) = B.
size2()-21;
408 B2 = B1; D(15) = B == B2 ? 0 : 1;
416 B1 = B; CM = B1.
i(); C = A - CM.
i();
Clean(C, 0.000000001);
Print(C);
417 B1 = B; B1.
release(); B1 = B1; B2 = B1;
418 D(19) = B == B1 ? 0 : 1; D(20) = B == B2 ? 0 : 1;
421 B1 = GM; D(23) = B == B1 ? 0 : 1;
422 B1 = A * 0; B2 = B1; D(24) = B2.
is_singular() ? 0 : 1;
432 D(25) = (d != d1 ? 0 : 1) + (d1 == d2 ? 0 : 1)
433 + (dd != dd1 ? 0 : 1) + (dd1 == dd2 ? 0 : 1)
434 + (i != i1 ? 0 : 1) + (i1 == i2 ? 0 : 1);
448 for (i = 1; i <= 100; ++i)
for (j = 1; j <= 10; ++j)
449 X(i, j) = 2.0 * (mwc.
Next() - 0.5);
463 for (j = 1; j <= 10; ++j) NewRow(j) = 2.0 * (mwc.
Next() - 0.5);
465 X = X1 & NewRow;
QRZ(X, U1);
466 Diff = U1 - U2;
Clean(Diff, 0.000000001);
Print(Diff);
471 X = X1.
Rows(1,19) & X1.
Rows(21,34) & X1.
Rows(36,100) & NewRow;
QRZ(X, U1);
472 Diff = U1 - U2;
Clean(Diff, 0.000000001);
Print(Diff);
507 X1.ReSize(
m, n1); X2.ReSize(
m, n2); X3.ReSize(
m, n3);
508 for (i = 1; i <=
m; ++i)
for (j = 1; j <= n1; ++j)
509 X1(i, j) = 2.0 * (mwc.Next() - 0.5);
510 for (i = 1; i <=
m; ++i)
for (j = 1; j <= n2; ++j)
511 X2(i, j) = 2.0 * (mwc.Next() - 0.5);
512 for (i = 1; i <=
m; ++i)
for (j = 1; j <= n3; ++j)
513 X3(i, j) = 2.0 * (mwc.Next() - 0.5);
518 Matrix XT1 = X1.
t(), XT2 = X2.
t(), XT3 = X3.
t();
528 Diff = SM - LX * LX.
t();
Clean(Diff, 0.000000001);
Print(Diff);
529 Diff = SM - LY * LY.
t();
Clean(Diff, 0.000000001);
Print(Diff);
533 Diff = U.
t() * U - SM;
Clean(Diff, 0.000000001);
Print(Diff);
537 Diff = SM - UY.
t() * UY;
Clean(Diff, 0.000000001);
Print(Diff);
539 Diff = SM - UX.
t() * UX;
Clean(Diff, 0.000000001);
Print(Diff);
void QRZT(Matrix &X, LowerTriangularMatrix &L)
void UpdateCholesky(UpperTriangularMatrix &chol, const RowVector &x)
TestUpdateQRZ(int mx, int n1x, int n2x=0, int n3x=0)
void UpdateQRZT(Matrix &X, LowerTriangularMatrix &L)
ReturnMatrix LU2(const Matrix &A)
ReturnMatrix Inverter2(CroutMatrix X)
void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector &x)
ReturnMatrix for_return() const
GetSubMatrix Columns(int f, int l) const
ReturnMatrix Inverter1(const CroutMatrix &X)
void release_and_delete()
const int * const_data_indx() const
virtual void ReSize(int m, int n)
const int * const_data_indx() const
void LeftCircularUpdateCholesky(UpperTriangularMatrix &chol, int k, int l)
ReturnMatrix LU1(const Matrix &A)
MatrixBandWidth bandwidth() const
void SetRow(int i, int j)
void RightCircularUpdateCholesky(UpperTriangularMatrix &chol, int k, int l)
ReturnMatrix LU3(const Matrix &A)
void Clean(Matrix &A, Real c)
GetSubMatrix Column(int f) const
TransposedMatrix t() const
const Real * const_data() const
The usual rectangular matrix.
ReturnMatrix Cholesky(const SymmetricMatrix &S)
void QRZ(Matrix &X, UpperTriangularMatrix &U)
LU decomposition of a band matrix.
void MatrixErrorNoSpace(const void *)
test for allocation fails
GetSubMatrix Row(int f) const
Lower triangular band matrix.
ReturnMatrix ForReturn() const
void UpdateQRZ(Matrix &X, UpperTriangularMatrix &U)
Real determinant(const BaseMatrix &B)
void CircularShift(const Matrix &X1, int first, int last)
void Print(const Matrix &X)
A matrix which can be of any GeneralMatrix type.
GetSubMatrix Rows(int f, int l) const