00001
00002
00003
00006
00007
00008
00009 #include "include.h"
00010 #include "newmatap.h"
00011
00012 #include "tmt.h"
00013
00014 #ifdef use_namespace
00015 using namespace NEWMAT;
00016 #endif
00017
00018
00019
00020
00021
00022 void trymat9()
00023 {
00024 Tracer et("Ninth test of Matrix package");
00025 Tracer::PrintTrace();
00026
00027
00028 int i; int j;
00029 Matrix A(7,7); Matrix X(7,3);
00030 for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
00031 for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
00032 Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
00033 {
00034 Tracer et1("Stage 1");
00035 Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
00036 Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
00037 Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
00038 }
00039 for (i=1;i<=7;i++) D(i,i)=i*i+1;
00040 DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
00041 {
00042 Tracer et1("Stage 2");
00043 Matrix DXE = D.i() * X * E;
00044 DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE);
00045 E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
00046 }
00047 DiagonalMatrix F=D;
00048 {
00049 Tracer et1("Stage 3");
00050 F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
00051 F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
00052 }
00053 {
00054 Tracer et1("Stage 4");
00055 F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
00056 }
00057 {
00058 Tracer et1("Stage 5");
00059
00060 ColumnVector A(18), B(18);
00061 Matrix X(3,3);
00062 X << 3 << 5 << 7 << 5 << 8 << 2 << 7 << 2 << 9;
00063 SymmetricMatrix S; S << X;
00064 B(1) = S == X; A(1) = true;
00065 B(2) = S == (X+1); A(2) = false;
00066 B(3) = (S+2) == (X+2); A(3) = true;
00067 Matrix Y = X;
00068 B(4) = X == Y; A(4) = true;
00069 B(5) = (X*2) == (Y*2); A(5) = true;
00070 Y(3,3) = 10;
00071 B(6) = X == Y; A(6) = false;
00072 B(7) = (X*2) == (Y*2); A(7) = false;
00073 B(8) = S == Y; A(8) = false;
00074 B(9) = S == S; A(9) = true;
00075 Matrix Z = X.SubMatrix(1,2,2,3);
00076 B(10) = X == Z; A(10) = false;
00077 GenericMatrix GS = S;
00078 GenericMatrix GX = X;
00079 GenericMatrix GY = Y;
00080 B(11) = GS == GX; A(11) = true;
00081 B(12) = GS == GY; A(12) = false;
00082 CroutMatrix CS = S;
00083 CroutMatrix CX = X;
00084 CroutMatrix CY = Y;
00085 B(13) = CS == CX; A(13) = true;
00086 B(14) = CS == CY; A(14) = false;
00087 B(15) = X == CX; A(15) = false;
00088 B(16) = X == A; A(16) = false;
00089 B(17) = X == (X | X); A(17) = false;
00090 B(18) = CX == X; A(18) = false;
00091 A = A - B; Print(A);
00092 }
00093 {
00094 Tracer et1("Stage 6");
00095
00096 ColumnVector A(22), B(22);
00097 BandMatrix X(6,2,1);
00098 X(1,1)=23; X(1,2)=21;
00099 X(2,1)=12; X(2,2)=17; X(2,3)=45;
00100 X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
00101 X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
00102 X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
00103 X(6,4)=49; X(6,5)=41; X(6,6)=17;
00104 SymmetricBandMatrix S1(6,2); S1.Inject(X);
00105 BandMatrix U(6,2,3); U = 0.0; U.Inject(X);
00106 B(1) = U == X; A(1) = true;
00107 B(2) = U == (X*3); A(2) = false;
00108 B(3) = (U*5) == (X*5); A(3) = true;
00109 Matrix Y = X;
00110 B(4) = X == Y; A(4) = true;
00111 B(5) = (X*2) == (Y*2); A(5) = true;
00112 Y(6,6) = 10;
00113 B(6) = X == Y; A(6) = false;
00114 B(7) = (X*2) == (Y*2); A(7) = false;
00115 B(8) = U == Y; A(8) = false;
00116 B(9) = U == U; A(9) = true;
00117 Matrix Z = X.SubMatrix(1,2,2,3);
00118 B(10) = X == Z; A(10) = false;
00119 GenericMatrix GU = U;
00120 GenericMatrix GX = X;
00121 GenericMatrix GY = Y;
00122 B(11) = GU == GX; A(11) = true;
00123 B(12) = GU == GY; A(12) = false;
00124 X = X + X.t(); U = U + U.t();
00125 SymmetricBandMatrix S(6,2); S.Inject(X);
00126 Matrix D = S-X; Print(D);
00127 BandLUMatrix BS = S;
00128 BandLUMatrix BX = X;
00129 BandLUMatrix BU = U;
00130 CroutMatrix CX = X;
00131 B(13) = BS == BX; A(13) = true;
00132 B(14) = BX == BU; A(14) = false;
00133 B(15) = X == BX; A(15) = false;
00134 B(16) = X != BX; A(16) = true;
00135 B(17) = BX != BS; A(17) = false;
00136 B(18) = (2*X) != (X*2);A(18) = false;
00137 B(19) = (X*2) != (X+2);A(19) = true;
00138 B(20) = BX == CX; A(20) = false;
00139 B(21) = CX == BX; A(21) = false;
00140 B(22) = BX == X; A(22) = false;
00141 A = A - B; Print(A);
00142 DiagonalMatrix I(6); I=1.0;
00143 D = BS.i() * X - I; Clean(D,0.00000001); Print(D);
00144 D = BX.i() * X - I; Clean(D,0.00000001); Print(D);
00145 D = BU.i() * X - I; Clean(D,0.00000001); Print(D);
00146
00147
00148 SymmetricBandMatrix X1(6,2);
00149 X1.Row(1) << 23;
00150 X1.Row(2) << 12 << 17;
00151 X1.Row(3) << 35 << 19 << 24;
00152 X1.Row(4) << 17 << 11 << 19;
00153 X1.Row(5) << 10 << 44 << 23;
00154 X1.Row(6) << 49 << 41 << 17;
00155 Matrix M = X1 - S1; Print(M);
00156
00157
00158 SymmetricBandMatrix X2(20,3); X2 = 0.0;
00159 X2.SubMatrix(2,7,2,7) = X1; X2.SymSubMatrix(11,16) = 2 * X1;
00160 Matrix MX1 = X1;
00161 Matrix MX2(20,20); MX2 = 0;
00162 MX2.SymSubMatrix(2,7) = MX1; MX2.SubMatrix(11,16,11,16) = MX1 * 2;
00163 MX2 -= X2; Print(MX2);
00164
00165 BandMatrix X4(20,3,3); X4 = 0.0;
00166 X4.SubMatrix(2,7,3,8) = X1; X4.SubMatrix(11,16,10,15) = 2 * X1;
00167 MX1 = X1;
00168 Matrix MX4(20,20); MX4 = 0;
00169 MX4.SubMatrix(2,7,3,8) = MX1; MX4.SubMatrix(11,16,10,15) = MX1 * 2;
00170 MX4 -= X4; Print(MX4);
00171
00172 MX1 = X1.i() * X1 - IdentityMatrix(6);
00173 Clean(MX1,0.00000001); Print(MX1);
00174
00175 }
00176
00177 {
00178 Tracer et1("Stage 7");
00179
00180 ColumnVector A(12), B(12);
00181 BandMatrix X(6,2,1);
00182 X(1,1)=23; X(1,2)=21;
00183 X(2,1)=12; X(2,2)=17; X(2,3)=45;
00184 X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
00185 X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
00186 X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
00187 X(6,4)=49; X(6,5)=41; X(6,6)=17;
00188 Matrix Y = X;
00189 LinearEquationSolver LX = X;
00190 LinearEquationSolver LY = Y;
00191 CroutMatrix CX = X;
00192 CroutMatrix CY = Y;
00193 BandLUMatrix BX = X;
00194 B(1) = LX == CX; A(1) = false;
00195 B(2) = LY == CY; A(2) = true;
00196 B(3) = X == Y; A(3) = true;
00197 B(4) = BX == LX; A(4) = true;
00198 B(5) = CX == CY; A(5) = true;
00199 B(6) = LX == LY; A(6) = false;
00200 B(7) = BX == BX; A(7) = true;
00201 B(8) = CX == CX; A(8) = true;
00202 B(9) = LX == LX; A(9) = true;
00203 B(10) = LY == LY; A(10) = true;
00204 CroutMatrix CX1 = X.SubMatrix(1,4,1,4);
00205 B(11) = CX == CX1; A(11) = false;
00206 BandLUMatrix BX1 = X.SymSubMatrix(1,4);
00207 B(12) = BX == BX1; A(12) = false;
00208 A = A - B; Print(A);
00209 DiagonalMatrix I(6); I=1.0; Matrix D;
00210 D = LX.i() * X - I; Clean(D,0.00000001); Print(D);
00211 D = LY.i() * X - I; Clean(D,0.00000001); Print(D);
00212 I.ReSize(4); I = 1;
00213 D = CX1.i() * X.SymSubMatrix(1,4) - I; Clean(D,0.00000001); Print(D);
00214 D = BX1.i() * X.SubMatrix(1,4,1,4) - I; Clean(D,0.00000001); Print(D);
00215 }
00216
00217 {
00218 Tracer et1("Stage 8");
00219
00220 MultWithCarry MWC;
00221 SymmetricBandMatrix SBM(50, 10);
00222 for (int i = 1; i <= 50; ++i) for (int j = 1; j <= i; ++j)
00223 if (i - j <= 10) SBM(i, j) = MWC.Next();
00224 CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
00225 CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
00226 CM1.release(); BM1.release();
00227 CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
00228 Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
00229 Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
00230 X *= SBM; X -= IdentityMatrix(50); Clean(X,0.00000001); Print(X);
00231 LogAndSign x = log_determinant(SBM);
00232 LogAndSign y = log_determinant(CM2);
00233 LogAndSign z = log_determinant(BM2);
00234 RowVector D(4);
00235 D(1) = y.value() - x.value();
00236 D(2) = z.value() - x.value();
00237 D(3) = y.sign() - x.sign();
00238 D(4) = z.sign() - x.sign();
00239 Clean(D,0.00000001); Print(D);
00240 }
00241
00242 {
00243 Tracer et1("Stage 9");
00244
00245 MultWithCarry MWC;
00246 SymmetricBandMatrix SBM(51, 10);
00247 for (int i = 1; i <= 51; ++i) for (int j = 1; j <= i; ++j)
00248 if (i - j <= 10) SBM(i, j) = MWC.Next();
00249 CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
00250 CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
00251 CM1.release(); BM1.release();
00252 CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
00253 Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
00254 Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
00255 X *= SBM; X -= IdentityMatrix(51); Clean(X,0.00000001); Print(X);
00256 LogAndSign x = log_determinant(SBM);
00257 LogAndSign y = log_determinant(CM2);
00258 LogAndSign z = log_determinant(BM2);
00259 RowVector D(4);
00260 D(1) = y.value() - x.value();
00261 D(2) = z.value() - x.value();
00262 D(3) = y.sign() - x.sign();
00263 D(4) = z.sign() - x.sign();
00264 Clean(D,0.00000001); Print(D);
00265 }
00266
00267
00268
00269 }
00270
00271