00001
00002
00003
00006
00007
00008
00009
00010
00011 #include "include.h"
00012
00013 #include "newmatap.h"
00014
00015 #include "tmt.h"
00016
00017 #ifdef use_namespace
00018 using namespace NEWMAT;
00019 #endif
00020
00021
00022
00023
00024
00025 static void process(const GeneralMatrix& A,
00026 const ColumnVector& X1, const ColumnVector& X2)
00027 {
00028 Matrix B = A;
00029 LinearEquationSolver L(A);
00030 Matrix Y(4,2);
00031 Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
00032 Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
00033 Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
00034 }
00035
00036
00037
00038 void trymata()
00039 {
00040
00041 Tracer et("Tenth test of Matrix package");
00042 Tracer::PrintTrace();
00043 int i; int j;
00044 UpperTriangularMatrix U(8);
00045 for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
00046 Matrix X(8,6);
00047 for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
00048 Matrix Y = U.i()*X; Matrix MU=U;
00049 Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
00050 Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
00051 UpperTriangularMatrix UX(8);
00052 for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
00053 UX(4,4)=0; UX(4,5)=0;
00054 UpperTriangularMatrix UY = U.i() * UX;
00055 { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
00056 LowerTriangularMatrix LY = U.t().i() * UX.t();
00057 { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
00058 DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
00059 { X=D.i()*MU; }
00060 { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
00061 { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
00062
00063
00064
00065
00066 U.ReSize(8);
00067 for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
00068 MU = U;
00069 MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
00070 MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
00071
00072
00073 {
00074 ColumnVector X1(4), X2(4);
00075 X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
00076 X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
00077
00078
00079 Matrix A(4,4);
00080 A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
00081 A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
00082 A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
00083 A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
00084 process(A,X1,X2);
00085
00086 BandMatrix B(4,1,1); B.Inject(A);
00087 process(B,X1,X2);
00088
00089 UpperTriangularMatrix U(4);
00090 U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
00091 U(2,2)=8; U(2,3)=7; U(2,4)=6;
00092 U(3,3)=2; U(3,4)=7;
00093 U(4,4)=1;
00094 process(U,X1,X2);
00095
00096
00097 UpperTriangularMatrix U1(4);
00098 U1.Row(1) << 1 << 2 << 3 << 4;
00099 U1.Row(2) << 8 << 7 << 6;
00100 U1.Row(3) << 2 << 7;
00101 U1.Row(4) << 1;
00102
00103 U1 -= U;
00104
00105 Print(U1);
00106
00107 LowerTriangularMatrix L = U.t();
00108 process(L,X1,X2);
00109 }
00110
00111
00112
00113 {
00114 Matrix M(4,4);
00115
00116 M << 8.613057e+00 << 8.693985e+00 << -2.322050e-01 << 0.000000e+00
00117 << 8.693985e+00 << 8.793868e+00 << -2.346310e-01 << 0.000000e+00
00118 << -2.322050e-01 << -2.346310e-01 << 6.264000e-03 << 0.000000e+00
00119 << 0.000000e+00 << 0.000000e+00 << 0.000000e+00 << 3.282806e+03 ;
00120 Matrix MI = M.i();
00121 DiagonalMatrix I(4); I = 1;
00122 Matrix Diff = MI * M - I;
00123 Clean(Diff,0.00000001); Print(Diff);
00124
00125 SymmetricMatrix SM; SM << M;
00126 LowerTriangularMatrix LT = Cholesky(SM).i();
00127 MI = LT.t() * LT; Diff = MI * M - I;
00128 Clean(Diff,0.00000001); Print(Diff);
00129 }
00130
00131
00132 }
00133
00134