tmtm.cpp
Go to the documentation of this file.
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 // test Kronecker Product
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");          // SquareMatrix
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13