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