Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/tmtm.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: 6.3 KB
Line 
1
2#define WANT_STREAM
3
4#define WANT_MATH
5
6#include "newmat.h"
7
8#include "tmt.h"
9
10#ifdef use_namespace
11using namespace NEWMAT;
12#endif
13
14
15
16// test Kronecker Product
17
18
19void trymatm()
20{
21   Tracer et("Twenty second test of Matrix package");
22   Tracer::PrintTrace();
23
24   {
25      Tracer et1("Stage 1");
26
27
28      Matrix A(2,3);
29      A << 3 << 5 << 2
30        << 4 << 1 << 6;
31
32      Matrix B(4,3);
33      B <<  7 <<  2 <<  9
34        <<  1 <<  3 <<  6
35        <<  4 << 10 <<  5
36        << 11 <<  8 << 12;
37
38      Matrix C(8, 9);
39
40      C.Row(1) << 21 <<  6 << 27  << 35 << 10 << 45  << 14 <<  4 << 18;
41      C.Row(2) <<  3 <<  9 << 18  <<  5 << 15 << 30  <<  2 <<  6 << 12;
42      C.Row(3) << 12 << 30 << 15  << 20 << 50 << 25  <<  8 << 20 << 10;
43      C.Row(4) << 33 << 24 << 36  << 55 << 40 << 60  << 22 << 16 << 24;
44
45      C.Row(5) << 28 <<  8 << 36  <<  7 <<  2 <<  9  << 42 << 12 << 54;
46      C.Row(6) <<  4 << 12 << 24  <<  1 <<  3 <<  6  <<  6 << 18 << 36;
47      C.Row(7) << 16 << 40 << 20  <<  4 << 10 <<  5  << 24 << 60 << 30;
48      C.Row(8) << 44 << 32 << 48  << 11 <<  8 << 12  << 66 << 48 << 72;
49
50      Matrix AB = KP(A,B) - C; Print(AB);
51
52      IdentityMatrix I1(10); IdentityMatrix I2(15); I2 *= 2;
53      DiagonalMatrix D = KP(I1, I2) - IdentityMatrix(150) * 2;
54      Print(D);
55   }
56
57   {
58      Tracer et1("Stage 2");
59
60      UpperTriangularMatrix A(3);
61      A << 3 << 8 << 5
62             << 7 << 2
63                  << 4;
64      UpperTriangularMatrix B(4);
65      B << 4 << 1 << 7 << 2
66             << 3 << 9 << 8
67                  << 1 << 5
68                       << 6;
69
70      UpperTriangularMatrix C(12);
71
72      C.Row(1) <<12<< 3<<21<< 6 <<32<< 8<<56<<16 <<20<< 5<<35<<10;
73      C.Row(2)     << 9<<27<<24 << 0<<24<<72<<64 << 0<<15<<45<<40;
74      C.Row(3)         << 3<<15 << 0<< 0<< 8<<40 << 0<< 0<< 5<<25;
75      C.Row(4)             <<18 << 0<< 0<< 0<<48 << 0<< 0<< 0<<30;
76
77      C.Row(5)                  <<28<< 7<<49<<14 << 8<< 2<<14<< 4;
78      C.Row(6)                      <<21<<63<<56 << 0<< 6<<18<<16;
79      C.Row(7)                          << 7<<35 << 0<< 0<< 2<<10;
80      C.Row(8)                              <<42 << 0<< 0<< 0<<12;
81
82      C.Row(9)                                   <<16<< 4<<28<< 8;
83      C.Row(10)                                      <<12<<36<<32;
84      C.Row(11)                                          << 4<<20;
85      C.Row(12)                                              <<24;
86
87
88      UpperTriangularMatrix AB = KP(A,B) - C; Print(AB);
89
90      LowerTriangularMatrix BT = B.t(); Matrix N(12,12);
91
92      N.Row(1) <<12 << 0<< 0<< 0 <<32<< 0<< 0<< 0 <<20<< 0<< 0<< 0;
93      N.Row(2) << 3 << 9<< 0<< 0 << 8<<24<< 0<< 0 << 5<<15<< 0<< 0;
94      N.Row(3) <<21 <<27<< 3<< 0 <<56<<72<< 8<< 0 <<35<<45<< 5<< 0;
95      N.Row(4) << 6 <<24<<15<<18 <<16<<64<<40<<48 <<10<<40<<25<<30;
96
97      N.Row(5) << 0 << 0<< 0<< 0 <<28<< 0<< 0<< 0 << 8<< 0<< 0<< 0;
98      N.Row(6) << 0 << 0<< 0<< 0 << 7<<21<< 0<< 0 << 2<< 6<< 0<< 0;
99      N.Row(7) << 0 << 0<< 0<< 0 <<49<<63<< 7<< 0 <<14<<18<< 2<< 0;
100      N.Row(8) << 0 << 0<< 0<< 0 <<14<<56<<35<<42 << 4<<16<<10<<12;
101
102      N.Row(9) << 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<16<< 0<< 0<< 0;
103      N.Row(10)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 4<<12<< 0<< 0;
104      N.Row(11)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<28<<36<< 4<< 0;
105      N.Row(12)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 8<<32<<20<<24;
106
107      Matrix N1 = KP(A, BT); N1 -= N; Print(N1);
108      AB << KP(A, BT); AB << (AB - N); Print(AB);
109      BT << KP(A, BT); BT << (BT - N); Print(BT);
110
111      LowerTriangularMatrix AT = A.t();
112      N1 = KP(AT, B); N1 -= N.t(); Print(N1);
113      AB << KP(AT, B); AB << (AB - N.t()); Print(AB);
114      BT << KP(AT, B); BT << (BT - N.t()); Print(BT);
115   }
116
117   {
118      Tracer et1("Stage 3");
119
120      BandMatrix BMA(6,2,3);
121      BMA.Row(1) << 5.25 << 4.75 << 2.25 << 1.75;
122      BMA.Row(2) << 1.25 << 9.75 << 4.50 << 0.25 << 1.50;
123      BMA.Row(3) << 7.75 << 1.50 << 3.00 << 4.25 << 0.50 << 5.50;
124      BMA.Row(4) << 2.75 << 9.00 << 8.00 << 3.25 << 3.50;
125      BMA.Row(5) << 8.75 << 6.25 << 5.00 << 5.75;
126      BMA.Row(6) << 3.75 << 6.75 << 6.00;
127
128      Matrix A = BMA;
129
130      BandMatrix BMB(4,2,1);
131      BMB.Row(1) << 4.5 << 9.5;
132      BMB.Row(2) << 1.5 << 6.0 << 2.0;
133      BMB.Row(3) << 0.5 << 2.5 << 8.5 << 7.5;
134      BMB.Row(4) << 3.0 << 4.0 << 6.5;
135
136      Matrix B = BMB;
137
138      BandMatrix BMC = KP(BMA, BMB);
139      BandMatrix BMC1(24,11,15);
140      BMC1.Inject(Matrix(KP(BMA, B)));  // not directly Band Matrix
141      Matrix C2 = KP(A, BMB);
142      Matrix C = KP(A, B);
143
144      Matrix M = C - BMC; Print(M);
145      M = C - BMC1; Print(M);
146      M = C - C2; Print(M);
147
148      RowVector X(4);
149      X(1) = BMC.BandWidth().Lower() - 10;
150      X(2) = BMC.BandWidth().Upper() - 13;
151      X(3) = BMC1.BandWidth().Lower() - 11;
152      X(4) = BMC1.BandWidth().Upper() - 15;
153      Print(X);
154
155      UpperTriangularMatrix UT;  UT << KP(BMA, BMB);
156      UpperTriangularMatrix UT1; UT1 << (C - UT); Print(UT1);
157      LowerTriangularMatrix LT;  LT << KP(BMA, BMB);
158      LowerTriangularMatrix LT1; LT1 << (C - LT); Print(LT1);
159   }
160
161   {
162      Tracer et1("Stage 4");
163
164      SymmetricMatrix SM1(4);
165      SM1.Row(1) << 2;
166      SM1.Row(2) << 4 << 5;
167      SM1.Row(3) << 9 << 2 << 1;
168      SM1.Row(4) << 3 << 6 << 8 << 2;
169
170      SymmetricMatrix SM2(3);
171      SM2.Row(1) <<  3;
172      SM2.Row(2) << -7 << -6;
173      SM2.Row(3) <<  4 << -2 << -1;
174
175      SymmetricMatrix SM = KP(SM1, SM2);
176      Matrix M1 = SM1; Matrix M2 = SM2;
177      Matrix M = KP(SM1, SM2); M -= SM; Print(M);
178      M = KP(SM1, SM2) - SM; Print(M);
179      M = KP(M1, SM2) - SM; Print(M);
180      M = KP(SM1, M2) - SM; Print(M);
181      M = KP(M1, M2); M -= SM; Print(M);
182   }
183
184   {
185      Tracer et1("Stage 5");
186
187      Matrix A(2,3);
188      A << 3 << 5 << 2
189        << 4 << 1 << 6;
190
191      Matrix B(3,4);
192      B <<  7 <<  2 <<  9 << 11
193        <<  1 <<  3 <<  6 <<  8
194        <<  4 << 10 <<  5 << 12;
195
196      RowVector C(2); C << 3 << 7;
197      ColumnVector D(4); D << 0 << 5 << 13 << 11;
198
199      Matrix M = KP(C * A, B * D) - KP(C, B) * KP(A, D); Print(M);
200   }
201
202   {
203      Tracer et1("Stage 6");
204
205      RowVector A(3), B(5), C(15);
206      A << 5 << 2 << 4;
207      B << 3 << 2 << 0 << 1 << 6;
208      C << 15 << 10 << 0 << 5 << 30
209        <<  6 <<  4 << 0 << 2 << 12
210        << 12 <<  8 << 0 << 4 << 24;
211      Matrix N = KP(A, B) - C;    Print(N);
212      N = KP(A.t(), B.t()) - C.t();    Print(N);
213      N = KP(A.AsDiagonal(), B.AsDiagonal()) - C.AsDiagonal();    Print(N);
214   }
215
216}
217
218
219
220
221
222
Note: See TracBrowser for help on using the repository browser.