tmt6.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 #define WANT_MATH
10 
11 
12 #include "include.h"
13 
14 #include "newmatap.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 // slow sort program
27 
28 static void SimpleSortDescending(Real* first, const int length)
29 {
30  int i = length;
31  while (--i)
32  {
33  Real x = *first; Real* f = first; Real* g = f;
34  int j = i;
35  while (j--) if (x < *(++f)) { g = f; x = *g; }
36  *g = *first; *first++ = x;
37  }
38 }
39 
40 static void TestSort(int n)
41 {
42  // make some data
43  RowVector X(n);
44  int i;
45  for (i = 1; i <= n; i++)
46  X(i) = sin((Real)i) + 0.3 * cos(i/5.0) - 0.6 * sin(i/7.0) + 0.2 * sin(2.0 * i);
47  RowVector X_Sorted = X; SimpleSortDescending(X_Sorted.Store(), n);
48  RowVector A = X + X.Reverse(); SimpleSortDescending(A.Store(), n);
49 
50  // test descending sort
51 
52  RowVector Y = X; SortDescending(Y); Y -= X_Sorted; Print(Y);
53  Y = X_Sorted; SortDescending(Y); Y -= X_Sorted; Print(Y);
54  Y = X_Sorted.Reverse(); SortDescending(Y); Y -= X_Sorted; Print(Y);
55  Y = X + X.Reverse(); SortDescending(Y); Y -= A; Print(Y);
56 
57  // test ascending sort
58 
59  Y = X; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
60  Y = X_Sorted; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
61  Y = X_Sorted.Reverse(); SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
62  Y = X + X.Reverse(); SortAscending(Y); Y -= A.Reverse(); Print(Y);
63 }
64 
65 
66 void trymat6()
67 {
68  Tracer et("Sixth test of Matrix package");
70 
71  int i,j;
72 
73 
74  DiagonalMatrix D(6);
76  for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; }
77  LowerTriangularMatrix L=(U*3.0).t();
78  SymmetricMatrix S(6);
79  for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
80  Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S;
81  Matrix M(6,6);
82  for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;
83  {
84  Tracer et1("Stage 1");
85  Print(Matrix(MS+(-MS)));
86  Print(Matrix((S+M)-(MS+M)));
87  Print(Matrix((M+U)-(M+MU)));
88  Print(Matrix((M+L)-(M+ML)));
89  }
90  {
91  Tracer et1("Stage 2");
92  Print(Matrix((M+D)-(M+MD)));
93  Print(Matrix((U+D)-(MU+MD)));
94  Print(Matrix((D+L)-(ML+MD)));
95  Print(Matrix((-U+D)+MU-MD));
96  Print(Matrix((-L+D)+ML-MD));
97  }
98  {
99  Tracer et1("Stage 3 - concatenate");
100  RowVector A(5);
101  A << 1 << 2 << 3 << 4 << 5;
102  RowVector B(5);
103  B << 3 << 1 << 4 << 1 << 5;
104  Matrix C(3,5);
105  C << 2 << 3 << 5 << 7 << 11
106  << 13 << 17 << 19 << 23 << 29
107  << 31 << 37 << 41 << 43 << 47;
108  Matrix X1 = A & B & C;
109  Matrix X2 = (A.t() | B.t() | C.t()).t();
110  Matrix X3(5,5);
111  X3.Row(1)=A; X3.Row(2)=B; X3.Rows(3,5)=C;
112  Print(Matrix(X1-X2));
113  Print(Matrix(X1-X3));
114  LowerTriangularMatrix LT1; LT1 << (A & B & C);
115  UpperTriangularMatrix UT1; UT1 << (A.t() | B.t() | C.t());
116  Print(LowerTriangularMatrix(LT1-UT1.t()));
117  DiagonalMatrix D1; D1 << (A.t() | B.t() | C.t());
118  ColumnVector At = A.t();
119  ColumnVector Bt = B.t();
120  Matrix Ct = C.t();
121  LowerTriangularMatrix LT2; LT2 << (At | Bt | Ct);
122  UpperTriangularMatrix UT2; UT2 << (At.t() & Bt.t() & Ct.t());
123  Matrix ABt = At | Bt;
124  DiagonalMatrix D2; D2 << (ABt | Ct);
125  Print(LowerTriangularMatrix(LT2-UT2.t()));
126  Print(DiagonalMatrix(D1-D2));
127  Print(Matrix(LT1+UT2-D2-X1));
128  Matrix M1 = LT1 | UT2; Matrix M2 = UT1 & LT2;
129  Print(Matrix(M1-M2.t()));
130  M1 = UT2 | LT1; M2 = LT2 & UT1;
131  Print(Matrix(M1-M2.t()));
132  M1 = (LT1 | UT2) & (UT2 | LT1);
133  M2 = (UT1 & LT2) | (LT2 & UT1);
134  Print(Matrix(M1-M2.t()));
135  SymmetricMatrix SM1; SM1 << (M1 + M1.t());
136  SymmetricMatrix SM2; SM2 << ((SM1 | M1) & (M1.t() | SM1));
137  Matrix M3(20,20);
138  M3.SubMatrix(1,10,1,10) = SM1;
139  M3.SubMatrix(1,10,11,20) = M1;
140  M3.SubMatrix(11,20,1,10) = M2;
141  M3.SubMatrix(11,20,11,20) = SM1;
142  Print(Matrix(M3-SM2));
143 
144  SymmetricMatrix SM(15); SM = 0; SM.SymSubMatrix(1,10) = SM1;
145  M3.ReSize(15,15); M3 = 0; M3.SubMatrix(1,10,1,10) = SM1;
146  M3 -= SM; Print(M3);
147  SM = 0; SM.SymSubMatrix(6,15) = SM1;
148  M3.ReSize(15,15); M3 = 0; M3.SubMatrix(6,15,6,15) = SM1;
149  M3 = M3.t() - SM; Print(M3);
150  }
151  {
152  Tracer et1("Stage 4 - sort");
153  TestSort(1); TestSort(2); TestSort(3); TestSort(4);
154  TestSort(15); TestSort(16); TestSort(17); TestSort(18);
155  TestSort(99); TestSort(100); TestSort(101);
156  }
157 
158 
159 // cout << "\nEnd of sixth test\n";
160 }
161 
162 
void trymat6()
Definition: tmt6.cpp:66
virtual void ReSize(int m, int n)
Definition: newmat.h:662
double Real
Definition: include.h:307
void SortAscending(GeneralMatrix &gm)
Definition: newmatap.h:137
ReversedMatrix Reverse() const
Definition: newmat.h:2140
static void TestSort(int n)
Definition: tmt6.cpp:40
Upper triangular matrix.
Definition: newmat.h:799
TransposedMatrix t() const
Definition: newmat6.cpp:320
Real * Store() const
Definition: newmat.h:497
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
void SortDescending(GeneralMatrix &gm)
Definition: newmatap.h:139
GetSubMatrix SymSubMatrix(int f, int l) const
Definition: newmat.h:2148
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const
Definition: newmat.h:2146
Row vector.
Definition: newmat.h:953
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
static void SimpleSortDescending(Real *first, const int length)
Definition: tmt6.cpp:28
GetSubMatrix Rows(int f, int l) const
Definition: newmat.h:2151
Symmetric matrix.
Definition: newmat.h:753


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