18 using namespace NEWMAT;
21 static int my_max(
int p,
int q) {
return (p > q) ? p : q; }
23 static int my_min(
int p,
int q) {
return (p < q) ? p : q; }
30 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
31 if (i - j <= l1 && i - j >= -u1) BM1(i, j) = 100 * i + j;
35 BM2.
ReSize(20, l2, u2); BM2 = 777777.0;
36 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
37 if (i - j <= l2 && i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
39 BandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP =
SP(BM1, BM2),
40 BMM = BM1 * BM2, BMN = -BM1;
43 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
44 if (i - j <= l1 && i - j >= -u1) M1(i, j) = 100 * i + j;
47 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
48 if (i - j <= l2 && i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
50 Matrix MA = M1 + M2, MS = M1 - M2, MSP =
SP(M1, M2), MM = M1 * M2, MN = -M1;
51 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
73 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
74 if (i - j <= l1 && i - j >= -u1) BM1E.
element(i-1, j-1) = 100 * i + j;
76 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
77 if (i - j <= l2 && i - j >= -u2)
78 BM2E.
element(i-1, j-1) = (100 * i + j) * 11;
79 M1 = BM1E - BM1;
Print(M1);
80 M2 = BM2E - BM2;
Print(M2);
83 BM1E = 444444.04; BM2E = 555555.0;
85 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
86 if (i - j <= l1 && i - j >= -u1)
88 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
89 if (i - j <= l2 && i - j >= -u2)
90 BM2E.
element(i-1, j-1) = BM2C.element(i-1, j-1);
91 M1 = BM1E - BM1;
Print(M1);
92 M2 = BM2E - BM2;
Print(M2);
95 BM1E = 444444.04; BM2E = 555555.0;
96 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
97 if (i - j <= l1 && i - j >= -u1) BM1E(i, j) = BM1C(i, j);
98 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
99 if (i - j <= l2 && i - j >= -u2) BM2E(i, j) = BM2C(i, j);
100 M1 = BM1E - BM1;
Print(M1);
101 M2 = BM2E - BM2;
Print(M2);
108 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
109 if (i - j <= l1) BM1(i, j) = 100 * i + j;
113 BM2.
ReSize(20, l2); BM2 = 777777.0;
114 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
115 if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
118 BMM = BM1 * BM2, BMN = -BM1;
121 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
122 if (i - j <= l1) M1(i, j) = 100 * i + j;
125 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
126 if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
128 Matrix MA = M1 + M2, MS = M1 - M2, MSP =
SP(M1, M2), MM = M1 * M2, MN = -M1;
129 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
151 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
152 if (i - j <= l1) BM1E.
element(i-1, j-1) = 100 * i + j;
154 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
155 if (i - j <= l2) BM2E.
element(i-1, j-1) = (100 * i + j) * 11;
156 M1 = BM1E - BM1;
Print(M1);
157 M2 = BM2E - BM2;
Print(M2);
160 BM1E = 444444.04; BM2E = 555555.0;
162 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
164 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
165 if (i - j <= l2) BM2E.
element(i-1, j-1) = BM2C.element(i-1, j-1);
166 M1 = BM1E - BM1;
Print(M1);
167 M2 = BM2E - BM2;
Print(M2);
170 BM1E = 444444.04; BM2E = 555555.0;
171 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
172 if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
173 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
174 if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
175 M1 = BM1E - BM1;
Print(M1);
176 M2 = BM2E - BM2;
Print(M2);
183 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
184 if (i - j >= -u1) BM1(i, j) = 100 * i + j;
188 BM2.
ReSize(20, u2); BM2 = 777777.0;
189 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
190 if (i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
193 BMM = BM1 * BM2, BMN = -BM1;
196 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
197 if (i - j >= -u1) M1(i, j) = 100 * i + j;
200 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
201 if (i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
203 Matrix MA = M1 + M2, MS = M1 - M2, MSP =
SP(M1, M2), MM = M1 * M2, MN = -M1;
204 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
226 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
227 if (i - j >= -u1) BM1E.
element(i-1, j-1) = 100 * i + j;
229 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
230 if (i - j >= -u2) BM2E.
element(i-1, j-1) = (100 * i + j) * 11;
231 M1 = BM1E - BM1;
Print(M1);
232 M2 = BM2E - BM2;
Print(M2);
235 BM1E = 444444.04; BM2E = 555555.0;
237 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
238 if (i - j >= -u1) BM1E.
element(i-1, j-1) = BM1C.
element(i-1, j-1);
239 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
240 if (i - j >= -u2) BM2E.
element(i-1, j-1) = BM2C.element(i-1, j-1);
241 M1 = BM1E - BM1;
Print(M1);
242 M2 = BM2E - BM2;
Print(M2);
245 BM1E = 444444.04; BM2E = 555555.0;
246 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
247 if (i - j >= -u1) BM1E(i, j) = BM1C(i, j);
248 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
249 if (i - j >= -u2) BM2E(i, j) = BM2C(i, j);
250 M1 = BM1E - BM1;
Print(M1);
251 M2 = BM2E - BM2;
Print(M2);
258 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
259 if (i - j <= l1) BM1(i, j) = 100 * i + j;
263 BM2.
ReSize(20, l2); BM2 = 777777.0;
264 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
265 if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
272 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
273 if (i - j <= l1) M1(i, j) = 100 * i + j;
276 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
277 if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
281 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
303 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
304 if (i - j <= l1) BM1E.
element(i-1, j-1) = 100 * i + j;
306 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
307 if (i - j <= l2) BM2E.
element(i-1, j-1) = (100 * i + j) * 11;
308 M1 = BM1E - BM1;
Print(M1);
309 M2 = BM2E - BM2;
Print(M2);
313 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
314 if (i - j <= l1) BM1E.
element(j-1, i-1) = 100 * i + j;
316 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
317 if (i - j <= l2) BM2E.
element(j-1, i-1) = (100 * i + j) * 11;
318 M1 = BM1E - BM1;
Print(M1);
319 M2 = BM2E - BM2;
Print(M2);
322 BM1E = 444444.04; BM2E = 555555.0;
324 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
326 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
327 if (i - j <= l2) BM2E.
element(i-1, j-1) = BM2C.element(i-1, j-1);
328 M1 = BM1E - BM1;
Print(M1);
329 M2 = BM2E - BM2;
Print(M2);
332 BM1E = 444444.0; BM2E = 555555.0;
333 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
335 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
336 if (i - j <= l2) BM2E.
element(j-1, i-1) = BM2C.element(j-1, i-1);
337 M1 = BM1E - BM1;
Print(M1);
338 M2 = BM2E - BM2;
Print(M2);
341 BM1E = 444444.0; BM2E = 555555.0;
342 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
343 if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
344 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
345 if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
346 M1 = BM1E - BM1;
Print(M1);
347 M2 = BM2E - BM2;
Print(M2);
350 BM1E = 444444.0; BM2E = 555555.0;
351 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
352 if (i - j <= l1) BM1E(j, i) = BM1C(j, i);
353 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
354 if (i - j <= l2) BM2E(j, i) = BM2C(j, i);
355 M1 = BM1E - BM1;
Print(M1);
356 M2 = BM2E - BM2;
Print(M2);
359 BM1E = 444444.0; BM2E = 555555.0;
360 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
361 if (i - j <= l1) BM1E(j, i) = BM1C(i, j);
362 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
363 if (i - j <= l2) BM2E(j, i) = BM2C(i, j);
364 M1 = BM1E - BM1;
Print(M1);
365 M2 = BM2E - BM2;
Print(M2);
375 Tracer et(
"Seventeenth test of Matrix package");
383 for (i=1; i<=8; i++)
for (j=-3; j<=1; j++)
384 {
int k = i+j;
if (k>0 && k<=8) B(i,k) = i + k/64.0; }
388 Print(
Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
390 Print(
Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
398 BI = BLU.
i(); A = X.
i()-BI;
Clean(A,0.000000001);
Print(A);
416 for (i=1; i<=7; i++)
for (j=1; j<=7; j++)
418 int k=i-j;
if (k<0) k = -k;
420 else if (k==1) A(i,j) = -4;
421 else if (k==2) A(i,j) = 1;
433 V = V / 64 - 1;
Clean(V,0.000000001);
Print(V);
436 Real a[] = {1,2,3,4,5,6,7};
438 M = (M.
i()*X).t() - (B.
i()*X).t() * 2.0 + (A.
i()*X).t();
443 for (i=1; i<=80; i++)
for (j=1; j<=80; j++)
444 {
int d = i-j;
if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
445 for (i=1; i<=80; i++) CX(i) = i*100.0;
450 V1 = P * V1; V2 = MP * V2; V = V1 - V2;
Clean(V,0.000000001);
Print(V);
456 for (i=1; i<=80; i++)
for (j=1; j<=80; j++)
457 {
int d = i-j;
if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
462 V1 = LP.
i() * CX; V2 = MP.
i() * CX;
466 for (i=1; i<=80; i++)
for (j=1; j<=80; j++)
467 {
int d = i-j;
if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
471 V1 = UP.
i() * CX; V2 = MP.
i() * CX;
480 for (i=1; i<=8; i++)
for (j=1; j<=8; j++)
481 if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
489 D = 10; SA = SA.
t() + D; MA1 = MA1 + D;
493 D << SA; UB << SA; LB << SA;
494 SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
498 for (i=1; i<=7; i++)
for (j=1; j<=7; j++)
502 else if (k==1) A(i,j) = -4;
503 else if (k==2) A(i,j) = 1;
518 Real a[] = {1,2,3,4,5,6,7};
521 X = M.
i()*X - C.
i()*X * 2 + A.
i()*X;
525 LB << A; UB << A; D << A;
529 for (i=1; i<=7; i++)
for (j=1; j<=7; j++)
533 else if (k==1) A(i,j) = -4;
534 else if (k==2) A(i,j) = 1;
539 M = LB.
i() * LB - D;
Clean(M,0.000000001);
Print(M);
540 M = UB.
i() * UB - D;
Clean(M,0.000000001);
Print(M);
541 M = XA.
i() * XA - D;
Clean(M,0.000000001);
Print(M);
543 M = LB.
i() * UB - LB.
i() * MUB;
Clean(M,0.000000001);
Print(M);
544 M = UB.
i() * LB - UB.
i() * MLB;
Clean(M,0.000000001);
Print(M);
590 BD = I; UBD = I; LBD = I; SBD = I;
610 Matrix Count(1, 1); Count = 0;
612 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 30; ++j)
613 M(i, j) = 100 * i + j;
615 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 30; ++j)
617 if (M(i, j) != CM(i, j)) ++Count(1,1);
618 if (M(i, j) != CM.
element(i-1, j-1)) ++Count(1,1);
619 if (M(i, j) != M.
element(i-1, j-1)) ++Count(1,1);
623 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
624 U(i, j) = 100 * i + j;
626 for (i = 1; i <= 20; ++i)
for (j = i; j <= 20; ++j)
628 if (U(i, j) != CU(i, j)) ++Count(1,1);
629 if (U(i, j) != CU.
element(i-1, j-1)) ++Count(1,1);
630 if (U(i, j) != U.
element(i-1, j-1)) ++Count(1,1);
634 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
635 L(i, j) = 100 * i + j;
637 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
639 if (L(i, j) != CL(i, j)) ++Count(1,1);
640 if (L(i, j) != CL.
element(i-1, j-1)) ++Count(1,1);
641 if (L(i, j) != L.
element(i-1, j-1)) ++Count(1,1);
645 for (i = 1; i <= 20; ++i)
for (j = 1; j <= i; ++j)
646 S(i, j) = 100 * i + j;
648 for (i = 1; i <= 20; ++i)
for (j = 1; j <= 20; ++j)
650 if (S(i, j) != CS(i, j)) ++Count(1,1);
651 if (S(i, j) != CS.
element(i-1, j-1)) ++Count(1,1);
652 if (S(i, j) != S.
element(i-1, j-1)) ++Count(1,1);
653 if (S(i, j) != S(j, i)) ++Count(1,1);
654 if (S(i, j) != CS(i, j)) ++Count(1,1);
655 if (S(i, j) != CS.
element(i-1, j-1)) ++Count(1,1);
656 if (S(i, j) != S.
element(i-1, j-1)) ++Count(1,1);
660 for (i = 1; i <= 20; ++i) D(i) = 100 * i + i * i;
662 for (i = 1; i <= 20; ++i)
664 if (D(i, i) != CD(i, i)) ++Count(1,1);
665 if (D(i, i) != CD.
element(i-1, i-1)) ++Count(1,1);
666 if (D(i, i) != D.
element(i-1, i-1)) ++Count(1,1);
667 if (D(i, i) != D(i)) ++Count(1,1);
668 if (D(i) != CD(i)) ++Count(1,1);
669 if (D(i) != CD.
element(i-1)) ++Count(1,1);
670 if (D(i) != D.
element(i-1)) ++Count(1,1);
674 for (i = 1; i <= 20; ++i) R(i) = 100 * i + i * i;
676 for (i = 1; i <= 20; ++i)
678 if (R(i) != CR(i)) ++Count(1,1);
679 if (R(i) != CR.
element(i-1)) ++Count(1,1);
680 if (R(i) != R.
element(i-1)) ++Count(1,1);
684 for (i = 1; i <= 20; ++i) C(i) = 100 * i + i * i;
686 for (i = 1; i <= 20; ++i)
688 if (C(i) != CC(i)) ++Count(1,1);
689 if (C(i) != CC.
element(i-1)) ++Count(1,1);
690 if (C(i) != C.
element(i-1)) ++Count(1,1);
748 UT << 3 << 7 << 3 << 9
760 D1 << 6 << 45 << 48 << 20;
784 for (i = 1; i <= 9; ++i)
for (j = 1; j <= 8; ++j)
795 for (i = 1; i <= 8; ++i)
for (j = 1; j <= 8; ++j)
809 for (i = 1; i <= 9; ++i) A.Column(i) =
Helmert(9, i,
true);
814 for (i = 1; i <= 9; ++i) A.Column(i) =
Helmert(9, i);
LogAndSign LogDeterminant(const BaseMatrix &B)
void ReSize(int m, int n, int b)
Real Sum(const BaseMatrix &B)
ReturnMatrix Helmert_transpose(const ColumnVector &Y, bool full=false)
void ReSize(int m, int n, int b)
LogAndSign LogDeterminant() const
void SymmetricBandFunctions(int l1, int l2)
ReturnMatrix Helmert(int n, bool full=false)
GetSubMatrix rows(int, int) const
void LowerBandFunctions(int l1, int l2)
SPMatrix SP(const BaseMatrix &, const BaseMatrix &)
virtual void ReSize(int m, int n)
void Inject(const GeneralMatrix &GM)
static int my_max(int p, int q)
GetSubMatrix column(int) const
virtual void ReSize(int m, int n, int b)
Real SumSquare(const BaseMatrix &B)
void UpperBandFunctions(int u1, int u2)
virtual MatrixBandWidth BandWidth() const
Upper triangular band matrix.
ReturnMatrix sum_columns() const
void Clean(Matrix &A, Real c)
GetSubMatrix Column(int f) const
RowedMatrix AsRow() const
TransposedMatrix t() const
Real Trace(const BaseMatrix &B)
Real MaximumAbsoluteValue(const BaseMatrix &B)
void ReSize(int m, int b)
The usual rectangular matrix.
static int my_min(int p, int q)
Real SumAbsoluteValue(const BaseMatrix &B)
FloatVector FloatVector * a
LU decomposition of a band matrix.
GetSubMatrix Row(int f) const
Lower triangular band matrix.
virtual void resize(int, int)
DiagedMatrix AsDiagonal() const
Real Determinant(const BaseMatrix &B)
void Print(const Matrix &X)
void BandFunctions(int l1, int u1, int l2, int u2)
GetSubMatrix row(int) const
A matrix which can be of any GeneralMatrix type.
void FillWithValues(MultWithCarry &MWC, Matrix &M)