tmt1.cpp
Go to the documentation of this file.
1 
6 
7 
8 #define WANT_STREAM
9 
10 
11 
12 #include "include.h"
13 
14 #include "newmat.h"
15 
16 #include "tmt.h"
17 
18 #ifdef use_namespace
19 using namespace NEWMAT;
20 #endif
21 
22 
23 /**************************** test program ******************************/
24 
25 
26 void trymat1()
27 {
28 // cout << "\nFirst test of Matrix package\n\n";
29  Tracer et("First test of Matrix package");
31  {
32  Tracer et1("Stage 1");
33  int i,j;
34 
36  for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
37  SymmetricMatrix S(10);
38  for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
39  SymmetricMatrix S1 = S / 2.0;
40  S = S1 * 2.0;
41  UpperTriangularMatrix U=L.t()*2.0;
42  Print(LowerTriangularMatrix(L-U.t()*0.5));
43  DiagonalMatrix D(10);
44  for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
45  Matrix M=(S+U-D+L)*(L+U-D+S);
46  DiagonalMatrix DD=D*D;
47  LowerTriangularMatrix LD=L*D;
48  // expressions split for Turbo C
49  Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
50  { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
51  { M1 = M1 - U*D + DD - LD + S*S; }
52  { M1 = M1 + U*S - D*S + L*S - 10.0; }
53  M=M1-M;
54  Print(M);
55  }
56  {
57  Tracer et1("Stage 2");
58  int i,j;
59 
61  for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
63  for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
65  for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
67  for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
68  {
69  L=L+LX/0.5; L=L-LX*3.0; L=LX*2.0+L;
70  U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
71  }
72 
73 
74  SymmetricMatrix S(9);
75  for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
76  {
77  SymmetricMatrix S1 = S;
78  S=S1+5.0;
79  S=S-3.0;
80  }
81 
82  DiagonalMatrix D(9);
83  for (i=1;i<=9;i++) D(i,i)=S(i,i);
84  UpperTriangularMatrix U=L.t()*2.0;
85  {
86  U1=U1*2.0 - U; Print(U1);
87  L=L*2.0-D; U=U-D;
88  }
89  Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
90  }
91  {
92  Tracer et1("Stage 3");
93  int i,j;
94  Matrix M(10,3), N(10,3);
95  for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
96  { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
97  Matrix MN = M + N, M1;
98 
99  M1 = M; M1 += N; M1 -= MN; Print(M1);
100  M1 = M; M1 += M1; M1 = M1 - M * 2; Print(M1);
101  M1 = M; M1 += N * 2; M1 -= (MN + N); Print(M1);
102  M1 = M; M1 -= M1; Print(M1);
103  M1 = M; M1 -= MN + M1; M1 += N + M; Print(M1);
104  M1 = M; M1 -= 5; M1 -= M; M1 *= 0.2; M1 = M1 + 1; Print(M1);
105  Matrix NT = N.t();
106  M1 = M; M1 *= NT; M1 -= M * N.t(); Print(M1);
107  M = M * M.t();
108  DiagonalMatrix D(10); D = 2;
109  M1 = M; M1 += D; M1 -= M; M1 = M1 - D; Print(M1);
110  M1 = M; M1 -= D; M1 -= M; M1 = M1 + D; Print(M1);
111  M1 = M; M1 *= D; M1 /= 2; M1 -= M; Print(M1);
112  SymmetricMatrix SM; SM << M;
113  // UpperTriangularMatrix SM; SM << M;
114  SM += 10; M1 = SM - M; M1 /=10; M1 = M1 - 1; Print(M1);
115  }
116  {
117  Tracer et1("Stage 4");
118  int i,j;
119  Matrix M(10,3), N(10,5);
120  for (i = 1; i<=10; i++) for (j = 1; j<=3; j++) M(i,j) = 2*i-j;
121  for (i = 1; i<=10; i++) for (j = 1; j<=5; j++) N(i,j) = i*j + 20;
122  Matrix M1;
123  M1 = M; M1 |= N; M1 &= N | M;
124  M1 -= (M | N) & (N | M); Print(M1);
125  M1 = M; M1 |= M1; M1 &= M1;
126  M1 -= (M | M) & (M | M); Print(M1);
127 
128  }
129  {
130  Tracer et1("Stage 5");
131  int i,j;
132  BandMatrix BM1(10,2,3), BM2(10,4,1); Matrix M1(10,10), M2(10,10);
133  for (i=1;i<=10;i++) for (j=1;j<=10;j++)
134  { M1(i,j) = 0.5*i+j*j-50; M2(i,j) = (i*101 + j*103) % 13; }
135  BM1.Inject(M1); BM2.Inject(M2);
136  BandMatrix BM = BM1; BM += BM2;
137  Matrix M1X = BM1; Matrix M2X = BM2; Matrix MX = BM;
138  MX -= M1X + M2X; Print(MX);
139  MX = BM1; MX += BM2; MX -= M1X; MX -= M2X; Print(MX);
140  SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t();
141  SymmetricBandMatrix SM2; SM2 << BM2 * BM2.t();
142  SM1 *= 5.5;
143  M1X *= M1X.t(); M1X *= 5.5; M2X *= M2X.t();
144  SM1 -= SM2; M1 = SM1 - M1X + M2X; Print(M1);
145  M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1);
146  M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1);
147  M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1);
148 
149  }
150  {
151  Tracer et1("Stage 6");
152  int i,j;
153  Matrix M(10,10), N(10,10);
154  for (i = 1; i<=10; i++) for (j = 1; j<=10; j++)
155  { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
156  GenericMatrix GM = M;
157  GM += N; Matrix M1 = GM - N - M; Print(M1);
158  DiagonalMatrix D(10); D = 3;
159  GM = D; GM += N; GM += M; GM += D;
160  M1 = D*2 - GM + M + N; Print(M1);
161  GM = D; GM *= 4; GM += 16; GM /= 8; GM -= 2;
162  GM -= D / 2; M1 = GM; Print(M1);
163  GM = D; GM *= M; GM *= N; GM /= 3; M1 = M*N - GM; Print(M1);
164  GM = D; GM |= M; GM &= N | D; M1 = GM - ((D | M) & (N | D));
165  Print(M1);
166  GM = M; M1 = M; GM += 5; GM *= 3; M *= 3; M += 15; M1 = GM - M;
167  Print(M1);
168  D.ReSize(10); for (i = 1; i<=10; i++) D(i) = i;
169  M1 = D + 10; GM = D; GM += 10; M1 -= GM; Print(M1);
170  GM = M; GM -= D; M1 = GM; GM = D; GM -= M; M1 += GM; Print(M1);
171  GM = M; GM *= D; M1 = GM; GM = D; GM *= M.t();
172  M1 -= GM.t(); Print(M1);
173  GM = M; GM += 2 * GM; GM -= 3 * M; M1 = GM; Print(M1);
174  GM = M; GM |= GM; GM -= (M | M); M1 = GM; Print(M1);
175  GM = M; GM &= GM; GM -= (M & M); M1 = GM; Print(M1);
176  M1 = M; M1 = (M1.t() & M.t()) - (M | M).t(); Print(M1);
177  M1 = M; M1 = (M1.t() | M.t()) - (M & M).t(); Print(M1);
178 
179  }
180 
181  {
182  Tracer et1("Stage 7");
183  // test for bug in MS VC5
184  int n = 3;
185  int k; int j;
186  Matrix A(n,n), B(n,n);
187 
188  //first version - MS VC++ 5 mis-compiles if optimisation is on
189  for (k=1; k<=n; k++)
190  {
191  for (j = 1; j <= n; j++) A(k,j) = ((k-1) * (2*j-1));
192  }
193 
194  //second version
195  for (k=1; k<=n; k++)
196  {
197  const int k1 = k-1; // otherwise Visual C++ 5 fails
198  for (j = 1; j <= n; j++) B(k,j) = (k1 * (2*j-1));
199  }
200 
201  if (A != B)
202  {
203  cout << "\nVisual C++ version 5 compiler error?";
204  cout << "\nTurn off optimisation";
205  }
206 
207  A -= B; Print(A);
208 
209  }
210 
211  {
212  Tracer et1("Stage 8");
213  // Cross product
214  ColumnVector i(3); i << 1 << 0 << 0;
215  ColumnVector j(3); j << 0 << 1 << 0;
216  ColumnVector k(3); k << 0 << 0 << 1;
217  ColumnVector X;
218  X = CrossProduct(i,j) - k; Print(X);
219  X = CrossProduct(j,k) - i; Print(X);
220  X = CrossProduct(k,i) - j; Print(X);
221  X = CrossProduct(j,i) + k; Print(X);
222  X = CrossProduct(k,j) + i; Print(X);
223  X = CrossProduct(i,k) + j; Print(X);
224  X = CrossProduct(i,i); Print(X);
225  X = CrossProduct(j,j); Print(X);
226  X = CrossProduct(k,k); Print(X);
227 
228  ColumnVector A(3); A << 2.25 << 1.75 << -7.5;
229  ColumnVector B(3); B << -0.5 << 4.75 << 3.25;
230  ColumnVector C(3); C << 9.25 << 3.5 << 1.25;
231 
232  Real d0 = (A | B | C).Determinant(); // Vector triple product
233  Real d1 = DotProduct(CrossProduct(A, B), C);
234  Real d2 = DotProduct(CrossProduct(B, C), A);
235  Real d3 = DotProduct(CrossProduct(C, A), B);
236  X << (d1 - d0) << (d2 - d0) << (d3 - d0);
237  Clean(X, 0.000000001); Print(X);
238 
239  X = CrossProduct(A, CrossProduct(B, C))
240  + CrossProduct(B, CrossProduct(C, A))
241  + CrossProduct(C, CrossProduct(A, B));
242  Print(X);
243 
244  RowVector XT = CrossProduct(A.AsRow(), B.AsRow());
245  XT -= CrossProduct(A, B).AsRow();
246  Print(XT);
247  }
248 
249  {
250  Tracer et1("Stage 9");
251  // More cross product
252  int i, j;
253  Matrix M(10,3), N(10,3);
254  for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
255  { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
256 
257  Matrix CP1 = CrossProductRows(M, N);
258  Matrix CP2(10,3);
259  for (i = 1; i<=10; i++)
260  CP2.Row(i) = CrossProduct(M.Row(i), N.Row(i));
261  CP2 -= CP1; Print(CP2);
262 
263  CP2 = CrossProductColumns(M.t(), N.t());
264  CP2 -= CP1.t(); Print(CP2);
265  }
266 
267  {
268  Tracer et1("Stage 10");
269  // Make sure RNG works
270  MultWithCarry mwc;
271  ColumnVector cv(10);
272  for (int i = 1; i <= 10; ++i) cv(i) = mwc.Next();
273  cv *= 100.0;
274  cv(1) -= 6.27874;
275  cv(2) -= 42.1718;
276  cv(3) -= 80.2854;
277  cv(4) -= 12.961;
278  cv(5) -= 17.7499;
279  cv(6) -= 13.2657;
280  cv(7) -= 50.4923;
281  cv(8) -= 26.095;
282  cv(9) -= 57.9147;
283  cv(10) -= 30.1778;
284  Clean(cv, 0.0001); Print(cv);
285  }
286 
287 
288 // cout << "\nEnd of first test\n";
289 }
290 
291 
ReturnMatrix CrossProductColumns(const Matrix &A, const Matrix &B)
Definition: newmat.h:2066
void trymat1()
Definition: tmt1.cpp:26
Real DotProduct(const Matrix &A, const Matrix &B)
Definition: newmat.h:2060
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
double Real
Definition: include.h:307
void ReSize(int m)
Definition: newmat.h:935
Upper triangular matrix.
Definition: newmat.h:799
Real Next()
Definition: tmt.cpp:187
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
RowedMatrix AsRow() const
Definition: newmat.h:2141
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
Matrix CrossProduct(const Matrix &A, const Matrix &B)
Definition: newmat.h:2062
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
Symmetric band matrix.
Definition: newmat.h:1245
Real Determinant(const BaseMatrix &B)
Definition: newmat.h:2092
Row vector.
Definition: newmat.h:953
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
ReturnMatrix CrossProductRows(const Matrix &A, const Matrix &B)
Definition: newmat.h:2064
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jun 7 2019 22:06:45