$search
00001 00002 00003 00006 00007 00008 //#define WANT_STREAM 00009 00010 00011 #include "include.h" 00012 00013 #include "newmat.h" 00014 00015 #include "tmt.h" 00016 00017 #ifdef use_namespace 00018 using namespace NEWMAT; 00019 #endif 00020 00021 00022 /**************************** test program ******************************/ 00023 00024 00025 void trymat4() 00026 { 00027 // cout << "\nFourth test of Matrix package\n"; 00028 Tracer et("Fourth test of Matrix package"); 00029 Tracer::PrintTrace(); 00030 00031 00032 { 00033 Tracer et1("Stage 1"); 00034 int i, j; 00035 Matrix M(10,10); 00036 UpperTriangularMatrix U(10); 00037 for (i=1;i<=10;i++) for (j=1;j<=10;j++) M(i,j) = 100*i+j; 00038 U << -M; 00039 Matrix X1 = M.Rows(2,4); 00040 Matrix Y1 = U.t().Rows(2,4); 00041 Matrix X = U; { Print(Matrix(X.Columns(2,4).t()-Y1)); } 00042 RowVector RV = M.Row(5); 00043 { 00044 X.ReSize(3,10); 00045 X.Row(1) << M.Row(2); X.Row(2) << M.Row(3); X.Row(3) << M.Row(4); 00046 Print(Matrix(X-X1)); 00047 } 00048 { 00049 UpperTriangularMatrix V = U.SymSubMatrix(3,5); 00050 Matrix MV = U.SubMatrix(3,5,3,5); { Print(Matrix(MV-V)); } 00051 Matrix X2 = M.t().Columns(2,4); { Print(Matrix(X2-X1.t())); } 00052 Matrix Y2 = U.Columns(2,4); { Print(Matrix(Y2-Y1.t())); } 00053 ColumnVector CV = M.t().Column(5); { Print(ColumnVector(CV-RV.t())); } 00054 X.ReSize(10,3); M = M.t(); 00055 X.Column(1) << M.Column(2); X.Column(2) << M.Column(3); 00056 X.Column(3) << M.Column(4); 00057 Print(Matrix(X-X2)); 00058 } 00059 } 00060 00061 { 00062 Tracer et1("Stage 2"); 00063 int i, j; 00064 Matrix M; Matrix X; M.ReSize(5,8); 00065 for (i=1;i<=5;i++) for (j=1;j<=8;j++) M(i,j) = 100*i+j; 00066 { 00067 X = M.Columns(5,8); M.Columns(5,8) << M.Columns(1,4); 00068 M.Columns(1,4) << X; 00069 X = M.Columns(3,4); M.Columns(3,4) << M.Columns(1,2); 00070 M.Columns(1,2) << X; 00071 X = M.Columns(7,8); M.Columns(7,8) << M.Columns(5,6); 00072 M.Columns(5,6) << X; 00073 } 00074 { 00075 X = M.Column(2); M.Column(2) = M.Column(1); M.Column(1) = X; 00076 X = M.Column(4); M.Column(4) = M.Column(3); M.Column(3) = X; 00077 X = M.Column(6); M.Column(6) = M.Column(5); M.Column(5) = X; 00078 X = M.Column(8); M.Column(8) = M.Column(7); M.Column(7) = X; 00079 X.ReSize(5,8); 00080 } 00081 for (i=1;i<=5;i++) for (j=1;j<=8;j++) X(i,9-j) = 100*i+j; 00082 Print(Matrix(X-M)); 00083 } 00084 { 00085 Tracer et1("Stage 3"); 00086 // try submatrices of zero dimension 00087 int i, j; 00088 Matrix A(4,5); Matrix B, C; 00089 for (i=1; i<=4; i++) for (j=1; j<=5; j++) 00090 A(i,j) = 100+i*10+j; 00091 B = A + 100; 00092 C = A | B.Columns(4,3); Print(Matrix(A - C)); 00093 C = A | B.Columns(1,0); Print(Matrix(A - C)); 00094 C = A | B.Columns(6,5); Print(Matrix(A - C)); 00095 C = A & B.Rows(2,1); Print(Matrix(A - C)); 00096 } 00097 { 00098 Tracer et1("Stage 4"); 00099 BandMatrix BM(5,3,2); 00100 BM(1,1) = 1; BM(1,2) = 2; BM(1,3) = 3; 00101 BM(2,1) = 4; BM(2,2) = 5; BM(2,3) = 6; BM(2,4) = 7; 00102 BM(3,1) = 8; BM(3,2) = 9; BM(3,3) =10; BM(3,4) =11; BM(3,5) =12; 00103 BM(4,1) =13; BM(4,2) =14; BM(4,3) =15; BM(4,4) =16; BM(4,5) =17; 00104 BM(5,2) =18; BM(5,3) =19; BM(5,4) =20; BM(5,5) =21; 00105 SymmetricBandMatrix SM(5,3); 00106 SM.Inject(BandMatrix(BM + BM.t())); 00107 Matrix A = BM + 1; 00108 Matrix M = A + A.t() - 2; 00109 Matrix C = A.i() * BM; 00110 C = A * C - BM; Clean(C, 0.000000001); Print(C); 00111 C = A.i() * SM; 00112 C = A * C - M; Clean(C, 0.000000001); Print(C); 00113 00114 // check row-wise load 00115 BandMatrix BM1(5,3,2); 00116 BM1.Row(1) << 1 << 2 << 3; 00117 BM1.Row(2) << 4 << 5 << 6 << 7; 00118 BM1.Row(3) << 8 << 9 << 10 << 11 << 12; 00119 BM1.Row(4) << 13 << 14 << 15 << 16 << 17; 00120 BM1.Row(5) << 18 << 19 << 20 << 21; 00121 Matrix M1 = BM1 - BM; Print(M1); 00122 } 00123 { 00124 Tracer et1("Stage 5"); 00125 Matrix X(4,4); 00126 X << 1 << 2 << 3 << 4 00127 << 5 << 6 << 7 << 8 00128 << 9 <<10 <<11 <<12 00129 <<13 <<14 <<15 <<16; 00130 Matrix Y(4,0); 00131 Y = X | Y; 00132 X -= Y; Print(X); 00133 00134 DiagonalMatrix D(1); 00135 D << 23; // matrix input with just one value 00136 D(1) -= 23; Print(D); 00137 00138 } 00139 { 00140 Tracer et1("Stage 6"); 00141 Matrix h (2,2); 00142 h << 1.0 << 2.0 << 0.0 << 1.0 ; 00143 RowVector c(2); 00144 c << 0.0 << 1.0; 00145 h -= c & c; 00146 h -= c.t().Reverse() | c.Reverse().t(); 00147 Print(h); 00148 } 00149 { 00150 Tracer et1("Stage 7"); 00151 // Check row-wise input for diagonal matrix 00152 DiagonalMatrix D(4); 00153 D << 18 << 23 << 31 << 17; 00154 DiagonalMatrix D1(4); 00155 D1.Row(1) << 18; D1.Row(2) << 23; D1.Row(3) << 31; D1.Row(4) << 17; 00156 D1 -= D; Print(D1); 00157 D1(1) = 18; D1(2) = 23; D1(3) = 31; D1(4) = 17; 00158 D1 -= D; Print(D1); 00159 } 00160 00161 { 00162 Tracer et1("Stage 8"); 00163 // test swap functions 00164 MultWithCarry MWC; 00165 Matrix A(3,4); Matrix B(5,6); 00166 FillWithValues(MWC, A); FillWithValues(MWC, B); 00167 Matrix A1 = A; Matrix B1 = B; A.Release(); B.Release(2); 00168 swap(A, B); 00169 int a = A.size() - B1.size(), b = B.size() - A1.size(); 00170 Matrix D = A - B1; Print(D); 00171 D = B - A1; Print(D); 00172 Print(B); // should be zero because of release 00173 D = A - B1; Print(D); 00174 Print(A); // now should be zero 00175 D.ReSize(1,2); D(1,1) = a; D(1,2) = b; Print(D); 00176 00177 A.ReSize(20,20); FillWithValues(MWC, A); 00178 00179 UpperTriangularMatrix UA; UA << A; UpperTriangularMatrix UA1 = UA; 00180 UpperTriangularMatrix UB; 00181 swap(UA, UB); Print(UA); UB -= UA1; Print(UB); 00182 00183 LowerTriangularMatrix LA; LA << A; LowerTriangularMatrix LA1 = LA; 00184 LowerTriangularMatrix LB; 00185 swap(LB, LA); Print(LA); LB -= LA1; Print(LB); 00186 00187 SymmetricMatrix SA; SA << A; SymmetricMatrix SA1 = SA; 00188 SymmetricMatrix SB; 00189 swap(SA, SB); Print(SA); SB -= SA1; Print(SB); 00190 00191 DiagonalMatrix DA; DA << A; DiagonalMatrix DA1 = DA; 00192 DiagonalMatrix DB; 00193 swap(DB, DA); Print(DA); DB -= DA1; Print(DB); 00194 00195 RowVector RVA = A.Row(1); RowVector RVA1 = RVA; 00196 RowVector RVB; 00197 swap(RVB, RVA); Print(RVA); RVB -= RVA1; Print(RVB); 00198 00199 ColumnVector CVA = A.Column(1); ColumnVector CVA1 = CVA; 00200 ColumnVector CVB; 00201 swap(CVA, CVB); Print(CVA); CVB -= CVA1; Print(CVB); 00202 00203 BandMatrix BA(20, 7, 4); BA.Inject(A); BandMatrix BA1 = BA; 00204 BandMatrix BB; 00205 swap(BA, BB); D = BA; Print(D); BB -= BA1; D = BB; Print(D); 00206 00207 LowerBandMatrix LBA(20, 6); LBA.Inject(A); LowerBandMatrix LBA1 = LBA; 00208 LowerBandMatrix LBB; 00209 swap(LBB, LBA); D = LBA; Print(D); LBB -= LBA1; D = LBB; Print(D); 00210 00211 UpperBandMatrix UBA(20, 9); UBA.Inject(A); UpperBandMatrix UBA1 = UBA; 00212 UpperBandMatrix UBB; 00213 swap(UBA, UBB); D = UBA; Print(D); UBB -= UBA1; D = UBB; Print(D); 00214 00215 SymmetricBandMatrix SBA(20, 4); SBA.Inject(A); 00216 SymmetricBandMatrix SBA1 = SBA; 00217 SymmetricBandMatrix SBB; 00218 00219 swap(SBB, SBA); D = SBA; Print(D); 00220 SBB -= SBA1; D = SBB; Print(D); 00221 00222 B.ReSize(10,10); FillWithValues(MWC, B); 00223 00224 CroutMatrix CA = A; IdentityMatrix IA(20); 00225 CroutMatrix CB = B; IdentityMatrix IB(10); 00226 swap(CA, CB); swap(IA, IB); 00227 D = CA.i() * B - IA; Clean(D,0.00000001); Print(D); 00228 D = CB.i() * A - IB; Clean(D,0.00000001); Print(D); 00229 00230 BA.ReSize(20, 5, 7); BA.Inject(A); BandLUMatrix BLUA = BA; 00231 BB.ReSize(10, 3, 4); BB.Inject(B); BandLUMatrix BLUB = BB; 00232 swap(BLUA, BLUB); 00233 D = BLUA.i() * BB - IA; Clean(D,0.00000001); Print(D); 00234 D = BLUB.i() * BA - IB; Clean(D,0.00000001); Print(D); 00235 00236 00237 SBA.ReSize(20, 5); SBA.Inject(A); BandLUMatrix SBLUA = SBA; 00238 SBB.ReSize(10, 3); SBB.Inject(B); BandLUMatrix SBLUB = SBB; 00239 swap(SBLUA, SBLUB); 00240 D = SBLUA.i() * SBB - IA; Clean(D,0.00000001); Print(D); 00241 D = SBLUB.i() * SBA - IB; Clean(D,0.00000001); Print(D); 00242 00243 UA << A; 00244 GenericMatrix GUA = UA; GenericMatrix GB = B; swap(GUA, GB); 00245 D = GB - UA; Print(D); D = B - GUA; Print(D); 00246 00247 } 00248 00249 // cout << "\nEnd of fourth test\n"; 00250 } 00251 00252