00001
00002
00003
00006
00007
00008
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
00023
00024
00025 void trymat4()
00026 {
00027
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
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
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;
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
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
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);
00173 D = A - B1; Print(D);
00174 Print(A);
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
00250 }
00251
00252