Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/util/newmat/tmt9.cpp @ 4565

Last change on this file since 4565 was 4565, checked in by patrick, 19 years ago

orxonox/trunk: added the newmat library to the project. needs some translation in directory, temp under util/newmat. is needed by the collision detection engine to perform lin alg operations such as eigenvector decomposition. perhaps we will make our own library to do that later.
File size: 7.2 KB
Line 
1//#define WANT_STREAM
2
3#include "include.h"
4#include "newmatap.h"
5
6#include "tmt.h"
7
8#ifdef use_namespace
9using namespace NEWMAT;
10#endif
11
12
13/**************************** test program ******************************/
14
15
16void trymat9()
17{
18   Tracer et("Ninth test of Matrix package");
19   Tracer::PrintTrace();
20
21
22   int i; int j;
23   Matrix A(7,7); Matrix X(7,3);
24   for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
25   for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
26   Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
27   {
28      Tracer et1("Stage 1");
29      Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
30      Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
31      Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
32   }
33   for (i=1;i<=7;i++) D(i,i)=i*i+1;
34   DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
35   {
36      Tracer et1("Stage 2");
37      Matrix DXE = D.i() * X * E;
38      DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE); 
39      E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
40   }
41   DiagonalMatrix F=D;
42   {
43      Tracer et1("Stage 3");
44      F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
45      F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
46   }
47   {
48      Tracer et1("Stage 4");
49      F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
50   }
51   {
52      Tracer et1("Stage 5");
53      // testing equal
54      ColumnVector A(18), B(18);
55      Matrix X(3,3);
56      X << 3 << 5  << 7 << 5 << 8 << 2 << 7 << 2 << 9;
57      SymmetricMatrix S; S << X;
58      B(1) = S == X;         A(1) = true;
59      B(2) = S == (X+1);     A(2) = false;
60      B(3) = (S+2) == (X+2); A(3) = true;
61      Matrix Y = X;
62      B(4) = X == Y;         A(4) = true;
63      B(5) = (X*2) == (Y*2); A(5) = true;
64      Y(3,3) = 10;
65      B(6) = X == Y;         A(6) = false;
66      B(7) = (X*2) == (Y*2); A(7) = false;
67      B(8) = S == Y;         A(8) = false;
68      B(9) = S == S;         A(9) = true;
69      Matrix Z = X.SubMatrix(1,2,2,3);
70      B(10) = X == Z;        A(10) = false;
71      GenericMatrix GS = S;
72      GenericMatrix GX = X;
73      GenericMatrix GY = Y;
74      B(11) = GS == GX;      A(11) = true;
75      B(12) = GS == GY;      A(12) = false;
76      CroutMatrix CS = S;
77      CroutMatrix CX = X;
78      CroutMatrix CY = Y;
79      B(13) = CS == CX;      A(13) = true;
80      B(14) = CS == CY;      A(14) = false;
81      B(15) = X == CX;       A(15) = false;
82      B(16) = X == A;        A(16) = false;
83      B(17) = X == (X | X);  A(17) = false;
84      B(18) = CX == X;       A(18) = false;
85      A = A - B; Print(A);
86   }
87   {
88      Tracer et1("Stage 6");
89      // testing equal
90      ColumnVector A(22), B(22);
91      BandMatrix X(6,2,1);
92      X(1,1)=23; X(1,2)=21;
93      X(2,1)=12; X(2,2)=17; X(2,3)=45;
94      X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
95                 X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
96                            X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
97                                       X(6,4)=49; X(6,5)=41; X(6,6)=17;
98      SymmetricBandMatrix S1(6,2); S1.Inject(X);
99      BandMatrix U(6,2,3); U = 0.0; U.Inject(X);
100      B(1) = U == X;         A(1) = true;
101      B(2) = U == (X*3);     A(2) = false;
102      B(3) = (U*5) == (X*5); A(3) = true;
103      Matrix Y = X;
104      B(4) = X == Y;         A(4) = true;
105      B(5) = (X*2) == (Y*2); A(5) = true;
106      Y(6,6) = 10;
107      B(6) = X == Y;         A(6) = false;
108      B(7) = (X*2) == (Y*2); A(7) = false;
109      B(8) = U == Y;         A(8) = false;
110      B(9) = U == U;         A(9) = true;
111      Matrix Z = X.SubMatrix(1,2,2,3);
112      B(10) = X == Z;        A(10) = false;
113      GenericMatrix GU = U;
114      GenericMatrix GX = X;
115      GenericMatrix GY = Y;
116      B(11) = GU == GX;      A(11) = true;
117      B(12) = GU == GY;      A(12) = false;
118      X = X + X.t(); U = U + U.t();
119      SymmetricBandMatrix S(6,2); S.Inject(X);
120      Matrix D = S-X; Print(D);
121      BandLUMatrix BS = S;
122      BandLUMatrix BX = X;
123      BandLUMatrix BU = U;
124      CroutMatrix CX = X;
125      B(13) = BS == BX;      A(13) = true;
126      B(14) = BX == BU;      A(14) = false;
127      B(15) = X == BX;       A(15) = false;
128      B(16) = X != BX;       A(16) = true;
129      B(17) = BX != BS;      A(17) = false;
130      B(18) = (2*X) != (X*2);A(18) = false;
131      B(19) = (X*2) != (X+2);A(19) = true;
132      B(20) = BX == CX;      A(20) = false;
133      B(21) = CX == BX;      A(21) = false;
134      B(22) = BX == X;       A(22) = false;
135      A = A - B; Print(A);
136      DiagonalMatrix I(6); I=1.0;
137      D = BS.i() * X - I;  Clean(D,0.00000001); Print(D);
138      D = BX.i() * X - I;  Clean(D,0.00000001); Print(D);
139      D = BU.i() * X - I;  Clean(D,0.00000001); Print(D);
140
141      // test row wise load
142      SymmetricBandMatrix X1(6,2);
143      X1.Row(1) << 23;
144      X1.Row(2) << 12 << 17;
145      X1.Row(3) << 35 << 19 << 24;
146      X1.Row(4)       << 17 << 11 << 19;
147      X1.Row(5)             << 10 << 44 << 23;
148      X1.Row(6)                   << 49 << 41 << 17;
149      Matrix M = X1 - S1; Print(M);
150
151      // check out submatrix
152      SymmetricBandMatrix X2(20,3); X2 = 0.0;
153      X2.SubMatrix(2,7,2,7) = X1; X2.SymSubMatrix(11,16) = 2 * X1;
154      Matrix MX1 = X1;
155      Matrix MX2(20,20); MX2 = 0;
156      MX2.SymSubMatrix(2,7) = MX1; MX2.SubMatrix(11,16,11,16) = MX1 * 2;
157      MX2 -= X2; Print(MX2);
158
159      BandMatrix X4(20,3,3); X4 = 0.0;
160      X4.SubMatrix(2,7,3,8) = X1; X4.SubMatrix(11,16,10,15) = 2 * X1;
161      MX1 = X1;
162      Matrix MX4(20,20); MX4 = 0;
163      MX4.SubMatrix(2,7,3,8) = MX1; MX4.SubMatrix(11,16,10,15) = MX1 * 2;
164      MX4 -= X4; Print(MX4);
165
166      MX1 = X1.i() * X1 - IdentityMatrix(6);
167      Clean(MX1,0.00000001); Print(MX1);
168
169   }
170
171   {
172      Tracer et1("Stage 7");
173      // testing equal
174      ColumnVector A(12), B(12);
175      BandMatrix X(6,2,1);
176      X(1,1)=23; X(1,2)=21;
177      X(2,1)=12; X(2,2)=17; X(2,3)=45;
178      X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
179                 X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
180                            X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
181                                       X(6,4)=49; X(6,5)=41; X(6,6)=17;
182      Matrix Y = X;
183      LinearEquationSolver LX = X;
184      LinearEquationSolver LY = Y;
185      CroutMatrix CX = X;
186      CroutMatrix CY = Y;
187      BandLUMatrix BX = X;
188      B(1) = LX == CX;       A(1) = false;
189      B(2) = LY == CY;       A(2) = true;
190      B(3) = X == Y;         A(3) = true;
191      B(4) = BX == LX;       A(4) = true;
192      B(5) = CX == CY;       A(5) = true;
193      B(6) = LX == LY;       A(6) = false;
194      B(7) = BX == BX;       A(7) = true;
195      B(8) = CX == CX;       A(8) = true;
196      B(9) = LX == LX;       A(9) = true;
197      B(10) = LY == LY;      A(10) = true;
198      CroutMatrix CX1 = X.SubMatrix(1,4,1,4);
199      B(11) = CX == CX1;     A(11) = false;
200      BandLUMatrix BX1 = X.SymSubMatrix(1,4);    // error with SubMatrix
201      B(12) = BX == BX1;     A(12) = false;
202      A = A - B; Print(A);
203      DiagonalMatrix I(6); I=1.0; Matrix D;
204      D = LX.i() * X - I;  Clean(D,0.00000001); Print(D);
205      D = LY.i() * X - I;  Clean(D,0.00000001); Print(D);
206      I.ReSize(4); I = 1;
207      D = CX1.i() * X.SymSubMatrix(1,4) - I;  Clean(D,0.00000001); Print(D);
208      D = BX1.i() * X.SubMatrix(1,4,1,4) - I;  Clean(D,0.00000001); Print(D);
209   }
210
211//   cout << "\nEnd of ninth test\n";
212}
213
Note: See TracBrowser for help on using the repository browser.