tmtc.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 
11 #include "include.h"
12 #include "newmatap.h"
13 
14 #include "tmt.h"
15 
16 #ifdef use_namespace
17 using namespace NEWMAT;
18 #endif
19 
20 
21 
22 
23 void trymatc()
24 {
25 // cout << "\nTwelfth test of Matrix package\n";
26  Tracer et("Twelfth test of Matrix package");
27  Tracer::PrintTrace();
28  DiagonalMatrix D(15); D=1.5;
29  Matrix A(15,15);
30  int i,j;
31  for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
32  { A = A + D; }
33  ColumnVector B(15);
34  for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
35  {
36  Tracer et1("Stage 1");
37  ColumnVector B1=B;
38  B=(A*2.0).i() * B1;
39  Matrix X = A*B-B1/2.0;
40  Clean(X, 0.000000001); Print(X);
41  A.ReSize(3,5);
42  for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
43 
44  B = A.AsColumn()+10000;
45  RowVector R = (A+10000).AsColumn().t();
46  Print( RowVector(R-B.t()) );
47  }
48 
49  {
50  Tracer et1("Stage 2");
51  B = A.AsColumn()+10000;
52  Matrix XR = (A+10000).AsMatrix(15,1).t();
53  Print( RowVector(XR-B.t()) );
54  }
55 
56  {
57  Tracer et1("Stage 3");
58  B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
59  Matrix MR = (A+10000).AsColumn().t();
60  Print( RowVector(MR-B.t()) );
61 
62  B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
63  MR = A.AsColumn().t();
64  Print( RowVector(MR-B.t()) );
65  }
66 
67  {
68  Tracer et1("Stage 4");
69  B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
70  RowVector R = A.AsColumn().t();
71  Print( RowVector(R-B.t()) );
72  }
73 
74  {
75  Tracer et1("Stage 5");
76  RowVector R = (A.AsColumn()-5000).t();
77  B = ((R.t()+10000) - A.AsColumn())-5000;
78  Print( RowVector(B.t()) );
79  }
80 
81  {
82  Tracer et1("Stage 6");
83  B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
84  Print(ColumnVector(B1-B));
85  }
86 
87  {
88  Tracer et1("Stage 7");
89  Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
90  for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
91  Print(B);
92  }
93 
94  {
95  Tracer et1("Stage 8");
96  A.ReSize(7,7); D.ReSize(7);
97  for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
98  for (i=1; i<=7; i++) D(i,i) = i;
99  UpperTriangularMatrix U; U << A;
100  Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
101  A.Inject(D); Print(Matrix(X-A));
102  X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
103  Print(Matrix(X-A));
104  }
105 
106  {
107  Tracer et1("Stage 9");
108  A.ReSize(7,5);
109  for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
110  Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
111  Matrix X = A;
112  Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
113  Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
114  }
115 
116  {
117  Tracer et1("Stage 10");
118  // some tests on submatrices
119  UpperTriangularMatrix U(20);
120  for (i=1; i<=20; i++) for (j=i; j<=20; j++) U(i,j)=100 * i + j;
121  UpperTriangularMatrix V = U.SymSubMatrix(1,5);
122  UpperTriangularMatrix U1 = U;
123  U1.SubMatrix(4,8,5,9) /= 2;
124  U1.SubMatrix(4,8,5,9) += 388 * V;
125  U1.SubMatrix(4,8,5,9) *= 2;
126  U1.SubMatrix(4,8,5,9) += V;
127  U1 -= U; UpperTriangularMatrix U2 = U1;
128  U1 << U1.SubMatrix(4,8,5,9);
129  U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
130  U1 -= (777*V); Print(U1);
131 
132  U1 = U; U1.SubMatrix(4,8,5,9) -= U.SymSubMatrix(1,5);
133  U1 -= U; U2 = U1; U1 << U1.SubMatrix(4,8,5,9);
134  U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
135  U1 += V; Print(U1);
136 
137  U1 = U;
138  U1.SubMatrix(3,10,15,19) += 29;
139  U1 -= U;
140  Matrix X = U1.SubMatrix(3,10,15,19); X -= 29; Print(X);
141  U1.SubMatrix(3,10,15,19) *= 0; Print(U1);
142 
143  LowerTriangularMatrix L = U.t();
144  LowerTriangularMatrix M = L.SymSubMatrix(1,5);
145  LowerTriangularMatrix L1 = L;
146  L1.SubMatrix(5,9,4,8) /= 2;
147  L1.SubMatrix(5,9,4,8) += 388 * M;
148  L1.SubMatrix(5,9,4,8) *= 2;
149  L1.SubMatrix(5,9,4,8) += M;
150  L1 -= L; LowerTriangularMatrix L2 = L1;
151  L1 << L1.SubMatrix(5,9,4,8);
152  L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
153  L1 -= (777*M); Print(L1);
154 
155  L1 = L; L1.SubMatrix(5,9,4,8) -= L.SymSubMatrix(1,5);
156  L1 -= L; L2 =L1; L1 << L1.SubMatrix(5,9,4,8);
157  L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
158  L1 += M; Print(L1);
159 
160  L1 = L;
161  L1.SubMatrix(15,19,3,10) -= 29;
162  L1 -= L;
163  X = L1.SubMatrix(15,19,3,10); X += 29; Print(X);
164  L1.SubMatrix(15,19,3,10) *= 0; Print(L1);
165  }
166 
167  {
168  Tracer et1("Stage 11");
169  // more tests on submatrices
170  Matrix M(20,30);
171  for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
172  Matrix M1 = M;
173 
174  for (j=1; j<=30; j++)
175  { ColumnVector CV = 3 * M1.Column(j); M.Column(j) += CV; }
176  for (i=1; i<=20; i++)
177  { RowVector RV = 5 * M1.Row(i); M.Row(i) -= RV; }
178 
179  M += M1; Print(M);
180 
181  }
182 
183  {
184  Tracer et1("Stage 12");
185  // more tests on Release
186  Matrix M(20,30);
187  for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
188  Matrix M1 = M;
189  M.Release();
190  Matrix M2 = M;
191  Matrix X = M; Print(X);
192  X = M1 - M2; Print(X);
193 
194 #ifndef DONT_DO_NRIC
195  nricMatrix N = M1;
196  nricMatrix N1 = N;
197  N.Release();
198  nricMatrix N2 = N;
199  nricMatrix Y = N; Print(Y);
200  Y = N1 - N2; Print(Y);
201 
202  N = M1 / 2; N1 = N * 2; N.Release(); N2 = N * 2; Y = N; Print(N);
203  Y = (N1 - M1) | (N2 - M1); Print(Y);
204 #endif
205 
206  }
207 
208  {
209  Tracer et("Stage 13");
210  // test sum of squares of rows or columns
211  MultWithCarry mwc;
212  DiagonalMatrix DM; Matrix X;
213  // rectangular matrix
214  Matrix A(20, 15);
215  FillWithValues(mwc, A);
216  // sum of squares of rows
217  DM << A * A.t();
218  ColumnVector CV = A.sum_square_rows();
219  X = CV - DM.AsColumn(); Clean(X, 0.000000001); Print(X);
220  DM << A.t() * A;
221  RowVector RV = A.sum_square_columns();
222  X = RV - DM.AsRow(); Clean(X, 0.000000001); Print(X);
223  X = RV - A.t().sum_square_rows().t(); Clean(X, 0.000000001); Print(X);
224  X = CV - A.t().sum_square_columns().t(); Clean(X, 0.000000001); Print(X);
225  // UpperTriangularMatrix
226  A.ReSize(17,17); FillWithValues(mwc, A);
227  UpperTriangularMatrix UT; UT << A;
228  Matrix A1 = UT;
229  X = UT.sum_square_rows() - A1.sum_square_rows(); Print(X);
230  X = UT.sum_square_columns() - A1.sum_square_columns(); Print(X);
231  // LowerTriangularMatrix
232  LowerTriangularMatrix LT; LT << A;
233  A1 = LT;
234  X = LT.sum_square_rows() - A1.sum_square_rows(); Print(X);
235  X = LT.sum_square_columns() - A1.sum_square_columns(); Print(X);
236  // SymmetricMatrix
237  SymmetricMatrix SM; SM << A;
238  A1 = SM;
239  X = SM.sum_square_rows() - A1.sum_square_rows(); Print(X);
240  X = SM.sum_square_columns() - A1.sum_square_columns(); Print(X);
241  // DiagonalMatrix
242  DM << A;
243  A1 = DM;
244  X = DM.sum_square_rows() - A1.sum_square_rows(); Print(X);
245  X = DM.sum_square_columns() - A1.sum_square_columns(); Print(X);
246  // BandMatrix
247  BandMatrix BM(17, 3, 5); BM.Inject(A);
248  A1 = BM;
249  X = BM.sum_square_rows() - A1.sum_square_rows(); Print(X);
250  X = BM.sum_square_columns() - A1.sum_square_columns(); Print(X);
251  // SymmetricBandMatrix
252  SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
253  A1 = SBM;
254  X = SBM.sum_square_rows() - A1.sum_square_rows(); Print(X);
255  X = SBM.sum_square_columns() - A1.sum_square_columns(); Print(X);
256  // IdentityMatrix
257  IdentityMatrix IM(29);
258  X = IM.sum_square_rows() - 1; Print(X);
259  X = IM.sum_square_columns() - 1; Print(X);
260  // Matrix with zero rows
261  A1.ReSize(0,10);
262  X.ReSize(1,10); X = 0; X -= A1.sum_square_columns(); Print(X);
263  X.ReSize(0,1); X -= A1.sum_square_rows(); Print(X);
264  // Matrix with zero columns
265  A1.ReSize(10,0);
266  X.ReSize(10,1); X = 0; X -= A1.sum_square_rows(); Print(X);
267  X.ReSize(1,0); X -= A1.sum_square_columns(); Print(X);
268 
269  }
270 
271  {
272  Tracer et("Stage 14");
273  // test extend orthonormal
274  MultWithCarry mwc;
275  Matrix A(20,5); FillWithValues(mwc, A);
276  // Orthonormalise
277  UpperTriangularMatrix R;
278  Matrix A_old = A;
279  QRZ(A,R);
280  // Check decomposition
281  Matrix X = A * R - A_old; Clean(X, 0.000000001); Print(X);
282  // Check orthogonality
283  X = A.t() * A - IdentityMatrix(5);
284  Clean(X, 0.000000001); Print(X);
285  // Try orthonality extend
286  SquareMatrix A1(20);
287  A1.Columns(1,5) = A;
288  extend_orthonormal(A1,5);
289  // check columns unchanged
290  X = A - A1.Columns(1,5); Print(X);
291  // Check orthogonality
292  X = A1.t() * A1 - IdentityMatrix(20);
293  Clean(X, 0.000000001); Print(X);
294  X = A1 * A1.t() - IdentityMatrix(20);
295  Clean(X, 0.000000001); Print(X);
296  // Test with smaller number of columns
297  Matrix A2(20,15);
298  A2.Columns(1,5) = A;
299  extend_orthonormal(A2,5);
300  // check columns unchanged
301  X = A - A2.Columns(1,5); Print(X);
302  // Check orthogonality
303  X = A2.t() * A2 - IdentityMatrix(15);
304  Clean(X, 0.000000001); Print(X);
305  // check it works with no columns to start with
306  A2.ReSize(100,100);
307  extend_orthonormal(A2,0);
308  // Check orthogonality
309  X = A2.t() * A2 - IdentityMatrix(100);
310  Clean(X, 0.000000001); Print(X);
311  X = A2 * A2.t() - IdentityMatrix(100);
312  Clean(X, 0.000000001); Print(X);
313  }
314 
315  {
316  Tracer et("Stage 15");
317  // test extend orthonormal for SVD
318  MultWithCarry mwc;
319  Matrix A(15,40); FillWithValues(mwc, A);
320  Matrix U, V; DiagonalMatrix D;
321  SVD(A.t(),D,V,U);
322  Matrix VE(40,40); VE.columns(1,15) = V;
323  extend_orthonormal(VE,15);
324  Matrix DE(15,40); DE = 0; DE.columns(1,15) = D;
325  Matrix B = U * DE * VE.t();
326  B -= A;
327  Clean(B, 0.000000001); Print(B);
328  }
329 
330  {
331  Tracer et("Stage 16");
332  // test sum of rows or columns
333  MultWithCarry mwc;
334  ColumnVector CVX; Matrix X;
335  // rectangular matrix
336  Matrix A(20, 15);
337  FillWithValues(mwc, A);
338  // sum of rows
339  ColumnVector Ones(15); Ones = 1;
340  CVX = A * Ones;
341  ColumnVector CV = A.sum_rows();
342  X = CV - CVX; Clean(X, 0.000000001); Print(X);
343  Ones.resize(20); Ones = 1;
344  CVX << A.t() * Ones;
345  RowVector RV = A.sum_columns();
346  X = RV - CVX.AsRow(); Clean(X, 0.000000001); Print(X);
347  X = RV - A.t().sum_rows().t(); Clean(X, 0.000000001); Print(X);
348  X = CV - A.t().sum_columns().t(); Clean(X, 0.000000001); Print(X);
349  // UpperTriangularMatrix
350  A.ReSize(17,17); FillWithValues(mwc, A);
351  UpperTriangularMatrix UT; UT << A;
352  Matrix A1 = UT;
353  X = UT.sum_rows() - A1.sum_rows(); Print(X);
354  X = UT.sum_columns() - A1.sum_columns(); Print(X);
355  // LowerTriangularMatrix
356  LowerTriangularMatrix LT; LT << A;
357  A1 = LT;
358  X = LT.sum_rows() - A1.sum_rows(); Print(X);
359  X = LT.sum_columns() - A1.sum_columns(); Print(X);
360  // SymmetricMatrix
361  SymmetricMatrix SM; SM << A;
362  A1 = SM;
363  X = SM.sum_rows() - A1.sum_rows(); Print(X);
364  X = SM.sum_columns() - A1.sum_columns(); Print(X);
365  // DiagonalMatrix
366  DiagonalMatrix DM; DM << A;
367  A1 = DM;
368  X = DM.sum_rows() - A1.sum_rows(); Print(X);
369  X = DM.sum_columns() - A1.sum_columns(); Print(X);
370  // BandMatrix
371  BandMatrix BM(17, 3, 5); BM.Inject(A);
372  A1 = BM;
373  X = BM.sum_rows() - A1.sum_rows(); Print(X);
374  X = BM.sum_columns() - A1.sum_columns(); Print(X);
375  // SymmetricBandMatrix
376  SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
377  A1 = SBM;
378  X = SBM.sum_rows() - A1.sum_rows(); Print(X);
379  X = SBM.sum_columns() - A1.sum_columns(); Print(X);
380  // IdentityMatrix
381  IdentityMatrix IM(29);
382  X = IM.sum_rows() - 1; Print(X);
383  X = IM.sum_columns() - 1; Print(X);
384  // Matrix with zero rows
385  A1.ReSize(0,10);
386  X.ReSize(1,10); X = 0; X -= A1.sum_columns(); Print(X);
387  X.ReSize(0,1); X -= A1.sum_rows(); Print(X);
388  // Matrix with zero columns
389  A1.ReSize(10,0);
390  X.ReSize(10,1); X = 0; X -= A1.sum_rows(); Print(X);
391  X.ReSize(1,0); X -= A1.sum_columns(); Print(X);
392 
393  }
394 
395 
396 
397 // cout << "\nEnd of twelfth test\n";
398 }
399 
400 
GeneralMatrix::Release
void Release()
Definition: newmat.h:514
BaseMatrix::Columns
GetSubMatrix Columns(int f, int l) const
Definition: newmat.h:2153
BaseMatrix::columns
GetSubMatrix columns(int, int) const
Definition: submat.cpp:77
Matrix::ReSize
virtual void ReSize(int m, int n)
Definition: newmat.h:662
GeneralMatrix::Inject
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
BaseMatrix::sum_rows
ReturnMatrix sum_rows() const
Definition: newmat8.cpp:781
DiagonalMatrix::ReSize
void ReSize(int m)
Definition: newmat.h:935
BaseMatrix::sum_columns
ReturnMatrix sum_columns() const
Definition: newmat8.cpp:805
UpperTriangularMatrix
Upper triangular matrix.
Definition: newmat.h:799
SquareMatrix
Square matrix.
Definition: newmat.h:679
Clean
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
BaseMatrix::Column
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
BaseMatrix::AsRow
RowedMatrix AsRow() const
Definition: newmat.h:2141
BandMatrix
Band matrix.
Definition: newmat.h:1096
BaseMatrix::t
TransposedMatrix t() const
Definition: newmat6.cpp:320
extend_orthonormal
void extend_orthonormal(Matrix &A, int n)
Definition: hholder.cpp:341
SVD
void SVD(const Matrix &, DiagonalMatrix &, Matrix &, Matrix &, bool=true, bool=true)
Definition: svd.cpp:30
BaseMatrix::sum_square_rows
ReturnMatrix sum_square_rows() const
Definition: newmat8.cpp:736
Tracer::PrintTrace
static void PrintTrace()
Definition: myexcept.cpp:109
Matrix
The usual rectangular matrix.
Definition: newmat.h:625
BaseMatrix::i
InvertedMatrix i() const
Definition: newmat6.cpp:329
BaseMatrix::SymSubMatrix
GetSubMatrix SymSubMatrix(int f, int l) const
Definition: newmat.h:2148
nricMatrix
Rectangular matrix for use with Numerical Recipes in C.
Definition: newmat.h:711
QRZ
void QRZ(Matrix &X, UpperTriangularMatrix &U)
Definition: hholder.cpp:117
BaseMatrix::AsMatrix
MatedMatrix AsMatrix(int m, int n) const
Definition: newmat.h:2144
DiagonalMatrix
Diagonal matrix.
Definition: newmat.h:896
LowerTriangularMatrix
Lower triangular matrix.
Definition: newmat.h:848
BaseMatrix::Row
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
BaseMatrix::SubMatrix
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const
Definition: newmat.h:2146
SymmetricBandMatrix
Symmetric band matrix.
Definition: newmat.h:1245
BaseMatrix::sum_square_columns
ReturnMatrix sum_square_columns() const
Definition: newmat8.cpp:760
RowVector
Row vector.
Definition: newmat.h:953
BaseMatrix::AsColumn
ColedMatrix AsColumn() const
Definition: newmat.h:2142
Print
void Print(const Matrix &X)
Definition: tmt.cpp:42
ColumnVector
Column vector.
Definition: newmat.h:1008
trymatc
void trymatc()
Definition: tmtc.cpp:23
IdentityMatrix
Identity matrix.
Definition: newmat.h:1350
FillWithValues
void FillWithValues(MultWithCarry &MWC, Matrix &M)
Definition: tmt.cpp:190
SymmetricMatrix
Symmetric matrix.
Definition: newmat.h:753

kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:17