00001
00002
00003
00006
00007
00008
00009
00010
00011 #include "include.h"
00012 #include "newmatap.h"
00013
00014 #include "tmt.h"
00015
00016 #ifdef use_namespace
00017 using namespace NEWMAT;
00018 #endif
00019
00020
00021
00022
00023 void trymatc()
00024 {
00025
00026 Tracer et("Twelfth test of Matrix package");
00027 Tracer::PrintTrace();
00028 DiagonalMatrix D(15); D=1.5;
00029 Matrix A(15,15);
00030 int i,j;
00031 for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
00032 { A = A + D; }
00033 ColumnVector B(15);
00034 for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
00035 {
00036 Tracer et1("Stage 1");
00037 ColumnVector B1=B;
00038 B=(A*2.0).i() * B1;
00039 Matrix X = A*B-B1/2.0;
00040 Clean(X, 0.000000001); Print(X);
00041 A.ReSize(3,5);
00042 for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
00043
00044 B = A.AsColumn()+10000;
00045 RowVector R = (A+10000).AsColumn().t();
00046 Print( RowVector(R-B.t()) );
00047 }
00048
00049 {
00050 Tracer et1("Stage 2");
00051 B = A.AsColumn()+10000;
00052 Matrix XR = (A+10000).AsMatrix(15,1).t();
00053 Print( RowVector(XR-B.t()) );
00054 }
00055
00056 {
00057 Tracer et1("Stage 3");
00058 B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
00059 Matrix MR = (A+10000).AsColumn().t();
00060 Print( RowVector(MR-B.t()) );
00061
00062 B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
00063 MR = A.AsColumn().t();
00064 Print( RowVector(MR-B.t()) );
00065 }
00066
00067 {
00068 Tracer et1("Stage 4");
00069 B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
00070 RowVector R = A.AsColumn().t();
00071 Print( RowVector(R-B.t()) );
00072 }
00073
00074 {
00075 Tracer et1("Stage 5");
00076 RowVector R = (A.AsColumn()-5000).t();
00077 B = ((R.t()+10000) - A.AsColumn())-5000;
00078 Print( RowVector(B.t()) );
00079 }
00080
00081 {
00082 Tracer et1("Stage 6");
00083 B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
00084 Print(ColumnVector(B1-B));
00085 }
00086
00087 {
00088 Tracer et1("Stage 7");
00089 Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
00090 for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
00091 Print(B);
00092 }
00093
00094 {
00095 Tracer et1("Stage 8");
00096 A.ReSize(7,7); D.ReSize(7);
00097 for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
00098 for (i=1; i<=7; i++) D(i,i) = i;
00099 UpperTriangularMatrix U; U << A;
00100 Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
00101 A.Inject(D); Print(Matrix(X-A));
00102 X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
00103 Print(Matrix(X-A));
00104 }
00105
00106 {
00107 Tracer et1("Stage 9");
00108 A.ReSize(7,5);
00109 for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
00110 Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
00111 Matrix X = A;
00112 Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
00113 Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
00114 }
00115
00116 {
00117 Tracer et1("Stage 10");
00118
00119 UpperTriangularMatrix U(20);
00120 for (i=1; i<=20; i++) for (j=i; j<=20; j++) U(i,j)=100 * i + j;
00121 UpperTriangularMatrix V = U.SymSubMatrix(1,5);
00122 UpperTriangularMatrix U1 = U;
00123 U1.SubMatrix(4,8,5,9) /= 2;
00124 U1.SubMatrix(4,8,5,9) += 388 * V;
00125 U1.SubMatrix(4,8,5,9) *= 2;
00126 U1.SubMatrix(4,8,5,9) += V;
00127 U1 -= U; UpperTriangularMatrix U2 = U1;
00128 U1 << U1.SubMatrix(4,8,5,9);
00129 U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
00130 U1 -= (777*V); Print(U1);
00131
00132 U1 = U; U1.SubMatrix(4,8,5,9) -= U.SymSubMatrix(1,5);
00133 U1 -= U; U2 = U1; U1 << U1.SubMatrix(4,8,5,9);
00134 U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
00135 U1 += V; Print(U1);
00136
00137 U1 = U;
00138 U1.SubMatrix(3,10,15,19) += 29;
00139 U1 -= U;
00140 Matrix X = U1.SubMatrix(3,10,15,19); X -= 29; Print(X);
00141 U1.SubMatrix(3,10,15,19) *= 0; Print(U1);
00142
00143 LowerTriangularMatrix L = U.t();
00144 LowerTriangularMatrix M = L.SymSubMatrix(1,5);
00145 LowerTriangularMatrix L1 = L;
00146 L1.SubMatrix(5,9,4,8) /= 2;
00147 L1.SubMatrix(5,9,4,8) += 388 * M;
00148 L1.SubMatrix(5,9,4,8) *= 2;
00149 L1.SubMatrix(5,9,4,8) += M;
00150 L1 -= L; LowerTriangularMatrix L2 = L1;
00151 L1 << L1.SubMatrix(5,9,4,8);
00152 L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
00153 L1 -= (777*M); Print(L1);
00154
00155 L1 = L; L1.SubMatrix(5,9,4,8) -= L.SymSubMatrix(1,5);
00156 L1 -= L; L2 =L1; L1 << L1.SubMatrix(5,9,4,8);
00157 L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
00158 L1 += M; Print(L1);
00159
00160 L1 = L;
00161 L1.SubMatrix(15,19,3,10) -= 29;
00162 L1 -= L;
00163 X = L1.SubMatrix(15,19,3,10); X += 29; Print(X);
00164 L1.SubMatrix(15,19,3,10) *= 0; Print(L1);
00165 }
00166
00167 {
00168 Tracer et1("Stage 11");
00169
00170 Matrix M(20,30);
00171 for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
00172 Matrix M1 = M;
00173
00174 for (j=1; j<=30; j++)
00175 { ColumnVector CV = 3 * M1.Column(j); M.Column(j) += CV; }
00176 for (i=1; i<=20; i++)
00177 { RowVector RV = 5 * M1.Row(i); M.Row(i) -= RV; }
00178
00179 M += M1; Print(M);
00180
00181 }
00182
00183 {
00184 Tracer et1("Stage 12");
00185
00186 Matrix M(20,30);
00187 for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
00188 Matrix M1 = M;
00189 M.Release();
00190 Matrix M2 = M;
00191 Matrix X = M; Print(X);
00192 X = M1 - M2; Print(X);
00193
00194 #ifndef DONT_DO_NRIC
00195 nricMatrix N = M1;
00196 nricMatrix N1 = N;
00197 N.Release();
00198 nricMatrix N2 = N;
00199 nricMatrix Y = N; Print(Y);
00200 Y = N1 - N2; Print(Y);
00201
00202 N = M1 / 2; N1 = N * 2; N.Release(); N2 = N * 2; Y = N; Print(N);
00203 Y = (N1 - M1) | (N2 - M1); Print(Y);
00204 #endif
00205
00206 }
00207
00208 {
00209 Tracer et("Stage 13");
00210
00211 MultWithCarry mwc;
00212 DiagonalMatrix DM; Matrix X;
00213
00214 Matrix A(20, 15);
00215 FillWithValues(mwc, A);
00216
00217 DM << A * A.t();
00218 ColumnVector CV = A.sum_square_rows();
00219 X = CV - DM.AsColumn(); Clean(X, 0.000000001); Print(X);
00220 DM << A.t() * A;
00221 RowVector RV = A.sum_square_columns();
00222 X = RV - DM.AsRow(); Clean(X, 0.000000001); Print(X);
00223 X = RV - A.t().sum_square_rows().t(); Clean(X, 0.000000001); Print(X);
00224 X = CV - A.t().sum_square_columns().t(); Clean(X, 0.000000001); Print(X);
00225
00226 A.ReSize(17,17); FillWithValues(mwc, A);
00227 UpperTriangularMatrix UT; UT << A;
00228 Matrix A1 = UT;
00229 X = UT.sum_square_rows() - A1.sum_square_rows(); Print(X);
00230 X = UT.sum_square_columns() - A1.sum_square_columns(); Print(X);
00231
00232 LowerTriangularMatrix LT; LT << A;
00233 A1 = LT;
00234 X = LT.sum_square_rows() - A1.sum_square_rows(); Print(X);
00235 X = LT.sum_square_columns() - A1.sum_square_columns(); Print(X);
00236
00237 SymmetricMatrix SM; SM << A;
00238 A1 = SM;
00239 X = SM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00240 X = SM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00241
00242 DM << A;
00243 A1 = DM;
00244 X = DM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00245 X = DM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00246
00247 BandMatrix BM(17, 3, 5); BM.Inject(A);
00248 A1 = BM;
00249 X = BM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00250 X = BM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00251
00252 SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
00253 A1 = SBM;
00254 X = SBM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00255 X = SBM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00256
00257 IdentityMatrix IM(29);
00258 X = IM.sum_square_rows() - 1; Print(X);
00259 X = IM.sum_square_columns() - 1; Print(X);
00260
00261 A1.ReSize(0,10);
00262 X.ReSize(1,10); X = 0; X -= A1.sum_square_columns(); Print(X);
00263 X.ReSize(0,1); X -= A1.sum_square_rows(); Print(X);
00264
00265 A1.ReSize(10,0);
00266 X.ReSize(10,1); X = 0; X -= A1.sum_square_rows(); Print(X);
00267 X.ReSize(1,0); X -= A1.sum_square_columns(); Print(X);
00268
00269 }
00270
00271 {
00272 Tracer et("Stage 14");
00273
00274 MultWithCarry mwc;
00275 Matrix A(20,5); FillWithValues(mwc, A);
00276
00277 UpperTriangularMatrix R;
00278 Matrix A_old = A;
00279 QRZ(A,R);
00280
00281 Matrix X = A * R - A_old; Clean(X, 0.000000001); Print(X);
00282
00283 X = A.t() * A - IdentityMatrix(5);
00284 Clean(X, 0.000000001); Print(X);
00285
00286 SquareMatrix A1(20);
00287 A1.Columns(1,5) = A;
00288 extend_orthonormal(A1,5);
00289
00290 X = A - A1.Columns(1,5); Print(X);
00291
00292 X = A1.t() * A1 - IdentityMatrix(20);
00293 Clean(X, 0.000000001); Print(X);
00294 X = A1 * A1.t() - IdentityMatrix(20);
00295 Clean(X, 0.000000001); Print(X);
00296
00297 Matrix A2(20,15);
00298 A2.Columns(1,5) = A;
00299 extend_orthonormal(A2,5);
00300
00301 X = A - A2.Columns(1,5); Print(X);
00302
00303 X = A2.t() * A2 - IdentityMatrix(15);
00304 Clean(X, 0.000000001); Print(X);
00305
00306 A2.ReSize(100,100);
00307 extend_orthonormal(A2,0);
00308
00309 X = A2.t() * A2 - IdentityMatrix(100);
00310 Clean(X, 0.000000001); Print(X);
00311 X = A2 * A2.t() - IdentityMatrix(100);
00312 Clean(X, 0.000000001); Print(X);
00313 }
00314
00315 {
00316 Tracer et("Stage 15");
00317
00318 MultWithCarry mwc;
00319 Matrix A(15,40); FillWithValues(mwc, A);
00320 Matrix U, V; DiagonalMatrix D;
00321 SVD(A.t(),D,V,U);
00322 Matrix VE(40,40); VE.columns(1,15) = V;
00323 extend_orthonormal(VE,15);
00324 Matrix DE(15,40); DE = 0; DE.columns(1,15) = D;
00325 Matrix B = U * DE * VE.t();
00326 B -= A;
00327 Clean(B, 0.000000001); Print(B);
00328 }
00329
00330 {
00331 Tracer et("Stage 16");
00332
00333 MultWithCarry mwc;
00334 ColumnVector CVX; Matrix X;
00335
00336 Matrix A(20, 15);
00337 FillWithValues(mwc, A);
00338
00339 ColumnVector Ones(15); Ones = 1;
00340 CVX = A * Ones;
00341 ColumnVector CV = A.sum_rows();
00342 X = CV - CVX; Clean(X, 0.000000001); Print(X);
00343 Ones.resize(20); Ones = 1;
00344 CVX << A.t() * Ones;
00345 RowVector RV = A.sum_columns();
00346 X = RV - CVX.AsRow(); Clean(X, 0.000000001); Print(X);
00347 X = RV - A.t().sum_rows().t(); Clean(X, 0.000000001); Print(X);
00348 X = CV - A.t().sum_columns().t(); Clean(X, 0.000000001); Print(X);
00349
00350 A.ReSize(17,17); FillWithValues(mwc, A);
00351 UpperTriangularMatrix UT; UT << A;
00352 Matrix A1 = UT;
00353 X = UT.sum_rows() - A1.sum_rows(); Print(X);
00354 X = UT.sum_columns() - A1.sum_columns(); Print(X);
00355
00356 LowerTriangularMatrix LT; LT << A;
00357 A1 = LT;
00358 X = LT.sum_rows() - A1.sum_rows(); Print(X);
00359 X = LT.sum_columns() - A1.sum_columns(); Print(X);
00360
00361 SymmetricMatrix SM; SM << A;
00362 A1 = SM;
00363 X = SM.sum_rows() - A1.sum_rows(); Print(X);
00364 X = SM.sum_columns() - A1.sum_columns(); Print(X);
00365
00366 DiagonalMatrix DM; DM << A;
00367 A1 = DM;
00368 X = DM.sum_rows() - A1.sum_rows(); Print(X);
00369 X = DM.sum_columns() - A1.sum_columns(); Print(X);
00370
00371 BandMatrix BM(17, 3, 5); BM.Inject(A);
00372 A1 = BM;
00373 X = BM.sum_rows() - A1.sum_rows(); Print(X);
00374 X = BM.sum_columns() - A1.sum_columns(); Print(X);
00375
00376 SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
00377 A1 = SBM;
00378 X = SBM.sum_rows() - A1.sum_rows(); Print(X);
00379 X = SBM.sum_columns() - A1.sum_columns(); Print(X);
00380
00381 IdentityMatrix IM(29);
00382 X = IM.sum_rows() - 1; Print(X);
00383 X = IM.sum_columns() - 1; Print(X);
00384
00385 A1.ReSize(0,10);
00386 X.ReSize(1,10); X = 0; X -= A1.sum_columns(); Print(X);
00387 X.ReSize(0,1); X -= A1.sum_rows(); Print(X);
00388
00389 A1.ReSize(10,0);
00390 X.ReSize(10,1); X = 0; X -= A1.sum_rows(); Print(X);
00391 X.ReSize(1,0); X -= A1.sum_columns(); Print(X);
00392
00393 }
00394
00395
00396
00397
00398 }
00399
00400