tmt8.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 #include "include.h"
11 
12 #include "newmatap.h"
13 
14 #include "tmt.h"
15 
16 #ifdef use_namespace
17 using namespace NEWMAT;
18 #endif
19 
20 
21 
22 // **************************** test program ******************************
23 
24 
25 void Transposer(const GenericMatrix& GM1, GenericMatrix&GM2)
26  { GM2 = GM1.t(); }
27 
28 // this is a routine in "Numerical Recipes in C" format
29 // if R is a row vector, C a column vector and D diagonal
30 // make matrix DCR
31 
32 static void DCR(Real d[], Real c[], int m, Real r[], int n, Real **dcr)
33 {
34  int i, j;
35  for (i = 1; i <= m; i++) for (j = 1; j <= n; j++)
36  dcr[i][j] = d[i] * c[i] * r[j];
37 }
38 
39 ReturnMatrix TestReturn(const GeneralMatrix& gm) { return gm; }
40 
41 void trymat8()
42 {
43 // cout << "\nEighth test of Matrix package\n";
44  Tracer et("Eighth test of Matrix package");
46 
47  int i;
48 
49 
50  DiagonalMatrix D(6);
51  for (i=1;i<=6;i++) D(i,i)=i*i+i-10;
52  DiagonalMatrix D2=D;
53  Matrix MD=D;
54 
55  DiagonalMatrix D1(6); for (i=1;i<=6;i++) D1(i,i)=-100+i*i*i;
56  Matrix MD1=D1;
57  Print(Matrix(D*D1-MD*MD1));
58  Print(Matrix((-D)*D1+MD*MD1));
59  Print(Matrix(D*(-D1)+MD*MD1));
60  DiagonalMatrix DX=D;
61  {
62  Tracer et1("Stage 1");
63  DX=(DX+D1)*DX; Print(Matrix(DX-(MD+MD1)*MD));
64  DX=D;
65  DX=-DX*DX+(DX-(-D1))*((-D1)+DX);
66  // Matrix MX = Matrix(MD1);
67  // MD1=DX+(MX.t())*(MX.t()); Print(MD1);
68  MD1=DX+(Matrix(MD1).t())*(Matrix(MD1).t()); Print(MD1);
69  DX=D; DX=DX; DX=D2-DX; Print(DiagonalMatrix(DX));
70  DX=D;
71  }
72  {
73  Tracer et1("Stage 2");
74  D.Release(2);
75  D1=D; D2=D;
76  Print(DiagonalMatrix(D1-DX));
77  Print(DiagonalMatrix(D2-DX));
78  MD1=1.0;
79  Print(Matrix(MD1-1.0));
80  }
81  {
82  Tracer et1("Stage 3");
83  //GenericMatrix
85  LT << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9 << 10;
86  UpperTriangularMatrix UT = LT.t() * 2.0;
87  GenericMatrix GM1 = LT;
88  LowerTriangularMatrix LT1 = GM1-LT; Print(LT1);
89  GenericMatrix GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
90  GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
91  GM2 = GM1*2; LT1 = GM2; LT1 = LT1-LT*2; Print(LT1);
92  GM1.Release();
93  GM1=GM1; LT1=GM1-LT; Print(LT1); LT1=GM1-LT; Print(LT1);
94  GM1.Release();
95  GM1=GM1*4; LT1=GM1-LT*4; Print(LT1);
96  LT1=GM1-LT*4; Print(LT1); GM1.CleanUp();
97  GM1=LT; GM2=UT; GM1=GM1*GM2; Matrix M=GM1; M=M-LT*UT; Print(M);
98  Transposer(LT,GM2); LT1 = LT - GM2.t(); Print(LT1);
99  GM1=LT; Transposer(GM1,GM2); LT1 = LT - GM2.t(); Print(LT1);
100  GM1 = LT; GM1 = GM1 + GM1; LT1 = LT*2-GM1; Print(LT1);
101  DiagonalMatrix D; D << LT; GM1 = D; LT1 = GM1; LT1 -= D; Print(LT1);
102  UpperTriangularMatrix UT1 = GM1; UT1 -= D; Print(UT1);
103  }
104  {
105  Tracer et1("Stage 4");
106  // Another test of SVD
107  Matrix M(12,12); M = 0;
108  M(1,1) = M(2,2) = M(4,4) = M(6,6) =
109  M(7,7) = M(8,8) = M(10,10) = M(12,12) = -1;
110  M(1,6) = M(1,12) = -5.601594;
111  M(3,6) = M(3,12) = -0.000165;
112  M(7,6) = M(7,12) = -0.008294;
113  DiagonalMatrix D;
114  SVD(M,D);
115  SortDescending(D);
116  // answer given by matlab
117  DiagonalMatrix DX(12);
118  DX(1) = 8.0461;
119  DX(2) = DX(3) = DX(4) = DX(5) = DX(6) = DX(7) = 1;
120  DX(8) = 0.1243;
121  DX(9) = DX(10) = DX(11) = DX(12) = 0;
122  D -= DX; Clean(D,0.0001); Print(D);
123  }
124 #ifndef DONT_DO_NRIC
125  {
126  Tracer et1("Stage 5");
127  // test numerical recipes in C interface
128  DiagonalMatrix D(10);
129  D << 1 << 4 << 6 << 2 << 1 << 6 << 4 << 7 << 3 << 1;
130  ColumnVector C(10);
131  C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
132  RowVector R(6);
133  R << 2 << 3 << 5 << 7 << 11 << 13;
134  nricMatrix M(10, 6);
135  DCR( D.nric(), C.nric(), 10, R.nric(), 6, M.nric() );
136  M -= D * C * R; Print(M);
137 
138  D.ReSize(5);
139  D << 1.25 << 4.75 << 9.5 << 1.25 << 3.75;
140  C.ReSize(5);
141  C << 1.5 << 7.5 << 4.25 << 0.0 << 7.25;
142  R.ReSize(9);
143  R << 2.5 << 3.25 << 5.5 << 7 << 11.25 << 13.5 << 0.0 << 1.5 << 3.5;
144  Matrix MX = D * C * R;
145  M.ReSize(MX);
146  DCR( D.nric(), C.nric(), 5, R.nric(), 9, M.nric() );
147  M -= MX; Print(M);
148 
149  // test swap
150  nricMatrix A(3,4); nricMatrix B(4,5);
151  A.Row(1) << 2 << 7 << 3 << 6;
152  A.Row(2) << 6 << 2 << 5 << 9;
153  A.Row(3) << 1 << 0 << 1 << 6;
154  B.Row(1) << 2 << 8 << 4 << 5 << 3;
155  B.Row(2) << 1 << 7 << 5 << 3 << 9;
156  B.Row(3) << 7 << 8 << 2 << 1 << 6;
157  B.Row(4) << 5 << 2 << 9 << 0 << 9;
158  nricMatrix A1(1,2); nricMatrix B1;
159  nricMatrix X(3,5); Matrix X1 = A * B;
160  swap(A, A1); swap(B1, B);
161  for (int i = 1; i <= 3; ++i) for (int j = 1; j <= 5; ++j)
162  {
163  X.nric()[i][j] = 0.0;
164  for (int k = 1; k <= 4; ++k)
165  X.nric()[i][j] += A1.nric()[i][k] * B1.nric()[k][j];
166  }
167  X1 -= X; Print(X1);
168  }
169 #endif
170  {
171  Tracer et1("Stage 6");
172  // test dotproduct
173  DiagonalMatrix test(5); test = 1;
174  ColumnVector C(10);
175  C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
176  RowVector R(10);
177  R << 2 << 3 << 5 << 7 << 11 << 13 << -3 << -4 << 2 << 4;
178  test(1) = (R * C).AsScalar() - DotProduct(C, R);
179  test(2) = C.SumSquare() - DotProduct(C, C);
180  test(3) = 6.0 * (C.t() * R.t()).AsScalar() - DotProduct(2.0 * C, 3.0 * R);
181  Matrix MC = C.AsMatrix(2,5), MR = R.AsMatrix(5,2);
182  test(4) = DotProduct(MC, MR) - (R * C).AsScalar();
183  UpperTriangularMatrix UT(5);
184  UT << 3 << 5 << 2 << 1 << 7
185  << 1 << 1 << 8 << 2
186  << 7 << 0 << 1
187  << 3 << 5
188  << 6;
189  LowerTriangularMatrix LT(5);
190  LT << 5
191  << 2 << 3
192  << 1 << 0 << 7
193  << 9 << 8 << 1 << 2
194  << 0 << 2 << 1 << 9 << 2;
195  test(5) = DotProduct(UT, LT) - Sum(SP(UT, LT));
196  Print(test);
197  // check row-wise load;
198  LowerTriangularMatrix LT1(5);
199  LT1.Row(1) << 5;
200  LT1.Row(2) << 2 << 3;
201  LT1.Row(3) << 1 << 0 << 7;
202  LT1.Row(4) << 9 << 8 << 1 << 2;
203  LT1.Row(5) << 0 << 2 << 1 << 9 << 2;
204  Matrix M = LT1 - LT; Print(M);
205  // check solution with identity matrix
206  IdentityMatrix IM(5); IM *= 2;
207  LinearEquationSolver LES1(IM);
208  LowerTriangularMatrix LTX = LES1.i() * LT;
209  M = LTX * 2 - LT; Print(M);
210  DiagonalMatrix D = IM;
211  LinearEquationSolver LES2(IM);
212  LTX = LES2.i() * LT;
213  M = LTX * 2 - LT; Print(M);
214  UpperTriangularMatrix UTX = LES1.i() * UT;
215  M = UTX * 2 - UT; Print(M);
216  UTX = LES2.i() * UT;
217  M = UTX * 2 - UT; Print(M);
218  }
219 
220  {
221  Tracer et1("Stage 7");
222  // Some more GenericMatrix stuff with *= |= &=
223  // but don't any additional checks
224  BandMatrix BM1(6,2,3);
225  BM1.Row(1) << 3 << 8 << 4 << 1;
226  BM1.Row(2) << 5 << 1 << 9 << 7 << 2;
227  BM1.Row(3) << 1 << 0 << 6 << 3 << 1 << 3;
228  BM1.Row(4) << 4 << 2 << 5 << 2 << 4;
229  BM1.Row(5) << 3 << 3 << 9 << 1;
230  BM1.Row(6) << 4 << 2 << 9;
231  BandMatrix BM2(6,1,1);
232  BM2.Row(1) << 2.5 << 7.5;
233  BM2.Row(2) << 1.5 << 3.0 << 8.5;
234  BM2.Row(3) << 6.0 << 6.5 << 7.0;
235  BM2.Row(4) << 2.5 << 2.0 << 8.0;
236  BM2.Row(5) << 0.5 << 4.5 << 3.5;
237  BM2.Row(6) << 9.5 << 7.5;
238  Matrix RM1 = BM1, RM2 = BM2;
239  Matrix X;
240  GenericMatrix GRM1 = RM1, GBM1 = BM1, GRM2 = RM2, GBM2 = BM2;
241  Matrix Z(6,0); Z = 5; Print(Z);
242  GRM1 |= Z; GBM1 |= Z; GRM2 &= Z.t(); GBM2 &= Z.t();
243  X = GRM1 - BM1; Print(X); X = GBM1 - BM1; Print(X);
244  X = GRM2 - BM2; Print(X); X = GBM2 - BM2; Print(X);
245 
246  GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
247  GRM1 *= GRM2; GBM1 *= GBM2;
248  X = GRM1 - BM1 * BM2; Print(X);
249  X = RM1 * RM2 - GBM1; Print(X);
250 
251  GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
252  GRM1 *= GBM2; GBM1 *= GRM2; // Bs and Rs swapped on LHS
253  X = GRM1 - BM1 * BM2; Print(X);
254  X = RM1 * RM2 - GBM1; Print(X);
255 
256  X = BM1.t(); BandMatrix BM1X = BM1.t();
257  GRM1 = RM1; X -= GRM1.t(); Print(X); X = BM1X - BM1.t(); Print(X);
258 
259  // check that linear equation solver works with Identity Matrix
260  IdentityMatrix IM(6); IM *= 2;
261  GBM1 = BM1; GBM1 *= 4; GRM1 = RM1; GRM1 *= 4;
262  DiagonalMatrix D = IM;
263  LinearEquationSolver LES1(D);
264  BandMatrix BX;
265  BX = LES1.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
266  LinearEquationSolver LES2(IM);
267  BX = LES2.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
268  BX = D.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
269  BX = IM.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
270  BX = IM.i(); BX *= GBM1; BX -= BM1 * 2; X = BX; Print(X);
271 
272  // try symmetric band matrices
273  SymmetricBandMatrix SBM; SBM << SP(BM1, BM1.t());
274  SBM << IM.i() * SBM;
275  X = 2 * SBM - SP(RM1, RM1.t()); Print(X);
276 
277  // Do this again with more general D
278  D << 2.5 << 7.5 << 2 << 5 << 4.5 << 7.5;
279  BX = D.i() * BM1; X = BX - D.i() * RM1;
280  Clean(X,0.00000001); Print(X);
281  BX = D.i(); BX *= BM1; X = BX - D.i() * RM1;
282  Clean(X,0.00000001); Print(X);
283  SBM << SP(BM1, BM1.t());
284  BX = D.i() * SBM; X = BX - D.i() * SP(RM1, RM1.t());
285  Clean(X,0.00000001); Print(X);
286 
287  // test return
288  BX = TestReturn(BM1); X = BX - BM1;
289  if (BX.BandWidth() != BM1.BandWidth()) X = 5;
290  Print(X);
291  }
292 
293 // cout << "\nEnd of eighth test\n";
294 }
295 
296 
297 
void Release()
Definition: newmat.h:514
Real Sum(const BaseMatrix &B)
Definition: newmat.h:2107
Real * nric() const
Definition: newmat.h:1041
void ReSize(int m)
Definition: newmat.h:1034
static void test(int n)
Definition: tmtf.cpp:115
static void DCR(Real d[], Real c[], int m, Real r[], int n, Real **dcr)
Definition: tmt8.cpp:32
SPMatrix SP(const BaseMatrix &, const BaseMatrix &)
Definition: newmat6.cpp:278
void swap(Matrix &A, Matrix &B)
Definition: newmat.h:2159
Real DotProduct(const Matrix &A, const Matrix &B)
Definition: newmat.h:2060
Real * nric() const
Definition: newmat.h:992
double Real
Definition: include.h:307
Real ** nric() const
Definition: newmat.h:741
void ReSize(int m)
Definition: newmat.h:935
Upper triangular matrix.
Definition: newmat.h:799
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
FloatVector * d
Band matrix.
Definition: newmat.h:1096
Real * nric() const
Definition: newmat.h:939
TransposedMatrix t() const
Definition: newmat6.cpp:320
void SVD(const Matrix &, DiagonalMatrix &, Matrix &, Matrix &, bool=true, bool=true)
Definition: svd.cpp:30
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
void SortDescending(GeneralMatrix &gm)
Definition: newmatap.h:139
Real SumSquare() const
Definition: newmat.h:346
Rectangular matrix for use with Numerical Recipes in C.
Definition: newmat.h:711
void Release()
Definition: newmat.h:1421
MatedMatrix AsMatrix(int m, int n) const
Definition: newmat.h:2144
void CleanUp()
Definition: newmat.h:396
Diagonal matrix.
Definition: newmat.h:896
void Transposer(const GenericMatrix &GM1, GenericMatrix &GM2)
Definition: tmt8.cpp:25
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
Symmetric band matrix.
Definition: newmat.h:1245
Row vector.
Definition: newmat.h:953
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
ReturnMatrix TestReturn(const GeneralMatrix &gm)
Definition: tmt8.cpp:39
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
void trymat8()
Definition: tmt8.cpp:41
Identity matrix.
Definition: newmat.h:1350


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