1 /// \ingroup newmat
2 ///@{
3
4 /// \file tmt5.cpp
5 /// Part of matrix library test program.
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
21 // **************************** test program ******************************
22
23
24
Returner0(const GenericMatrix & GM)25 ReturnMatrix Returner0(const GenericMatrix& GM)
26 { Matrix M = GM; M.Release(); return M; }
27
Returner1(const GenericMatrix & GM)28 ReturnMatrix Returner1(const GenericMatrix& GM)
29 { Matrix M = GM+1; M.Release(); return M; }
30
Returner2(const GenericMatrix & GM)31 ReturnMatrix Returner2(const GenericMatrix& GM)
32 { UpperBandMatrix M = GM*2; M.Release(); return M; }
33
Returner3(const GenericMatrix & GM)34 ReturnMatrix Returner3(const GenericMatrix& GM)
35 { LowerBandMatrix M = GM*3; M.Release(); return M; }
36
Returner4(const GenericMatrix & GM)37 ReturnMatrix Returner4(const GenericMatrix& GM)
38 { SymmetricMatrix M = GM+4; M.Release(); return M; }
39
Returner5(const GenericMatrix & GM)40 ReturnMatrix Returner5(const GenericMatrix& GM)
41 { SymmetricBandMatrix M = GM*5; M.Release(); return M; }
42
Returner6(const GenericMatrix & GM)43 ReturnMatrix Returner6(const GenericMatrix& GM)
44 { BandMatrix M = GM*6; M.Release(); return M; }
45
Returner7(const GenericMatrix & GM)46 ReturnMatrix Returner7(const GenericMatrix& GM)
47 { DiagonalMatrix M = GM*7; M.Release(); return M; }
48
trymat5()49 void trymat5()
50 {
51 // cout << "\nFifth test of Matrix package\n";
52 Tracer et("Fifth test of Matrix package");
53 Tracer::PrintTrace();
54
55 int i,j;
56
57 Matrix A(5,6);
58 for (i=1;i<=5;i++) for (j=1;j<=6;j++) A(i,j)=1+i*j+i*i+j*j;
59 ColumnVector CV(6);
60 for (i=1;i<=6;i++) CV(i)=i*i+3;
61 ColumnVector CV2(5); for (i=1;i<=5;i++) CV2(i)=1.0;
62 ColumnVector CV1=CV;
63
64 {
65 CV=A*CV;
66 RowVector RV=CV.t(); // RowVector RV; RV=CV.t();
67 RV=RV-1.0;
68 CV=(RV*A).t()+A.t()*CV2; CV1=(A.t()*A)*CV1 - CV;
69 Print(CV1);
70 }
71
72 CV1.ReSize(6);
73 CV2.ReSize(6);
74 CV.ReSize(6);
75 for (i=1;i<=6;i++) { CV1(i)=i*3+1; CV2(i)=10-i; CV(i)=11+i*2; }
76 ColumnVector CX=CV2-CV; { CX=CX+CV1; Print(CX); }
77 Print(ColumnVector(CV1+CV2-CV));
78 RowVector RV=CV.t(); RowVector RV1=CV1.t();
79 RowVector R=RV-RV1; Print(RowVector(R-CV2.t()));
80
81 // test loading of list
82
83 RV.ReSize(10);
84 for (i=1;i<=10;i++) RV(i) = i*i;
85 RV1.ReSize(10);
86 RV1 << 1 << 4 << 9 << 16 << 25 << 36 << 49 << 64 << 81 << 100; // << 121;
87 Print(RowVector(RV-RV1));
88
89 et.ReName("Fifth test of Matrix package - almost at end");
90
91 Matrix X(2,3);
92 X << 11 << 12 << 13
93 << 21 << 22 << 23;
94
95 Matrix Y = X.t(); // check simple transpose
96
97 X(1,1) -= 11; X(1,2) -= 12; X(1,3) -= 13;
98 X(2,1) -= 21; X(2,2) -= 22; X(2,3) -= 23;
99 Print(X);
100
101 Y(1,1) -= 11; Y(2,1) -= 12; Y(3,1) -= 13;
102 Y(1,2) -= 21; Y(2,2) -= 22; Y(3,2) -= 23;
103 Print(Y);
104
105 et.ReName("Fifth test of Matrix package - at end");
106
107 RV = Returner1(RV)-1; Print(RowVector(RV-RV1));
108 CV1 = Returner1(RV.t())-1; Print(ColumnVector(RV.t()-CV1));
109 #ifndef DONT_DO_NRIC
110 nricMatrix AA = A;
111 X = Returner1(AA)-A-1; Print(X);
112 #endif
113 UpperTriangularMatrix UT(31);
114 for (i=1; i<=31; i++) for (j=i; j<=31; j++) UT(i,j) = i+j+(i-j)*(i-2*j);
115 UpperBandMatrix UB(31,5); UB.Inject(UT);
116 LowerTriangularMatrix LT = UT.t();
117 LowerBandMatrix LB(31,5); LB.Inject(LT);
118 A = Returner0(UB)-LB.t(); Print(A);
119 A = Returner2(UB).t()-LB*2; Print(A);
120 A = Returner3(LB).t()-UB*3; Print(A);
121 SymmetricMatrix SM; SM << (UT+LT);
122 A = Returner4(SM)-UT-LT-4; Print(A);
123 SymmetricBandMatrix SB(31,5); SB.Inject(SM);
124 A = Returner5(SB)/5-UB-LB; Print(A);
125 BandMatrix B = UB+LB*LB; A = LB;
126 A = Returner6(B)/6 - UB - A*A; Print(A);
127 DiagonalMatrix D; D << UT;
128 D << (Returner7(D)/7 - UT); Print(D);
129
130 // cout << "\nEnd of fifth test\n";
131 }
132
133
134
135 ///@}
136