00001
00002
00003
00006
00007
00008 #define WANT_STREAM
00009
00010 #define WANT_MATH
00011
00012 #include "newmat.h"
00013 #include "newmatio.h"
00014
00015 #include "tmt.h"
00016
00017 #ifdef use_namespace
00018 using namespace NEWMAT;
00019 #endif
00020
00021
00022
00023
00024
00025
00026 void trymatm()
00027 {
00028 Tracer et("Twenty second test of Matrix package");
00029 Tracer::PrintTrace();
00030
00031 {
00032 Tracer et1("Stage 1");
00033
00034
00035 Matrix A(2,3);
00036 A << 3 << 5 << 2
00037 << 4 << 1 << 6;
00038
00039 Matrix B(4,3);
00040 B << 7 << 2 << 9
00041 << 1 << 3 << 6
00042 << 4 << 10 << 5
00043 << 11 << 8 << 12;
00044
00045 Matrix C(8, 9);
00046
00047 C.Row(1) << 21 << 6 << 27 << 35 << 10 << 45 << 14 << 4 << 18;
00048 C.Row(2) << 3 << 9 << 18 << 5 << 15 << 30 << 2 << 6 << 12;
00049 C.Row(3) << 12 << 30 << 15 << 20 << 50 << 25 << 8 << 20 << 10;
00050 C.Row(4) << 33 << 24 << 36 << 55 << 40 << 60 << 22 << 16 << 24;
00051
00052 C.Row(5) << 28 << 8 << 36 << 7 << 2 << 9 << 42 << 12 << 54;
00053 C.Row(6) << 4 << 12 << 24 << 1 << 3 << 6 << 6 << 18 << 36;
00054 C.Row(7) << 16 << 40 << 20 << 4 << 10 << 5 << 24 << 60 << 30;
00055 C.Row(8) << 44 << 32 << 48 << 11 << 8 << 12 << 66 << 48 << 72;
00056
00057 Matrix AB = KP(A,B) - C; Print(AB);
00058
00059 IdentityMatrix I1(10); IdentityMatrix I2(15); I2 *= 2;
00060 DiagonalMatrix D = KP(I1, I2) - IdentityMatrix(150) * 2;
00061 Print(D);
00062 }
00063
00064 {
00065 Tracer et1("Stage 2");
00066
00067 UpperTriangularMatrix A(3);
00068 A << 3 << 8 << 5
00069 << 7 << 2
00070 << 4;
00071 UpperTriangularMatrix B(4);
00072 B << 4 << 1 << 7 << 2
00073 << 3 << 9 << 8
00074 << 1 << 5
00075 << 6;
00076
00077 UpperTriangularMatrix C(12);
00078
00079 C.Row(1) <<12<< 3<<21<< 6 <<32<< 8<<56<<16 <<20<< 5<<35<<10;
00080 C.Row(2) << 9<<27<<24 << 0<<24<<72<<64 << 0<<15<<45<<40;
00081 C.Row(3) << 3<<15 << 0<< 0<< 8<<40 << 0<< 0<< 5<<25;
00082 C.Row(4) <<18 << 0<< 0<< 0<<48 << 0<< 0<< 0<<30;
00083
00084 C.Row(5) <<28<< 7<<49<<14 << 8<< 2<<14<< 4;
00085 C.Row(6) <<21<<63<<56 << 0<< 6<<18<<16;
00086 C.Row(7) << 7<<35 << 0<< 0<< 2<<10;
00087 C.Row(8) <<42 << 0<< 0<< 0<<12;
00088
00089 C.Row(9) <<16<< 4<<28<< 8;
00090 C.Row(10) <<12<<36<<32;
00091 C.Row(11) << 4<<20;
00092 C.Row(12) <<24;
00093
00094
00095 UpperTriangularMatrix AB = KP(A,B) - C; Print(AB);
00096
00097 LowerTriangularMatrix BT = B.t(); Matrix N(12,12);
00098
00099 N.Row(1) <<12 << 0<< 0<< 0 <<32<< 0<< 0<< 0 <<20<< 0<< 0<< 0;
00100 N.Row(2) << 3 << 9<< 0<< 0 << 8<<24<< 0<< 0 << 5<<15<< 0<< 0;
00101 N.Row(3) <<21 <<27<< 3<< 0 <<56<<72<< 8<< 0 <<35<<45<< 5<< 0;
00102 N.Row(4) << 6 <<24<<15<<18 <<16<<64<<40<<48 <<10<<40<<25<<30;
00103
00104 N.Row(5) << 0 << 0<< 0<< 0 <<28<< 0<< 0<< 0 << 8<< 0<< 0<< 0;
00105 N.Row(6) << 0 << 0<< 0<< 0 << 7<<21<< 0<< 0 << 2<< 6<< 0<< 0;
00106 N.Row(7) << 0 << 0<< 0<< 0 <<49<<63<< 7<< 0 <<14<<18<< 2<< 0;
00107 N.Row(8) << 0 << 0<< 0<< 0 <<14<<56<<35<<42 << 4<<16<<10<<12;
00108
00109 N.Row(9) << 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<16<< 0<< 0<< 0;
00110 N.Row(10)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 4<<12<< 0<< 0;
00111 N.Row(11)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<28<<36<< 4<< 0;
00112 N.Row(12)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 8<<32<<20<<24;
00113
00114 Matrix N1 = KP(A, BT); N1 -= N; Print(N1);
00115 AB << KP(A, BT); AB << (AB - N); Print(AB);
00116 BT << KP(A, BT); BT << (BT - N); Print(BT);
00117
00118 LowerTriangularMatrix AT = A.t();
00119 N1 = KP(AT, B); N1 -= N.t(); Print(N1);
00120 AB << KP(AT, B); AB << (AB - N.t()); Print(AB);
00121 BT << KP(AT, B); BT << (BT - N.t()); Print(BT);
00122 }
00123
00124 {
00125 Tracer et1("Stage 3");
00126
00127 BandMatrix BMA(6,2,3);
00128 BMA.Row(1) << 5.25 << 4.75 << 2.25 << 1.75;
00129 BMA.Row(2) << 1.25 << 9.75 << 4.50 << 0.25 << 1.50;
00130 BMA.Row(3) << 7.75 << 1.50 << 3.00 << 4.25 << 0.50 << 5.50;
00131 BMA.Row(4) << 2.75 << 9.00 << 8.00 << 3.25 << 3.50;
00132 BMA.Row(5) << 8.75 << 6.25 << 5.00 << 5.75;
00133 BMA.Row(6) << 3.75 << 6.75 << 6.00;
00134
00135 Matrix A = BMA;
00136
00137 BandMatrix BMB(4,2,1);
00138 BMB.Row(1) << 4.5 << 9.5;
00139 BMB.Row(2) << 1.5 << 6.0 << 2.0;
00140 BMB.Row(3) << 0.5 << 2.5 << 8.5 << 7.5;
00141 BMB.Row(4) << 3.0 << 4.0 << 6.5;
00142
00143 SquareMatrix B = BMB;
00144
00145 BandMatrix BMC = KP(BMA, BMB);
00146 BandMatrix BMC1 = KP(BMA, B);
00147 Matrix C2 = KP(A, BMB);
00148 Matrix C = KP(A, B);
00149
00150 Matrix M = C - BMC; Print(M);
00151 M = C - BMC1; Print(M);
00152 M = C - C2; Print(M);
00153
00154 RowVector X(4);
00155 X(1) = BMC.BandWidth().Lower() - 10;
00156 X(2) = BMC.BandWidth().Upper() - 13;
00157 X(3) = BMC1.BandWidth().Lower() - 11;
00158 X(4) = BMC1.BandWidth().Upper() - 15;
00159 Print(X);
00160
00161 UpperTriangularMatrix UT; UT << KP(BMA, BMB);
00162 UpperTriangularMatrix UT1; UT1 << (C - UT); Print(UT1);
00163 LowerTriangularMatrix LT; LT << KP(BMA, BMB);
00164 LowerTriangularMatrix LT1; LT1 << (C - LT); Print(LT1);
00165 }
00166
00167 {
00168 Tracer et1("Stage 4");
00169
00170 SymmetricMatrix SM1(4);
00171 SM1.Row(1) << 2;
00172 SM1.Row(2) << 4 << 5;
00173 SM1.Row(3) << 9 << 2 << 1;
00174 SM1.Row(4) << 3 << 6 << 8 << 2;
00175
00176 SymmetricMatrix SM2(3);
00177 SM2.Row(1) << 3;
00178 SM2.Row(2) << -7 << -6;
00179 SM2.Row(3) << 4 << -2 << -1;
00180
00181 SymmetricMatrix SM = KP(SM1, SM2);
00182 Matrix M1 = SM1; Matrix M2 = SM2;
00183 Matrix M = KP(SM1, SM2); M -= SM; Print(M);
00184 M = KP(SM1, SM2) - SM; Print(M);
00185 M = KP(M1, SM2) - SM; Print(M);
00186 M = KP(SM1, M2) - SM; Print(M);
00187 M = KP(M1, M2); M -= SM; Print(M);
00188 }
00189
00190 {
00191 Tracer et1("Stage 5");
00192
00193 Matrix A(2,3);
00194 A << 3 << 5 << 2
00195 << 4 << 1 << 6;
00196
00197 Matrix B(3,4);
00198 B << 7 << 2 << 9 << 11
00199 << 1 << 3 << 6 << 8
00200 << 4 << 10 << 5 << 12;
00201
00202 RowVector C(2); C << 3 << 7;
00203 ColumnVector D(4); D << 0 << 5 << 13 << 11;
00204
00205 Matrix M = KP(C * A, B * D) - KP(C, B) * KP(A, D); Print(M);
00206 }
00207
00208 {
00209 Tracer et1("Stage 6");
00210
00211 RowVector A(3), B(5), C(15);
00212 A << 5 << 2 << 4;
00213 B << 3 << 2 << 0 << 1 << 6;
00214 C << 15 << 10 << 0 << 5 << 30
00215 << 6 << 4 << 0 << 2 << 12
00216 << 12 << 8 << 0 << 4 << 24;
00217 Matrix N = KP(A, B) - C; Print(N);
00218 N = KP(A.t(), B.t()) - C.t(); Print(N);
00219 N = KP(A.AsDiagonal(), B.AsDiagonal()) - C.AsDiagonal(); Print(N);
00220 }
00221
00222 {
00223 Tracer et1("Stage 7");
00224 IdentityMatrix I(3);
00225 ColumnVector CV(4); CV << 4 << 3 << 1 << 7;
00226 Matrix A = KP(I, CV) + 5;
00227 Matrix B(3,12);
00228 B.Row(1) << 9 << 8 << 6 << 12 << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 5;
00229 B.Row(2) << 5 << 5 << 5 << 5 << 9 << 8 << 6 << 12 << 5 << 5 << 5 << 5;
00230 B.Row(3) << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 9 << 8 << 6 << 12;
00231 B -= A.t(); Print(B);
00232
00233 }
00234
00235 {
00236 Tracer et1("Stage 8");
00237 Matrix A(2,3), B(3,2);
00238 A << 2 << 6 << 7
00239 << 4 << 3 << 9;
00240 B << 1 << 3
00241 << 4 << 8
00242 << 0 << 6;
00243 SquareMatrix AB = A * B;
00244 Matrix M = (B.t() * A.t()).t(); M -= AB; Print(M);
00245 AB = B * A;
00246 M = (A.t() * B.t()).t(); M -= AB; Print(M);
00247 AB.ReSize(5,5); AB = 0;
00248 AB.SubMatrix(1,2,1,3) = A; AB.SubMatrix(4,5,3,5) = A;
00249 AB.SubMatrix(1,3,4,5) = B; AB.SubMatrix(3,5,1,2) = B;
00250 SquareMatrix C(5);
00251 C.Row(1) << 2 << 6 << 7 << 1 << 3;
00252 C.Row(2) << 4 << 3 << 9 << 4 << 8;
00253 C.Row(3) << 1 << 3 << 0 << 0 << 6;
00254 C.Row(4) << 4 << 8 << 2 << 6 << 7;
00255 C.Row(5) << 0 << 6 << 4 << 3 << 9;
00256 C -= AB; Print(C);
00257 AB = A.SymSubMatrix(1,2);
00258 AB = (AB | AB) & (AB | AB);
00259 C.ReSize(4);
00260 C.Row(1) << 2 << 6 << 2 << 6;
00261 C.Row(2) << 4 << 3 << 4 << 3;
00262 C.Row(3) << 2 << 6 << 2 << 6;
00263 C.Row(4) << 4 << 3 << 4 << 3;
00264 M = AB;
00265 C -= M; Print(C);
00266 C << M; C += -M; Print(C);
00267
00268 }
00269
00270
00271 }
00272
00273
00275
00276
00277
00278
00279