Go to the documentation of this file.00001
00002
00003
00004 #include "include.h"
00005
00006 #include "newmat.h"
00007
00008 #include "tmt.h"
00009
00010 #ifdef use_namespace
00011 using namespace NEWMAT;
00012 #endif
00013
00014
00015
00016
00017 void trymat3()
00018 {
00019 Tracer et("Third test of Matrix package");
00020 Tracer::PrintTrace();
00021
00022 {
00023 Tracer et1("Stage 1");
00024 int i,j;
00025 SymmetricMatrix S(7);
00026 for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
00027 S=-S+2.0;
00028
00029 DiagonalMatrix D(7);
00030 for (i=1;i<=7;i++) D(i,i)=S(i,i);
00031
00032 Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0; M4=M4-4.0; Print(M4); }
00033 SymmetricMatrix S2=D; Matrix M2=S2; { M2=-D+M2; Print(M2); }
00034 UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
00035 LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
00036 M2=D; M2=M2-D; Print(M2);
00037 for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
00038 U2=L2.t(); D=D.t(); S=S.t();
00039 M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
00040 M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
00041 }
00042
00043 {
00044 Tracer et1("Stage 2");
00045 int i,j;
00046 DiagonalMatrix D(6);
00047 for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
00048 UpperTriangularMatrix U2(7); LowerTriangularMatrix L2(7);
00049 for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
00050 { U2=L2.t(); }
00051 DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
00052 Matrix M2(6,7);
00053 for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
00054 Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
00055 Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
00056 MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
00057 {
00058 D.ReSize(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
00059 LowerTriangularMatrix LD(1); LD=D;
00060 UpperTriangularMatrix UD(1); UD=D;
00061 M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
00062 M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
00063 M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
00064 M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
00065 }
00066 }
00067
00068 {
00069 Tracer et1("Stage 3");
00070
00071 DiagonalMatrix D(6);
00072 for (int i=1;i<=6;i++) D(i)=i*i;
00073 DiagonalMatrix E = D.i() * 4.0;
00074 DiagonalMatrix I(6); I = 1.0;
00075 E=D*E-I*4.0; Print(E);
00076 E = D.i() / 0.25; E=D*E-I*4.0; Print(E);
00077 }
00078 {
00079 Tracer et1("Stage 4");
00080 Matrix sigma(3,3); Matrix sigmaI(3,3);
00081 sigma = 0; sigma(1,1) = 1.0; sigma(2,2) = 1.0; sigma(3,3) = 1.0;
00082 sigmaI = sigma.i();
00083 sigmaI -= sigma; Clean(sigmaI, 0.000000001); Print(sigmaI);
00084 }
00085 {
00086 Tracer et1("Stage 5");
00087 Matrix X(5,5); DiagonalMatrix DM(5);
00088 for (int i=1; i<=5; i++) for (int j=1; j<=5; j++)
00089 X(i,j) = (23*i+59*j) % 43;
00090 DM << 1 << 8 << -7 << 2 << 3;
00091 Matrix Y = X.i() * DM; Y = X * Y - DM;
00092 Clean(Y, 0.000000001); Print(Y);
00093 }
00094 {
00095 Tracer et1("Stage 6");
00096 ColumnVector CV(10), RCV(10);
00097 CV << 2 << 7 << 1 << 6 << -3 << 1 << 8 << -4 << 0 << 17;
00098 RCV << 17 << 0 << -4 << 8 << 1 << -3 << 6 << 1 << 7 << 2;
00099 ColumnVector X = CV - RCV.Reverse(); Print(X);
00100 RowVector Y = CV.t() - RCV.t().Reverse(); Print(Y);
00101 DiagonalMatrix D = CV.AsDiagonal() - RCV.AsDiagonal().Reverse();
00102 Print(D);
00103 X = CV & CV.Rows(1,9).Reverse();
00104 ColumnVector Z(19);
00105 Z.Rows(1,10) = RCV.Reverse(); Z.Rows(11,19) = RCV.Rows(2,10);
00106 X -= Z; Print(X); Z -= Z.Reverse(); Print(Z);
00107 Matrix A(3,3); A << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9;
00108 Matrix B(3,3); B << 9 << 8 << 7 << 6 << 5 << 4 << 3 << 2 << 1;
00109 Matrix Diff = A - B.Reverse(); Print(Diff);
00110 Diff = (-A).Reverse() + B; Print(Diff);
00111 UpperTriangularMatrix U;
00112 U << A.Reverse(); Diff = U; U << B; Diff -= U; Print(Diff);
00113 U << (-A).Reverse(); Diff = U; U << B; Diff += U; Print(Diff);
00114 }
00115 {
00116 Tracer et1("Stage 7");
00117 ColumnVector XX(4);
00118 Matrix A(3,3);
00119 A = 0;
00120 CroutMatrix B1 = A;
00121 XX(1) = B1.IsSingular() ? 0 : 1;
00122 A << 1 << 3 << 6
00123 << 7 << 11 << 13
00124 << 2 << 4 << 1;
00125 CroutMatrix B2(A);
00126 XX(2) = B2.IsSingular() ? 1 : 0;
00127 BandMatrix C(3,1,1); C.Inject(A);
00128 BandLUMatrix B3(C);
00129 XX(3) = B3.IsSingular() ? 1 : 0;
00130 C = 0;
00131 BandLUMatrix B4(C);
00132 XX(4) = B4.IsSingular() ? 0 : 1;
00133 Print(XX);
00134 }
00135 {
00136 Tracer et1("Stage 8");
00137 Matrix A(3,3); Matrix Z(3,3); ColumnVector X(6);
00138 A << 1 << 3 << 6
00139 << 7 << 11 << 13
00140 << 2 << 4 << 1;
00141 Z = 0;
00142 Matrix B = (A | Z) & (Z | A);
00143 X = 0.0;
00144 X = B.i() * X;
00145 Print(X);
00146
00147 Matrix Y(3,3);
00148 Y << 0.0 << 1.0 << 1.0
00149 << 5.0 << 0.0 << 5.0
00150 << 3.0 << 3.0 << 0.0;
00151 Matrix YY = Y & Y;
00152 YY = B.i() * YY;
00153 Matrix Y1 = A.i() * Y;
00154 YY -= Y1 & Y1; Clean(YY, 0.000000001); Print(YY);
00155 Y1 = A * Y1 - Y; Clean(Y1, 0.000000001); Print(Y1);
00156 }
00157
00158
00159 }
00160