tmt7.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 #include "include.h"
11 
12 #include "newmat.h"
13 
14 #include "tmt.h"
15 
16 #ifdef use_namespace
17 using namespace NEWMAT;
18 #endif
19 
20 // C matrix mulitply - for testing RealStarStar
21 
22 void c_matrix_multiply(int p, int q, int r,
23  const Real** a, const Real** b, Real** c)
24 {
25  for (int i = 0; i < p; ++i) for (int k = 0; k < r; ++k)
26  {
27  Real sum = 0.0;
28  for (int j = 0; j < q; ++j) sum += a[i][j] * b[j][k];
29  c[i][k] = sum;
30  }
31 }
32 
33 
34 
35 void trymat7()
36 {
37 // cout << "\nSeventh test of Matrix package\n";
38  Tracer et("Seventh test of Matrix package");
40 
41  int i,j;
42 
43 
44  DiagonalMatrix D(6);
46  for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*j-50; D(i,i)=i*i+i-10; }
47  LowerTriangularMatrix L=(U*3.0).t();
48  SymmetricMatrix S(6);
49  for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
50  Matrix MD=D; Matrix ML=L; Matrix MU=U;
51  Matrix MS=S;
52  Matrix M(6,6);
53  for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;
54  {
55  Tracer et1("Stage 1");
56  Print(Matrix((S-M)-(MS-M)));
57  Print(Matrix((-M-S)+(MS+M)));
58  Print(Matrix((U-M)-(MU-M)));
59  }
60  {
61  Tracer et1("Stage 2");
62  Print(Matrix((L-M)+(M-ML)));
63  Print(Matrix((D-M)+(M-MD)));
64  Print(Matrix((D-S)+(MS-MD)));
65  Print(Matrix((D-L)+(ML-MD)));
66  }
67 
68  { M=MU.t(); }
69  LowerTriangularMatrix LY=D.i()*U.t();
70  {
71  Tracer et1("Stage 3");
72  MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
73  L=U.t();
74  LY=D.i()*L; MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
75  }
76  {
77  Tracer et1("Stage 4");
78  UpperTriangularMatrix UT(11);
79  int i, j;
80  for (i=1;i<=11;i++) for (j=i;j<=11;j++) UT(i,j)=i*i+j*3;
81  GenericMatrix GM; Matrix X;
82  UpperBandMatrix UB(11,3); UB.Inject(UT); UT = UB;
83  UpperBandMatrix UB2 = UB / 8;
84  GM = UB2-UT/8; X = GM; Print(X);
85  SymmetricBandMatrix SB(11,4); SB << (UB + UB.t());
86  X = SB - UT - UT.t(); Print(X);
87  BandMatrix B = UB + UB.t()*2;
88  DiagonalMatrix D; D << B;
89  X.ReSize(1,1); X(1,1) = Trace(B)-Sum(D); Print(X);
90  X = SB + 5; Matrix X1=X; X = SP(UB,X); Matrix X2 =UB;
91  X1 = (X1.AsDiagonal() * X2.AsDiagonal()).AsRow()-X.AsColumn().t();
92  Print(X1);
93  X1=SB.t(); X2 = B.t(); X = SB.i() * B - X1.i() * X2.t();
94  Clean(X,0.00000001); Print(X);
95  X = SB.i(); X = X * B - X1.i() * X2.t();
96  Clean(X,0.00000001); Print(X);
97  D = 1; X = SB.i() * SB - D; Clean(X,0.00000001); Print(X);
98  ColumnVector CV(11);
99  CV << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
100  D << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
101  X = CV.AsDiagonal(); X = X-D; Print(X);
102  SymmetricBandMatrix SB1(11,7); SB1 = 5;
103  SymmetricBandMatrix SB2 = SB1 + D;
104  X.ReSize(11,11); X=0;
105  for (i=1;i<=11;i++) for (j=1;j<=11;j++)
106  {
107  if (abs(i-j)<=7) X(i,j)=5;
108  if (i==j) X(i,j)+=CV(i);
109  }
110  SymmetricMatrix SM; SM.ReSize(11);
111  SM=SB; SB = SB+SB2; X1 = SM+X-SB; Print(X1);
112  SB2=0; X2=SB2; X1=SB; Print(X2);
113  for (i=1;i<=11;i++) SB2.Column(i)<<SB.Column(i);
114  X1=X1-SB2; Print(X1);
115  X = SB; SB2.ReSize(11,4); SB2 = SB*5; SB2 = SB + SB2;
116  X1 = X*6 - SB2; Print(X1);
117  X1 = SP(SB,SB2/3); X1=X1-SP(X,X*2); Print(X1);
118  X1 = SP(SB2/6,X*2); X1=X1-SP(X*2,X); Print(X1);
119  }
120 
121  {
122  // test the simple integer array class
123  Tracer et("Stage 5");
124  ColumnVector Test(10); Test = 0.0;
125  int i;
126  SimpleIntArray A(100);
127  for (i = 0; i < 100; i++) A[i] = i*i+1;
128  SimpleIntArray B(100), C(50), D;
129  B = A; A.ReSize(50, true); C = A; A.ReSize(150, true); D = A;
130  for (i = 0; i < 100; i++) if (B[i] != i*i+1) Test(1)=1;
131  for (i = 0; i < 50; i++) if (C[i] != i*i+1) Test(2)=1;
132  for (i = 0; i < 50; i++) if (D[i] != i*i+1) Test(3)=1;
133  for (i = 50; i < 150; i++) if (D[i] != 0) Test(3)=1;
134  A.resize(75); A = A.size();
135  for (i = 0; i < 75; i++) if (A[i] != 75) Test(4)=1;
136  A.resize(25); A = A.size();
137  for (i = 0; i < 25; i++) if (A[i] != 25) Test(5)=1;
138  A.ReSize(25); A = 23;
139  for (i = 0; i < 25; i++) if (A[i] != 23) Test(6)=1;
140  A.ReSize(0); A.ReSize(15); A = A.Size();
141  for (i = 0; i < 15; i++) if (A[i] != 15) Test(7)=1;
142  const SimpleIntArray E = B;
143  for (i = 0; i < 100; i++) if (E[i] != i*i+1) Test(8)=1;
144  SimpleIntArray F; F.resize_keep(5);
145  for (i = 0; i < 5; i++) if (F[i] != 0) Test(9)=1;
146  Print(Test);
147  }
148 
149  {
150  // testing RealStarStar
151  Tracer et("Stage 6");
152  MultWithCarry MWC;
153 
154  Matrix A(10, 12), B(12, 15), C(10, 15);
155  FillWithValues(MWC, A); FillWithValues(MWC, B);
156  ConstRealStarStar a(A);
157  ConstRealStarStar b(B);
158  RealStarStar c(C);
159  c_matrix_multiply(10,12,15,a,b,c);
160  Matrix X = C - A * B; Clean(X,0.00000001); Print(X);
161  A.ReSize(11, 10); B.ReSize(10,8); C.ReSize(11,8);
162  FillWithValues(MWC, A); FillWithValues(MWC, B);
163  C = -1;
164  c_matrix_multiply(11,10,8,
166  X = C - A * B; Clean(X,0.00000001); Print(X);
167  }
168 
169  {
170  // testing resize_keep
171  Tracer et("Stage 7");
172  Matrix X, Y;
173  MultWithCarry MWC;
174 
175  X.resize(20,35); FillWithValues(MWC, X);
176  Matrix M(20,35); M = X;
177  X = M.submatrix(1,15,1,25);
178  M.resize_keep(15,25); Y = X - M; Print(Y);
179  M.resize_keep(15,25); Y = X - M; Print(Y);
180  Y.resize(29,27); Y = 0; Y.submatrix(1,15,1,25) = X;
181  M.resize_keep(29,27); Y -= M; Print(Y);
182  M.resize_keep(0,5); M.resize_keep(10,10); Print(M);
183  M.resize_keep(15,0); M.resize_keep(10,10); Print(M);
184 
185  X.resize(20,35); FillWithValues(MWC, X);
186  M = X;
187  M.resize_keep(38,17);
188  Y.resize(38,17); Y = 0;
189  Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17);
190  Y -= M; Print(Y);
191 
192  X.resize(40,12); FillWithValues(MWC, X);
193  M = X;
194  M.resize_keep(38,17);
195  Y.resize(38,17); Y = 0;
196  Y.submatrix(1,38,1,12) = X.submatrix(1,38,1,12);
197  Y -= M; Print(Y);
198 
199 #ifndef DONT_DO_NRIC
200 
201  X.resize(20,35); FillWithValues(MWC, X);
202  nricMatrix nM(20,35); nM = X;
203  X = nM.submatrix(1,15,1,25);
204  nM.resize_keep(15,25); Y = X - nM; Print(Y);
205  nM.resize_keep(15,25); Y = X - nM; Print(Y);
206  Y.resize(29,27); Y = 0; Y.submatrix(1,15,1,25) = X;
207  nM.resize_keep(29,27); Y -= nM; Print(Y);
208  nM.resize_keep(0,5); nM.resize_keep(10,10); Print(nM);
209  nM.resize_keep(15,0); nM.resize_keep(10,10); Print(nM);
210 
211  X.resize(20,35); FillWithValues(MWC, X);
212  nM = X;
213  nM.resize_keep(38,17);
214  Y.resize(38,17); Y = 0;
215  Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17);
216  Y -= nM; Print(Y);
217 
218  X.resize(40,12); FillWithValues(MWC, X);
219  nM = X;
220  nM.resize_keep(38,17);
221  Y.resize(38,17); Y = 0;
222  Y.submatrix(1,38,1,12) = X.submatrix(1,38,1,12);
223  Y -= nM; Print(Y);
224 
225 #endif
226 
227  X.resize(20,20); FillWithValues(MWC, X);
228  SquareMatrix SQM(20); SQM << X;
229  X = SQM.sym_submatrix(1,13);
230  SQM.resize_keep(13); Y = X - SQM; Print(Y);
231  SQM.resize_keep(13); Y = X - SQM; Print(Y);
232  Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
233  SQM.resize_keep(23,23); Y -= SQM; Print(Y);
234  SQM.resize_keep(0); SQM.resize_keep(50); Print(SQM);
235 
236  X.resize(20,20); FillWithValues(MWC, X);
237  SymmetricMatrix SM(20); SM << X;
238  X = SM.sym_submatrix(1,13);
239  SM.resize_keep(13); Y = X - SM; Print(Y);
240  SM.resize_keep(13); Y = X - SM; Print(Y);
241  Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
242  SM.resize_keep(23); Y -= SM; Print(Y);
243  SM.resize_keep(0); SM.resize_keep(50); Print(SM);
244 
245  X.resize(20,20); FillWithValues(MWC, X);
246  LowerTriangularMatrix LT(20); LT << X;
247  X = LT.sym_submatrix(1,13);
248  LT.resize_keep(13); Y = X - LT; Print(Y);
249  LT.resize_keep(13); Y = X - LT; Print(Y);
250  Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
251  LT.resize_keep(23); Y -= LT; Print(Y);
252  LT.resize_keep(0); LT.resize_keep(50); Print(LT);
253 
254  X.resize(20,20); FillWithValues(MWC, X);
255  UpperTriangularMatrix UT(20); UT << X;
256  X = UT.sym_submatrix(1,13);
257  UT.resize_keep(13); Y = X - UT; Print(Y);
258  UT.resize_keep(13); Y = X - UT; Print(Y);
259  Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
260  UT.resize_keep(23); Y -= UT; Print(Y);
261  UT.resize_keep(0); UT.resize_keep(50); Print(UT);
262 
263  X.resize(20,20); FillWithValues(MWC, X);
264  DiagonalMatrix DM(20); DM << X;
265  X = DM.sym_submatrix(1,13);
266  DM.resize_keep(13); Y = X - DM; Print(Y);
267  DM.resize_keep(13); Y = X - DM; Print(Y);
268  Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
269  DM.resize_keep(23); Y -= DM; Print(Y);
270  DM.resize_keep(0); DM.resize_keep(50); Print(DM);
271 
272  X.resize(1,20); FillWithValues(MWC, X);
273  RowVector RV(20); RV << X;
274  X = RV.columns(1,13);
275  RV.resize_keep(13); Y = X - RV; Print(Y);
276  RV.resize_keep(13); Y = X - RV; Print(Y);
277  Y.resize(1,23); Y = 0; Y.columns(1,13) = X;
278  RV.resize_keep(1,23); Y -= RV; Print(Y);
279  RV.resize_keep(0); RV.resize_keep(50); Print(RV);
280 
281  X.resize(20,1); FillWithValues(MWC, X);
282  ColumnVector CV(20); CV << X;
283  X = CV.rows(1,13);
284  CV.resize_keep(13); Y = X - CV; Print(Y);
285  CV.resize_keep(13); Y = X - CV; Print(Y);
286  Y.resize(23,1); Y = 0; Y.rows(1,13) = X;
287  CV.resize_keep(23,1); Y -= CV; Print(Y);
288  CV.resize_keep(0); CV.resize_keep(50); Print(CV);
289 
290 
291  }
292 
293 
294 // cout << "\nEnd of seventh test\n";
295 }
296 
297 
298 
Real Sum(const BaseMatrix &B)
Definition: newmat.h:2107
void resize_keep(int i)
change length, keep data
Definition: newmat.h:1882
void ReSize(int m)
Definition: newmat.h:786
GetSubMatrix rows(int, int) const
Definition: submat.cpp:58
void resize_keep(int)
Definition: newmat4.cpp:524
SPMatrix SP(const BaseMatrix &, const BaseMatrix &)
Definition: newmat6.cpp:278
GetSubMatrix columns(int, int) const
Definition: submat.cpp:77
virtual void ReSize(int m, int n)
Definition: newmat.h:662
void resize_keep(int)
Definition: newmat4.cpp:560
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
double Real
Definition: include.h:307
void c_matrix_multiply(int p, int q, int r, const Real **a, const Real **b, Real **c)
Definition: tmt7.cpp:22
virtual void ReSize(int m, int n, int b)
Definition: newmat.h:1145
Upper triangular band matrix.
Definition: newmat.h:1167
Upper triangular matrix.
Definition: newmat.h:799
void resize_keep(int)
Definition: newmat4.cpp:506
Square matrix.
Definition: newmat.h:679
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
void resize_keep(int)
Definition: newmat4.cpp:461
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
GetSubMatrix sym_submatrix(int, int) const
Definition: submat.cpp:39
Real Trace(const BaseMatrix &B)
Definition: newmat.h:2100
void resize_keep(int)
Definition: newmat4.cpp:542
int Size() const
return the size of the array
Definition: newmat.h:1869
void resize_keep(int)
Definition: newmat4.cpp:488
void ReSize(int m, int b)
Definition: newmat.h:1289
Let matrix simulate a C type two dimensional array.
Definition: newmat.h:1892
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
Rectangular matrix for use with Numerical Recipes in C.
Definition: newmat.h:711
FloatVector FloatVector * a
void resize(int i, bool keep=false)
change length, keep data if keep = true
Definition: newmat4.cpp:1228
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
int size() const
return the size of the array
Definition: newmat.h:1871
void trymat7()
Definition: tmt7.cpp:35
Symmetric band matrix.
Definition: newmat.h:1245
virtual void resize(int, int)
Definition: newmat4.cpp:289
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
DiagedMatrix AsDiagonal() const
Definition: newmat.h:2143
Row vector.
Definition: newmat.h:953
ColedMatrix AsColumn() const
Definition: newmat.h:2142
void resize_keep(int)
Definition: newmat4.cpp:586
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
virtual void resize_keep(int, int)
Definition: newmat4.cpp:432
GetSubMatrix submatrix(int, int, int, int) const
Definition: submat.cpp:27
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
void resize_keep(int m, int n)
Definition: newmat.h:734
Let matrix simulate a C type const two dimensional array.
Definition: newmat.h:1902
void FillWithValues(MultWithCarry &MWC, Matrix &M)
Definition: tmt.cpp:190
Symmetric matrix.
Definition: newmat.h:753
void ReSize(int i, bool keep=false)
change length, keep data if keep = true
Definition: newmat.h:1880


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