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