1 /// \ingroup newmat
2 ///@{
3 
4 /// \file tmtg.cpp
5 /// Part of matrix library test program.
6 
7 
8 //#define WANT_STREAM
9 
10 #define WANT_MATH                   // for sqrt
11 
12 #include "include.h"
13 
14 #include "newmatap.h"
15 
16 #include "tmt.h"
17 
18 #ifdef use_namespace
19 using namespace NEWMAT;
20 #endif
21 
22 
23 
trymatg()24 void trymatg()
25 {
26 //   cout << "\nSixteenth test of Matrix package\n";
27 //   cout << "\n";
28    Tracer et("Sixteenth test of Matrix package");
29    Tracer::PrintTrace();
30 
31    int i,j;
32    Matrix M(4,7);
33    for (i=1; i<=4; i++) for (j=1; j<=7; j++) M(i,j) = 100 * i + j;
34    ColumnVector CV = M.AsColumn();
35    {
36       Tracer et1("Stage 1");
37       RowVector test(7);
38       test(1) = SumSquare(M);
39       test(2) = SumSquare(CV);
40       test(3) = SumSquare(CV.t());
41       test(4) = SumSquare(CV.AsDiagonal());
42       test(5) = SumSquare(M.AsColumn());
43       test(6) = Matrix(CV.t()*CV)(1,1);
44       test(7) = (CV.t()*CV).AsScalar();
45       test = test - 2156560.0; Print(test);
46    }
47 
48    UpperTriangularMatrix U(6);
49    for (i=1; i<=6; i++) for (j=i; j<=6; j++) U(i,j) = i + (i-j) * (i-j);
50    M = U; DiagonalMatrix D; D << U;
51    LowerTriangularMatrix L = U.t(); SymmetricMatrix S; S << (L+U)/2.0;
52    {
53       Tracer et1("Stage 2");
54       RowVector test(6);
55       test(1) = U.Trace();
56       test(2) = L.Trace();
57       test(3) = D.Trace();
58       test(4) = S.Trace();
59       test(5) = M.Trace();
60       test(6) = ((L+U)/2.0).Trace();
61       test = test - 21; Print(test);
62       test(1) = LogDeterminant(U).Value();
63       test(2) = LogDeterminant(L).Value();
64       test(3) = LogDeterminant(D).Value();
65       test(4) = LogDeterminant(D).Value();
66       test(5) = LogDeterminant((L+D)/2.0).Value();
67       test(6) = Determinant((L+D)/2.0);
68       test = test - 720; Clean(test,0.000000001); Print(test);
69    }
70 
71    {
72       Tracer et1("Stage 3");
73       S << L*U; M = S;
74       RowVector test(8);
75       test(1) = LogDeterminant(S).Value();
76       test(2) = LogDeterminant(M).Value();
77       test(3) = LogDeterminant(L*U).Value();
78       test(4) = LogDeterminant(Matrix(L*L)).Value();
79       test(5) = Determinant(S);
80       test(6) = Determinant(M);
81       test(7) = Determinant(L*U);
82       test(8) = Determinant(Matrix(L*L));
83       test = test - 720.0 * 720.0; Clean(test,0.00000001); Print(test);
84    }
85 
86    {
87       Tracer et1("Stage 4");
88       M = S * S;
89       Matrix SX = S;
90       RowVector test(3);
91       test(1) = SumSquare(S);
92       test(2) = SumSquare(SX);
93       test(3) = Trace(M);
94 		test = test - 3925961.0; Print(test);
95    }
96    {
97       Tracer et1("Stage 5");
98       SymmetricMatrix SM(10), SN(10);
99       Real S = 0.0;
100       for (i=1; i<=10; i++) for (j=i; j<=10; j++)
101       {
102          SM(i,j) = 1.5 * i - j; SN(i,j) = SM(i,j) * SM(i,j);
103          if (SM(i,j) < 0.0)   SN(i,j) = - SN(i,j);
104          S += SN(i,j) * ((i==j) ? 1.0 : 2.0);
105       }
106       Matrix M = SM, N = SN; RowVector test(4);
107       test(1) = SumAbsoluteValue(SN);
108       test(2) = SumAbsoluteValue(-SN);
109       test(3) = SumAbsoluteValue(N);
110       test(4) = SumAbsoluteValue(-N);
111       test = test - 1168.75; Print(test);
112       test(1) = Sum(SN);
113       test(2) = -Sum(-SN);
114       test(3) = Sum(N);
115       test(4) = -Sum(-N);
116       test = test - S; Print(test);
117       test(1) = MaximumAbsoluteValue(SM);
118       test(2) = MaximumAbsoluteValue(-SM);
119       test(3) = MaximumAbsoluteValue(M);
120       test(4) = MaximumAbsoluteValue(-M);
121       test = test - 8.5; Print(test);
122    }
123    {
124       Tracer et1("Stage 6");
125       Matrix M(15,20); Real value = 0.0;
126       for (i=1; i<=15; i++) { for (j=1; j<=20; j++) M(i,j) = 1.5 * i - j; }
127       for (i=1; i<=20; i++)
128       { Real v = SumAbsoluteValue(M.Column(i)); if (value<v) value = v; }
129       RowVector test(3);
130       test(1) = value;
131       test(2) = Norm1(M);
132       test(3) = NormInfinity(M.t());
133       test = test - 165; Print(test);
134       test.ReSize(5);
135       ColumnVector CV = M.AsColumn(); M = CV.t();
136       test(1) = Norm1(CV.t());
137       test(2) = MaximumAbsoluteValue(M);
138       test(3) = NormInfinity(CV);
139       test(4) = Norm1(M);
140       test(5) = NormInfinity(M.t());
141       test = test - 21.5; Print(test);
142    }
143    {
144       Tracer et1("Stage 7");
145       Matrix M(15,20);
146       for (i=1; i<=15; i++) { for (j=1; j<=20; j++) M(i,j) = 2.5 * i - j; }
147       SymmetricMatrix SM; SM << M * M.t();
148       ColumnVector test(6);
149       test(1) = sqrt(SumSquare(M)) - NormFrobenius(M);
150       Real a = sqrt(SumSquare(SM));
151       test(2) = NormFrobenius(M * M.t()) - a;
152       test(3) = SM.NormFrobenius() - a;
153       test(4) = (M * M.t()).NormFrobenius() - a;
154       test(5) = NormFrobenius(SM) - a;
155       test(6) = Trace(SM) - M.SumSquare();
156       Clean(test, 0.00000001); Print(test);
157   }
158 
159 //   cout << "\nEnd of Sixteenth test\n";
160 }
161 
162 
163 ///@}
164