$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 trymat2() 00026 { 00027 // cout << "\nSecond test of Matrix package\n\n"; 00028 Tracer et("Second test of Matrix package"); 00029 Tracer::PrintTrace(); 00030 00031 int i,j; 00032 00033 Matrix M(3,5); 00034 for (i=1; i<=3; i++) for (j=1; j<=5; j++) M(i,j) = 100*i + j; 00035 Matrix X(8,10); 00036 for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j; 00037 Matrix Y = X; Matrix Z = X; 00038 { X.SubMatrix(2,4,3,7) << M; } 00039 for (i=1; i<=3; i++) for (j=1; j<=5; j++) Y(i+1,j+2) = 100*i + j; 00040 Print(Matrix(X-Y)); 00041 00042 00043 Real a[15]; Real* r = a; 00044 for (i=1; i<=3; i++) for (j=1; j<=5; j++) *r++ = 100*i + j; 00045 { Z.SubMatrix(2,4,3,7) << a; } 00046 Print(Matrix(Z-Y)); 00047 00048 { M=33; X.SubMatrix(2,4,3,7) << M; } 00049 { Z.SubMatrix(2,4,3,7) = 33; } 00050 Print(Matrix(Z-X)); 00051 00052 for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j; 00053 Y = X; 00054 UpperTriangularMatrix U(5); 00055 for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j; 00056 { X.SubMatrix(3,7,5,9) << U; } 00057 for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j; 00058 for (i=1; i<=5; i++) for (j=1; j<i; j++) Y(i+2,j+4) = 0.0; 00059 Print(Matrix(X-Y)); 00060 for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j; 00061 Y = X; 00062 for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j; 00063 { X.SubMatrix(3,7,5,9).Inject(U); } 00064 for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j; 00065 Print(Matrix(X-Y)); 00066 00067 00068 // test growing and shrinking a vector 00069 { 00070 ColumnVector V(100); 00071 for (i=1;i<=100;i++) V(i) = i*i+i; 00072 V = V.Rows(1,50); // to get first 50 vlaues. 00073 00074 { 00075 V.Release(); ColumnVector VX=V; 00076 V.ReSize(100); V = 0.0; V.Rows(1,50)=VX; 00077 } // V now length 100 00078 00079 M=V; M=100; // to make sure V will hold its values 00080 for (i=1;i<=50;i++) V(i) -= i*i+i; 00081 Print(V); 00082 00083 00084 // test redimensioning vectors with two dimensions given 00085 ColumnVector CV1(10); CV1 = 10; 00086 ColumnVector CV2(5); CV2.ReSize(10,1); CV2 = 10; 00087 V = CV1-CV2; Print(V); 00088 00089 RowVector RV1(20); RV1 = 100; 00090 RowVector RV2; RV2.ReSize(1,20); RV2 = 100; 00091 V = (RV1-RV2).t(); Print(V); 00092 00093 X.ReSize(4,7); 00094 for (i=1; i<=4; i++) for (j=1; j<=7; j++) X(i,j) = 1000*i + 10*j; 00095 Y = 10.5 * X; 00096 Z = 7.25 - Y; 00097 M = Z + X * 10.5 - 7.25; 00098 Print(M); 00099 Y = 2.5 * X; 00100 Z = 9.25 + Y; 00101 M = Z - X * 2.5 - 9.25; 00102 Print(M); 00103 U.ReSize(8); 00104 for (i=1; i<=8; i++) for (j=i; j<=8; j++) U(i,j) = 100*i + j; 00105 Y = 100 - U; 00106 M = Y + U - 100; 00107 Print(M); 00108 } 00109 00110 { 00111 SymmetricMatrix S,T; 00112 00113 S << (U + U.t()); 00114 T = 100 - S; M = T + S - 100; Print(M); 00115 T = 100 - 2 * S; M = T + S * 2 - 100; Print(M); 00116 X = 100 - 2 * S; M = X + S * 2 - 100; Print(M); 00117 T = S; T = 100 - T; M = T + S - 100; Print(M); 00118 } 00119 00120 // test new 00121 { 00122 ColumnVector CV1; RowVector RV1; 00123 Matrix* MX; MX = new Matrix; if (!MX) Throw(Bad_alloc("New fails ")); 00124 MX->ReSize(10,20); 00125 for (i = 1; i <= 10; i++) for (j = 1; j <= 20; j++) 00126 (*MX)(i,j) = 100 * i + j; 00127 ColumnVector* CV = new ColumnVector(10); 00128 if (!CV) Throw(Bad_alloc("New fails ")); 00129 *CV << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9 << 10; 00130 RowVector* RV = new RowVector(CV->t() | (*CV + 10).t()); 00131 if (!RV) Throw(Bad_alloc("New fails ")); 00132 CV1 = ColumnVector(10); CV1 = 1; RV1 = RowVector(20); RV1 = 1; 00133 *MX -= 100 * *CV * RV1 + CV1 * *RV; 00134 Print(*MX); 00135 delete MX; delete CV; delete RV; 00136 } 00137 00138 00139 // test copying of vectors and matrices with no elements 00140 { 00141 ColumnVector dims(16); 00142 Matrix M1; Matrix M2 = M1; Print(M2); 00143 dims(1) = M2.Nrows(); dims(2) = M2.Ncols(); 00144 dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage(); 00145 M2 = M1; 00146 dims(5) = M2.Nrows(); dims(6) = M2.Ncols(); 00147 dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage(); 00148 M2.ReSize(10,20); M2.CleanUp(); 00149 dims(9) = M2.Nrows(); dims(10) = M2.Ncols(); 00150 dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage(); 00151 M2.ReSize(20,10); M2.ReSize(0,0); 00152 dims(13) = M2.Nrows(); dims(14) = M2.Ncols(); 00153 dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage(); 00154 Print(dims); 00155 } 00156 00157 { 00158 ColumnVector dims(16); 00159 ColumnVector M1; ColumnVector M2 = M1; Print(M2); 00160 dims(1) = M2.Nrows(); dims(2) = M2.Ncols()-1; 00161 dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage(); 00162 M2 = M1; 00163 dims(5) = M2.Nrows(); dims(6) = M2.Ncols()-1; 00164 dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage(); 00165 M2.ReSize(10); M2.CleanUp(); 00166 dims(9) = M2.Nrows(); dims(10) = M2.Ncols()-1; 00167 dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage(); 00168 M2.ReSize(10); M2.ReSize(0); 00169 dims(13) = M2.Nrows(); dims(14) = M2.Ncols()-1; 00170 dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage(); 00171 Print(dims); 00172 } 00173 00174 { 00175 ColumnVector dims(16); 00176 RowVector M1; RowVector M2 = M1; Print(M2); 00177 dims(1) = M2.Nrows()-1; dims(2) = M2.Ncols(); 00178 dims(3) = (Real)(unsigned long)M2.Store(); dims(4) = M2.Storage(); 00179 M2 = M1; 00180 dims(5) = M2.Nrows()-1; dims(6) = M2.Ncols(); 00181 dims(7) = (Real)(unsigned long)M2.Store(); dims(8) = M2.Storage(); 00182 M2.ReSize(10); M2.CleanUp(); 00183 dims(9) = M2.Nrows()-1; dims(10) = M2.Ncols(); 00184 dims(11) = (Real)(unsigned long)M2.Store(); dims(12) = M2.Storage(); 00185 M2.ReSize(10); M2.ReSize(0); 00186 dims(13) = M2.Nrows()-1; dims(14) = M2.Ncols(); 00187 dims(15) = (Real)(unsigned long)M2.Store(); dims(16) = M2.Storage(); 00188 Print(dims); 00189 } 00190 00191 // test identity matrix 00192 { 00193 Matrix M; 00194 IdentityMatrix I(10); DiagonalMatrix D(10); D = 1; 00195 M = I; M -= D; Print(M); 00196 D -= I; Print(D); 00197 ColumnVector X(8); 00198 D = 1; 00199 X(1) = Sum(D) - Sum(I); 00200 X(2) = SumAbsoluteValue(D) - SumAbsoluteValue(I); 00201 X(3) = SumSquare(D) - SumSquare(I); 00202 X(4) = Trace(D) - Trace(I); 00203 X(5) = Maximum(D) - Maximum(I); 00204 X(6) = Minimum(D) - Minimum(I); 00205 X(7) = LogDeterminant(D).LogValue() - LogDeterminant(I).LogValue(); 00206 X(8) = LogDeterminant(D).Sign() - LogDeterminant(I).Sign(); 00207 Clean(X,0.00000001); Print(X); 00208 00209 for (i = 1; i <= 10; i++) for (j = 1; j <= 10; j++) 00210 M(i,j) = 100 * i + j; 00211 Matrix N; 00212 N = M * I - M; Print(N); 00213 N = I * M - M; Print(N); 00214 N = M * I.i() - M; Print(N); 00215 N = I.i() * M - M; Print(N); 00216 N = I.i(); N -= I; Print(N); 00217 N = I.t(); N -= I; Print(N); 00218 N = I.t(); N += (-I); Print(N); // <---------------- 00219 D = I; N = D; D = 1; N -= D; Print(N); 00220 N = I; D = 1; N -= D; Print(N); 00221 N = M + 2 * IdentityMatrix(10); N -= (M + 2 * D); Print(N); 00222 00223 I *= 4; 00224 00225 D = 4; 00226 00227 X.ReSize(14); 00228 X(1) = Sum(D) - Sum(I); 00229 X(2) = SumAbsoluteValue(D) - SumAbsoluteValue(I); 00230 X(3) = SumSquare(D) - SumSquare(I); 00231 X(4) = Trace(D) - Trace(I); 00232 X(5) = Maximum(D) - Maximum(I); 00233 X(6) = Minimum(D) - Minimum(I); 00234 X(7) = LogDeterminant(D).LogValue() - LogDeterminant(I).LogValue(); // <-- 00235 X(8) = LogDeterminant(D).Sign() - LogDeterminant(I).Sign(); 00236 int i,j; 00237 X(9) = I.Maximum1(i) - 4; X(10) = i-1; 00238 X(11) = I.Maximum2(i,j) - 4; X(12) = i-10; X(13) = j-10; 00239 X(14) = I.Nrows() - 10; 00240 Clean(X,0.00000001); Print(X); 00241 00242 00243 N = D.i(); 00244 N += I / (-16); 00245 Print(N); 00246 N = M * I - 4 * M; Print(N); 00247 N = I * M - 4 * M; Print(N); 00248 N = M * I.i() - 0.25 * M; Print(N); 00249 N = I.i() * M - 0.25 * M; Print(N); 00250 N = I.i(); N -= I * 0.0625; Print(N); 00251 N = I.i(); N = N - 0.0625 * I; Print(N); 00252 N = I.t(); N -= I; Print(N); 00253 D = I * 2; N = D; D = 1; N -= 8 * D; Print(N); 00254 N = I * 2; N -= 8 * D; Print(N); 00255 N = 0.5 * I + M; N -= M; N -= 2.0 * D; Print(N); 00256 00257 IdentityMatrix J(10); J = 8; 00258 D = 4; 00259 DiagonalMatrix E(10); E = 8; 00260 N = (I + J) - (D + E); Print(N); 00261 N = (5*I + 3*J) - (5*D + 3*E); Print(N); 00262 N = (-I + J) - (-D + E); Print(N); 00263 N = (I - J) - (D - E); Print(N); 00264 N = (I | J) - (D | E); Print(N); 00265 N = (I & J) - (D & E); Print(N); 00266 N = SP(I,J) - SP(D,E); Print(N); 00267 N = D.SubMatrix(2,5,3,8) - I.SubMatrix(2,5,3,8); Print(N); 00268 00269 N = M; N.Inject(I); D << M; N -= (M + I); N += D; Print(N); 00270 D = 4; 00271 00272 IdentityMatrix K = I.i()*7 - J.t()/4; 00273 N = D.i() * 7 - E / 4 - K; Print(N); 00274 K = I * J; N = K - D * E; Print(N); 00275 N = I * J; N -= D * E; Print(N); 00276 K = 5*I - 3*J; 00277 N = K - (5*D - 3*E); Print(N); 00278 K = I.i(); N = K - 0.0625 * I; Print(N); 00279 K = I.t(); N = K - I; Print(N); 00280 00281 00282 K.ReSize(20); D.ReSize(20); D = 1; 00283 D -= K; Print(D); 00284 00285 I.ReSize(3); J.ReSize(3); K = I * J; N = K - I; Print(N); 00286 K << D; N = K - D; Print(N); 00287 } 00288 00289 // test add integer 00290 { 00291 Matrix X(2,3); 00292 X << 5.25 << 7.75 << 1.25 00293 << 9.00 << 1.00 << 2.50; 00294 Matrix Y = X; 00295 X = 10 + X; 00296 X += (-10); 00297 X -= Y; 00298 Print(X); 00299 00300 // also test f suffix 00301 X << 5.25f << 7.75f << 1.25f 00302 << 9.00f << 1.00f << 2.50f; 00303 X -= Y; Print(X); 00304 00305 } 00306 00307 00308 00309 00310 // cout << "\nEnd of second test\n"; 00311 } 00312 00313