Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/tmt2.cpp @ 4617

Last change on this file since 4617 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: 9.0 KB
Line 
1
2//#define WANT_STREAM
3
4
5#include "include.h"
6
7#include "newmat.h"
8
9#include "tmt.h"
10
11#ifdef use_namespace
12using namespace NEWMAT;
13#endif
14
15
16/**************************** test program ******************************/
17
18
19void trymat2()
20{
21//   cout << "\nSecond test of Matrix package\n\n";
22   Tracer et("Second test of Matrix package");
23   Tracer::PrintTrace();
24
25   int i,j;
26
27   Matrix M(3,5);
28   for (i=1; i<=3; i++) for (j=1; j<=5; j++) M(i,j) = 100*i + j;
29   Matrix X(8,10);
30   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
31   Matrix Y = X; Matrix Z = X;
32   { X.SubMatrix(2,4,3,7) << M; }
33   for (i=1; i<=3; i++) for (j=1; j<=5; j++) Y(i+1,j+2) = 100*i + j;
34   Print(Matrix(X-Y));
35
36
37   Real a[15]; Real* r = a;
38   for (i=1; i<=3; i++) for (j=1; j<=5; j++) *r++ = 100*i + j;
39   { Z.SubMatrix(2,4,3,7) << a; }
40   Print(Matrix(Z-Y));
41
42   { M=33; X.SubMatrix(2,4,3,7) << M; }
43   { Z.SubMatrix(2,4,3,7) = 33; }
44   Print(Matrix(Z-X));
45
46   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
47   Y = X;
48   UpperTriangularMatrix U(5);
49   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
50   { X.SubMatrix(3,7,5,9) << U; }
51   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
52   for (i=1; i<=5; i++) for (j=1; j<i; j++) Y(i+2,j+4) = 0.0;
53   Print(Matrix(X-Y));
54   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
55   Y = X;
56   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
57   { X.SubMatrix(3,7,5,9).Inject(U); }
58   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
59   Print(Matrix(X-Y));
60
61
62   // test growing and shrinking a vector
63   {
64      ColumnVector V(100);
65      for (i=1;i<=100;i++) V(i) = i*i+i;
66      V = V.Rows(1,50);               // to get first 50 vlaues.
67
68      {
69         V.Release(); ColumnVector VX=V;
70         V.ReSize(100); V = 0.0; V.Rows(1,50)=VX;
71      }                               // V now length 100
72
73      M=V; M=100;                     // to make sure V will hold its values
74      for (i=1;i<=50;i++) V(i) -= i*i+i;
75      Print(V);
76
77
78           // test redimensioning vectors with two dimensions given
79      ColumnVector CV1(10); CV1 = 10;
80      ColumnVector CV2(5); CV2.ReSize(10,1); CV2 = 10;
81      V = CV1-CV2; Print(V);
82
83      RowVector RV1(20); RV1 = 100;
84      RowVector RV2; RV2.ReSize(1,20); RV2 = 100;
85      V = (RV1-RV2).t(); Print(V);
86
87      X.ReSize(4,7);
88      for (i=1; i<=4; i++) for (j=1; j<=7; j++) X(i,j) = 1000*i + 10*j;
89      Y = 10.5 * X;
90      Z = 7.25 - Y;
91      M = Z + X * 10.5 - 7.25;
92      Print(M);
93      Y = 2.5 * X;
94      Z = 9.25 + Y;
95      M = Z - X * 2.5 - 9.25;
96      Print(M);
97      U.ReSize(8);
98      for (i=1; i<=8; i++) for (j=i; j<=8; j++) U(i,j) = 100*i + j;
99      Y = 100 - U;
100      M = Y + U - 100;
101      Print(M);
102   }
103
104   {
105      SymmetricMatrix S,T;
106
107      S << (U + U.t());
108      T = 100 - S; M = T + S - 100; Print(M);
109      T = 100 - 2 * S; M = T + S * 2 - 100; Print(M);
110      X = 100 - 2 * S; M = X + S * 2 - 100; Print(M);
111      T = S; T = 100 - T; M = T + S - 100; Print(M);
112   }
113
114   // test new
115   {
116      ColumnVector CV1; RowVector RV1;
117      Matrix* MX; MX = new Matrix; if (!MX) Throw(Bad_alloc("New fails "));
118      MX->ReSize(10,20);
119      for (i = 1; i <= 10; i++) for (j = 1; j <= 20; j++)
120         (*MX)(i,j) = 100 * i + j;
121      ColumnVector* CV = new ColumnVector(10);
122      if (!CV) Throw(Bad_alloc("New fails "));
123      *CV << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9 << 10;
124      RowVector* RV =  new RowVector(CV->t() | (*CV + 10).t());
125      if (!RV) Throw(Bad_alloc("New fails "));
126      CV1 = ColumnVector(10); CV1 = 1; RV1 = RowVector(20); RV1 = 1;
127      *MX -= 100 * *CV * RV1 + CV1 * *RV;
128      Print(*MX);
129      delete MX; delete CV; delete RV;
130   }
131
132
133   // test copying of vectors and matrices with no elements
134   {
135      ColumnVector dims(16);
136      Matrix M1; Matrix M2 = M1; Print(M2);
137      dims(1) = M2.Nrows(); dims(2) = M2.Ncols();
138      dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage();
139      M2 = M1;
140      dims(5) = M2.Nrows(); dims(6) = M2.Ncols();
141      dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage();
142      M2.ReSize(10,20); M2.CleanUp();
143      dims(9) = M2.Nrows(); dims(10) = M2.Ncols();
144      dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage();
145      M2.ReSize(20,10); M2.ReSize(0,0);
146      dims(13) = M2.Nrows(); dims(14) = M2.Ncols();
147      dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage();
148      Print(dims);
149   }
150
151   {
152      ColumnVector dims(16);
153      ColumnVector M1; ColumnVector M2 = M1; Print(M2);
154      dims(1) = M2.Nrows(); dims(2) = M2.Ncols()-1;
155      dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage();
156      M2 = M1;
157      dims(5) = M2.Nrows(); dims(6) = M2.Ncols()-1;
158      dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage();
159      M2.ReSize(10); M2.CleanUp();
160      dims(9) = M2.Nrows(); dims(10) = M2.Ncols()-1;
161      dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage();
162      M2.ReSize(10); M2.ReSize(0);
163      dims(13) = M2.Nrows(); dims(14) = M2.Ncols()-1;
164      dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage();
165      Print(dims);
166   }
167
168   {
169      ColumnVector dims(16);
170      RowVector M1; RowVector M2 = M1; Print(M2);
171      dims(1) = M2.Nrows()-1; dims(2) = M2.Ncols();
172      dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage();
173      M2 = M1;
174      dims(5) = M2.Nrows()-1; dims(6) = M2.Ncols();
175      dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage();
176      M2.ReSize(10); M2.CleanUp();
177      dims(9) = M2.Nrows()-1; dims(10) = M2.Ncols();
178      dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage();
179      M2.ReSize(10); M2.ReSize(0);
180      dims(13) = M2.Nrows()-1; dims(14) = M2.Ncols();
181      dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage();
182      Print(dims);
183   }
184
185   // test identity matrix
186   {
187      Matrix M;
188      IdentityMatrix I(10); DiagonalMatrix D(10); D = 1;
189      M = I; M -= D; Print(M);
190      D -= I; Print(D);
191      ColumnVector X(8);
192      D = 1;
193      X(1) = Sum(D) - Sum(I);
194      X(2) = SumAbsoluteValue(D) - SumAbsoluteValue(I);
195      X(3) = SumSquare(D) - SumSquare(I);
196      X(4) = Trace(D) - Trace(I);
197      X(5) = Maximum(D) - Maximum(I);
198      X(6) = Minimum(D) - Minimum(I);
199      X(7) = LogDeterminant(D).LogValue() - LogDeterminant(I).LogValue();
200      X(8) = LogDeterminant(D).Sign() - LogDeterminant(I).Sign();
201      Clean(X,0.00000001); Print(X);
202
203      for (i = 1; i <= 10; i++) for (j = 1; j <= 10; j++)
204         M(i,j) = 100 * i + j;
205      Matrix N;
206      N = M * I - M; Print(N);
207      N = I * M - M; Print(N);
208      N = M * I.i() - M; Print(N);
209      N = I.i() * M - M; Print(N);
210      N = I.i(); N -= I; Print(N);
211      N = I.t(); N -= I; Print(N);
212      N = I.t(); N += (-I); Print(N); // <----------------
213      D = I; N = D; D = 1; N -= D; Print(N);
214      N = I; D = 1; N -= D; Print(N);
215      N = M + 2 * IdentityMatrix(10); N -= (M + 2 * D); Print(N);
216
217      I *= 4;
218
219      D = 4;
220
221      X.ReSize(14);
222      X(1) = Sum(D) - Sum(I);
223      X(2) = SumAbsoluteValue(D) - SumAbsoluteValue(I);
224      X(3) = SumSquare(D) - SumSquare(I);
225      X(4) = Trace(D) - Trace(I);
226      X(5) = Maximum(D) - Maximum(I);
227      X(6) = Minimum(D) - Minimum(I);
228      X(7) = LogDeterminant(D).LogValue() - LogDeterminant(I).LogValue();  // <--
229      X(8) = LogDeterminant(D).Sign() - LogDeterminant(I).Sign();
230      int i,j;
231      X(9) = I.Maximum1(i) - 4; X(10) = i-1;
232      X(11) = I.Maximum2(i,j) - 4; X(12) = i-10; X(13) = j-10;
233      X(14) = I.Nrows() - 10;
234      Clean(X,0.00000001); Print(X);
235
236
237      N = D.i();
238      N += I / (-16);
239      Print(N);
240      N = M * I - 4 * M; Print(N);
241      N = I * M - 4 * M; Print(N);
242      N = M * I.i() - 0.25 * M; Print(N);
243      N = I.i() * M - 0.25 * M; Print(N);
244      N = I.i(); N -= I * 0.0625; Print(N);
245      N = I.i(); N = N - 0.0625 * I; Print(N);
246      N = I.t(); N -= I; Print(N);
247      D = I * 2; N = D; D = 1; N -= 8 * D; Print(N);
248      N = I * 2; N -= 8 * D; Print(N);
249      N = 0.5 * I + M; N -= M; N -= 2.0 * D; Print(N);
250
251      IdentityMatrix J(10); J = 8;
252      D = 4;
253      DiagonalMatrix E(10); E = 8;
254      N = (I + J) - (D + E); Print(N);
255      N = (5*I + 3*J) - (5*D + 3*E); Print(N);
256      N = (-I + J) - (-D + E); Print(N);
257      N = (I - J) - (D - E); Print(N);
258      N = (I | J) - (D | E); Print(N);
259      N = (I & J) - (D & E); Print(N);
260      N = SP(I,J) - SP(D,E); Print(N);
261      N = D.SubMatrix(2,5,3,8) - I.SubMatrix(2,5,3,8); Print(N);
262
263      N = M; N.Inject(I); D << M; N -= (M + I); N += D; Print(N);
264      D = 4;
265
266      IdentityMatrix K = I.i()*7 - J.t()/4;
267      N = D.i() * 7 - E / 4 - K; Print(N);
268      K = I * J; N = K - D * E; Print(N);
269      N = I * J; N -= D * E; Print(N);
270      K = 5*I - 3*J;
271      N = K - (5*D - 3*E); Print(N);
272      K = I.i(); N = K - 0.0625 * I; Print(N);
273      K = I.t(); N = K - I; Print(N);
274
275
276      K.ReSize(20); D.ReSize(20); D = 1;
277      D -= K; Print(D);
278
279      I.ReSize(3); J.ReSize(3); K = I * J; N = K - I; Print(N);
280      K << D; N = K - D; Print(N);
281
282
283   }
284
285
286//   cout << "\nEnd of second test\n";
287}
Note: See TracBrowser for help on using the repository browser.