00001
00002
00003
00006
00007
00008
00009
00010 #include "include.h"
00011
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
00024
00025 void Transposer(const GenericMatrix& GM1, GenericMatrix&GM2)
00026 { GM2 = GM1.t(); }
00027
00028
00029
00030
00031
00032 static void DCR(Real d[], Real c[], int m, Real r[], int n, Real **dcr)
00033 {
00034 int i, j;
00035 for (i = 1; i <= m; i++) for (j = 1; j <= n; j++)
00036 dcr[i][j] = d[i] * c[i] * r[j];
00037 }
00038
00039 ReturnMatrix TestReturn(const GeneralMatrix& gm) { return gm; }
00040
00041 void trymat8()
00042 {
00043
00044 Tracer et("Eighth test of Matrix package");
00045 Tracer::PrintTrace();
00046
00047 int i;
00048
00049
00050 DiagonalMatrix D(6);
00051 for (i=1;i<=6;i++) D(i,i)=i*i+i-10;
00052 DiagonalMatrix D2=D;
00053 Matrix MD=D;
00054
00055 DiagonalMatrix D1(6); for (i=1;i<=6;i++) D1(i,i)=-100+i*i*i;
00056 Matrix MD1=D1;
00057 Print(Matrix(D*D1-MD*MD1));
00058 Print(Matrix((-D)*D1+MD*MD1));
00059 Print(Matrix(D*(-D1)+MD*MD1));
00060 DiagonalMatrix DX=D;
00061 {
00062 Tracer et1("Stage 1");
00063 DX=(DX+D1)*DX; Print(Matrix(DX-(MD+MD1)*MD));
00064 DX=D;
00065 DX=-DX*DX+(DX-(-D1))*((-D1)+DX);
00066
00067
00068 MD1=DX+(Matrix(MD1).t())*(Matrix(MD1).t()); Print(MD1);
00069 DX=D; DX=DX; DX=D2-DX; Print(DiagonalMatrix(DX));
00070 DX=D;
00071 }
00072 {
00073 Tracer et1("Stage 2");
00074 D.Release(2);
00075 D1=D; D2=D;
00076 Print(DiagonalMatrix(D1-DX));
00077 Print(DiagonalMatrix(D2-DX));
00078 MD1=1.0;
00079 Print(Matrix(MD1-1.0));
00080 }
00081 {
00082 Tracer et1("Stage 3");
00083
00084 LowerTriangularMatrix LT(4);
00085 LT << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9 << 10;
00086 UpperTriangularMatrix UT = LT.t() * 2.0;
00087 GenericMatrix GM1 = LT;
00088 LowerTriangularMatrix LT1 = GM1-LT; Print(LT1);
00089 GenericMatrix GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
00090 GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
00091 GM2 = GM1*2; LT1 = GM2; LT1 = LT1-LT*2; Print(LT1);
00092 GM1.Release();
00093 GM1=GM1; LT1=GM1-LT; Print(LT1); LT1=GM1-LT; Print(LT1);
00094 GM1.Release();
00095 GM1=GM1*4; LT1=GM1-LT*4; Print(LT1);
00096 LT1=GM1-LT*4; Print(LT1); GM1.CleanUp();
00097 GM1=LT; GM2=UT; GM1=GM1*GM2; Matrix M=GM1; M=M-LT*UT; Print(M);
00098 Transposer(LT,GM2); LT1 = LT - GM2.t(); Print(LT1);
00099 GM1=LT; Transposer(GM1,GM2); LT1 = LT - GM2.t(); Print(LT1);
00100 GM1 = LT; GM1 = GM1 + GM1; LT1 = LT*2-GM1; Print(LT1);
00101 DiagonalMatrix D; D << LT; GM1 = D; LT1 = GM1; LT1 -= D; Print(LT1);
00102 UpperTriangularMatrix UT1 = GM1; UT1 -= D; Print(UT1);
00103 }
00104 {
00105 Tracer et1("Stage 4");
00106
00107 Matrix M(12,12); M = 0;
00108 M(1,1) = M(2,2) = M(4,4) = M(6,6) =
00109 M(7,7) = M(8,8) = M(10,10) = M(12,12) = -1;
00110 M(1,6) = M(1,12) = -5.601594;
00111 M(3,6) = M(3,12) = -0.000165;
00112 M(7,6) = M(7,12) = -0.008294;
00113 DiagonalMatrix D;
00114 SVD(M,D);
00115 SortDescending(D);
00116
00117 DiagonalMatrix DX(12);
00118 DX(1) = 8.0461;
00119 DX(2) = DX(3) = DX(4) = DX(5) = DX(6) = DX(7) = 1;
00120 DX(8) = 0.1243;
00121 DX(9) = DX(10) = DX(11) = DX(12) = 0;
00122 D -= DX; Clean(D,0.0001); Print(D);
00123 }
00124 #ifndef DONT_DO_NRIC
00125 {
00126 Tracer et1("Stage 5");
00127
00128 DiagonalMatrix D(10);
00129 D << 1 << 4 << 6 << 2 << 1 << 6 << 4 << 7 << 3 << 1;
00130 ColumnVector C(10);
00131 C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
00132 RowVector R(6);
00133 R << 2 << 3 << 5 << 7 << 11 << 13;
00134 nricMatrix M(10, 6);
00135 DCR( D.nric(), C.nric(), 10, R.nric(), 6, M.nric() );
00136 M -= D * C * R; Print(M);
00137
00138 D.ReSize(5);
00139 D << 1.25 << 4.75 << 9.5 << 1.25 << 3.75;
00140 C.ReSize(5);
00141 C << 1.5 << 7.5 << 4.25 << 0.0 << 7.25;
00142 R.ReSize(9);
00143 R << 2.5 << 3.25 << 5.5 << 7 << 11.25 << 13.5 << 0.0 << 1.5 << 3.5;
00144 Matrix MX = D * C * R;
00145 M.ReSize(MX);
00146 DCR( D.nric(), C.nric(), 5, R.nric(), 9, M.nric() );
00147 M -= MX; Print(M);
00148
00149
00150 nricMatrix A(3,4); nricMatrix B(4,5);
00151 A.Row(1) << 2 << 7 << 3 << 6;
00152 A.Row(2) << 6 << 2 << 5 << 9;
00153 A.Row(3) << 1 << 0 << 1 << 6;
00154 B.Row(1) << 2 << 8 << 4 << 5 << 3;
00155 B.Row(2) << 1 << 7 << 5 << 3 << 9;
00156 B.Row(3) << 7 << 8 << 2 << 1 << 6;
00157 B.Row(4) << 5 << 2 << 9 << 0 << 9;
00158 nricMatrix A1(1,2); nricMatrix B1;
00159 nricMatrix X(3,5); Matrix X1 = A * B;
00160 swap(A, A1); swap(B1, B);
00161 for (int i = 1; i <= 3; ++i) for (int j = 1; j <= 5; ++j)
00162 {
00163 X.nric()[i][j] = 0.0;
00164 for (int k = 1; k <= 4; ++k)
00165 X.nric()[i][j] += A1.nric()[i][k] * B1.nric()[k][j];
00166 }
00167 X1 -= X; Print(X1);
00168 }
00169 #endif
00170 {
00171 Tracer et1("Stage 6");
00172
00173 DiagonalMatrix test(5); test = 1;
00174 ColumnVector C(10);
00175 C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
00176 RowVector R(10);
00177 R << 2 << 3 << 5 << 7 << 11 << 13 << -3 << -4 << 2 << 4;
00178 test(1) = (R * C).AsScalar() - DotProduct(C, R);
00179 test(2) = C.SumSquare() - DotProduct(C, C);
00180 test(3) = 6.0 * (C.t() * R.t()).AsScalar() - DotProduct(2.0 * C, 3.0 * R);
00181 Matrix MC = C.AsMatrix(2,5), MR = R.AsMatrix(5,2);
00182 test(4) = DotProduct(MC, MR) - (R * C).AsScalar();
00183 UpperTriangularMatrix UT(5);
00184 UT << 3 << 5 << 2 << 1 << 7
00185 << 1 << 1 << 8 << 2
00186 << 7 << 0 << 1
00187 << 3 << 5
00188 << 6;
00189 LowerTriangularMatrix LT(5);
00190 LT << 5
00191 << 2 << 3
00192 << 1 << 0 << 7
00193 << 9 << 8 << 1 << 2
00194 << 0 << 2 << 1 << 9 << 2;
00195 test(5) = DotProduct(UT, LT) - Sum(SP(UT, LT));
00196 Print(test);
00197
00198 LowerTriangularMatrix LT1(5);
00199 LT1.Row(1) << 5;
00200 LT1.Row(2) << 2 << 3;
00201 LT1.Row(3) << 1 << 0 << 7;
00202 LT1.Row(4) << 9 << 8 << 1 << 2;
00203 LT1.Row(5) << 0 << 2 << 1 << 9 << 2;
00204 Matrix M = LT1 - LT; Print(M);
00205
00206 IdentityMatrix IM(5); IM *= 2;
00207 LinearEquationSolver LES1(IM);
00208 LowerTriangularMatrix LTX = LES1.i() * LT;
00209 M = LTX * 2 - LT; Print(M);
00210 DiagonalMatrix D = IM;
00211 LinearEquationSolver LES2(IM);
00212 LTX = LES2.i() * LT;
00213 M = LTX * 2 - LT; Print(M);
00214 UpperTriangularMatrix UTX = LES1.i() * UT;
00215 M = UTX * 2 - UT; Print(M);
00216 UTX = LES2.i() * UT;
00217 M = UTX * 2 - UT; Print(M);
00218 }
00219
00220 {
00221 Tracer et1("Stage 7");
00222
00223
00224 BandMatrix BM1(6,2,3);
00225 BM1.Row(1) << 3 << 8 << 4 << 1;
00226 BM1.Row(2) << 5 << 1 << 9 << 7 << 2;
00227 BM1.Row(3) << 1 << 0 << 6 << 3 << 1 << 3;
00228 BM1.Row(4) << 4 << 2 << 5 << 2 << 4;
00229 BM1.Row(5) << 3 << 3 << 9 << 1;
00230 BM1.Row(6) << 4 << 2 << 9;
00231 BandMatrix BM2(6,1,1);
00232 BM2.Row(1) << 2.5 << 7.5;
00233 BM2.Row(2) << 1.5 << 3.0 << 8.5;
00234 BM2.Row(3) << 6.0 << 6.5 << 7.0;
00235 BM2.Row(4) << 2.5 << 2.0 << 8.0;
00236 BM2.Row(5) << 0.5 << 4.5 << 3.5;
00237 BM2.Row(6) << 9.5 << 7.5;
00238 Matrix RM1 = BM1, RM2 = BM2;
00239 Matrix X;
00240 GenericMatrix GRM1 = RM1, GBM1 = BM1, GRM2 = RM2, GBM2 = BM2;
00241 Matrix Z(6,0); Z = 5; Print(Z);
00242 GRM1 |= Z; GBM1 |= Z; GRM2 &= Z.t(); GBM2 &= Z.t();
00243 X = GRM1 - BM1; Print(X); X = GBM1 - BM1; Print(X);
00244 X = GRM2 - BM2; Print(X); X = GBM2 - BM2; Print(X);
00245
00246 GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
00247 GRM1 *= GRM2; GBM1 *= GBM2;
00248 X = GRM1 - BM1 * BM2; Print(X);
00249 X = RM1 * RM2 - GBM1; Print(X);
00250
00251 GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
00252 GRM1 *= GBM2; GBM1 *= GRM2;
00253 X = GRM1 - BM1 * BM2; Print(X);
00254 X = RM1 * RM2 - GBM1; Print(X);
00255
00256 X = BM1.t(); BandMatrix BM1X = BM1.t();
00257 GRM1 = RM1; X -= GRM1.t(); Print(X); X = BM1X - BM1.t(); Print(X);
00258
00259
00260 IdentityMatrix IM(6); IM *= 2;
00261 GBM1 = BM1; GBM1 *= 4; GRM1 = RM1; GRM1 *= 4;
00262 DiagonalMatrix D = IM;
00263 LinearEquationSolver LES1(D);
00264 BandMatrix BX;
00265 BX = LES1.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00266 LinearEquationSolver LES2(IM);
00267 BX = LES2.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00268 BX = D.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00269 BX = IM.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00270 BX = IM.i(); BX *= GBM1; BX -= BM1 * 2; X = BX; Print(X);
00271
00272
00273 SymmetricBandMatrix SBM; SBM << SP(BM1, BM1.t());
00274 SBM << IM.i() * SBM;
00275 X = 2 * SBM - SP(RM1, RM1.t()); Print(X);
00276
00277
00278 D << 2.5 << 7.5 << 2 << 5 << 4.5 << 7.5;
00279 BX = D.i() * BM1; X = BX - D.i() * RM1;
00280 Clean(X,0.00000001); Print(X);
00281 BX = D.i(); BX *= BM1; X = BX - D.i() * RM1;
00282 Clean(X,0.00000001); Print(X);
00283 SBM << SP(BM1, BM1.t());
00284 BX = D.i() * SBM; X = BX - D.i() * SP(RM1, RM1.t());
00285 Clean(X,0.00000001); Print(X);
00286
00287
00288 BX = TestReturn(BM1); X = BX - BM1;
00289 if (BX.BandWidth() != BM1.BandWidth()) X = 5;
00290 Print(X);
00291 }
00292
00293
00294 }
00295
00296
00297