$search
00001 00002 00003 00006 00007 00008 //#define WANT_STREAM 00009 #define WANT_MATH 00010 00011 00012 #include "include.h" 00013 00014 #include "newmatap.h" 00015 00016 #include "tmt.h" 00017 00018 #ifdef use_namespace 00019 using namespace NEWMAT; 00020 #endif 00021 00022 00023 /**************************** test program ******************************/ 00024 00025 00026 // slow sort program 00027 00028 static void SimpleSortDescending(Real* first, const int length) 00029 { 00030 int i = length; 00031 while (--i) 00032 { 00033 Real x = *first; Real* f = first; Real* g = f; 00034 int j = i; 00035 while (j--) if (x < *(++f)) { g = f; x = *g; } 00036 *g = *first; *first++ = x; 00037 } 00038 } 00039 00040 static void TestSort(int n) 00041 { 00042 // make some data 00043 RowVector X(n); 00044 int i; 00045 for (i = 1; i <= n; i++) 00046 X(i) = sin((Real)i) + 0.3 * cos(i/5.0) - 0.6 * sin(i/7.0) + 0.2 * sin(2.0 * i); 00047 RowVector X_Sorted = X; SimpleSortDescending(X_Sorted.Store(), n); 00048 RowVector A = X + X.Reverse(); SimpleSortDescending(A.Store(), n); 00049 00050 // test descending sort 00051 00052 RowVector Y = X; SortDescending(Y); Y -= X_Sorted; Print(Y); 00053 Y = X_Sorted; SortDescending(Y); Y -= X_Sorted; Print(Y); 00054 Y = X_Sorted.Reverse(); SortDescending(Y); Y -= X_Sorted; Print(Y); 00055 Y = X + X.Reverse(); SortDescending(Y); Y -= A; Print(Y); 00056 00057 // test ascending sort 00058 00059 Y = X; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y); 00060 Y = X_Sorted; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y); 00061 Y = X_Sorted.Reverse(); SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y); 00062 Y = X + X.Reverse(); SortAscending(Y); Y -= A.Reverse(); Print(Y); 00063 } 00064 00065 00066 void trymat6() 00067 { 00068 Tracer et("Sixth test of Matrix package"); 00069 Tracer::PrintTrace(); 00070 00071 int i,j; 00072 00073 00074 DiagonalMatrix D(6); 00075 UpperTriangularMatrix U(6); 00076 for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; } 00077 LowerTriangularMatrix L=(U*3.0).t(); 00078 SymmetricMatrix S(6); 00079 for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j; 00080 Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S; 00081 Matrix M(6,6); 00082 for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0; 00083 { 00084 Tracer et1("Stage 1"); 00085 Print(Matrix(MS+(-MS))); 00086 Print(Matrix((S+M)-(MS+M))); 00087 Print(Matrix((M+U)-(M+MU))); 00088 Print(Matrix((M+L)-(M+ML))); 00089 } 00090 { 00091 Tracer et1("Stage 2"); 00092 Print(Matrix((M+D)-(M+MD))); 00093 Print(Matrix((U+D)-(MU+MD))); 00094 Print(Matrix((D+L)-(ML+MD))); 00095 Print(Matrix((-U+D)+MU-MD)); 00096 Print(Matrix((-L+D)+ML-MD)); 00097 } 00098 { 00099 Tracer et1("Stage 3 - concatenate"); 00100 RowVector A(5); 00101 A << 1 << 2 << 3 << 4 << 5; 00102 RowVector B(5); 00103 B << 3 << 1 << 4 << 1 << 5; 00104 Matrix C(3,5); 00105 C << 2 << 3 << 5 << 7 << 11 00106 << 13 << 17 << 19 << 23 << 29 00107 << 31 << 37 << 41 << 43 << 47; 00108 Matrix X1 = A & B & C; 00109 Matrix X2 = (A.t() | B.t() | C.t()).t(); 00110 Matrix X3(5,5); 00111 X3.Row(1)=A; X3.Row(2)=B; X3.Rows(3,5)=C; 00112 Print(Matrix(X1-X2)); 00113 Print(Matrix(X1-X3)); 00114 LowerTriangularMatrix LT1; LT1 << (A & B & C); 00115 UpperTriangularMatrix UT1; UT1 << (A.t() | B.t() | C.t()); 00116 Print(LowerTriangularMatrix(LT1-UT1.t())); 00117 DiagonalMatrix D1; D1 << (A.t() | B.t() | C.t()); 00118 ColumnVector At = A.t(); 00119 ColumnVector Bt = B.t(); 00120 Matrix Ct = C.t(); 00121 LowerTriangularMatrix LT2; LT2 << (At | Bt | Ct); 00122 UpperTriangularMatrix UT2; UT2 << (At.t() & Bt.t() & Ct.t()); 00123 Matrix ABt = At | Bt; 00124 DiagonalMatrix D2; D2 << (ABt | Ct); 00125 Print(LowerTriangularMatrix(LT2-UT2.t())); 00126 Print(DiagonalMatrix(D1-D2)); 00127 Print(Matrix(LT1+UT2-D2-X1)); 00128 Matrix M1 = LT1 | UT2; Matrix M2 = UT1 & LT2; 00129 Print(Matrix(M1-M2.t())); 00130 M1 = UT2 | LT1; M2 = LT2 & UT1; 00131 Print(Matrix(M1-M2.t())); 00132 M1 = (LT1 | UT2) & (UT2 | LT1); 00133 M2 = (UT1 & LT2) | (LT2 & UT1); 00134 Print(Matrix(M1-M2.t())); 00135 SymmetricMatrix SM1; SM1 << (M1 + M1.t()); 00136 SymmetricMatrix SM2; SM2 << ((SM1 | M1) & (M1.t() | SM1)); 00137 Matrix M3(20,20); 00138 M3.SubMatrix(1,10,1,10) = SM1; 00139 M3.SubMatrix(1,10,11,20) = M1; 00140 M3.SubMatrix(11,20,1,10) = M2; 00141 M3.SubMatrix(11,20,11,20) = SM1; 00142 Print(Matrix(M3-SM2)); 00143 00144 SymmetricMatrix SM(15); SM = 0; SM.SymSubMatrix(1,10) = SM1; 00145 M3.ReSize(15,15); M3 = 0; M3.SubMatrix(1,10,1,10) = SM1; 00146 M3 -= SM; Print(M3); 00147 SM = 0; SM.SymSubMatrix(6,15) = SM1; 00148 M3.ReSize(15,15); M3 = 0; M3.SubMatrix(6,15,6,15) = SM1; 00149 M3 = M3.t() - SM; Print(M3); 00150 } 00151 { 00152 Tracer et1("Stage 4 - sort"); 00153 TestSort(1); TestSort(2); TestSort(3); TestSort(4); 00154 TestSort(15); TestSort(16); TestSort(17); TestSort(18); 00155 TestSort(99); TestSort(100); TestSort(101); 00156 } 00157 00158 00159 // cout << "\nEnd of sixth test\n"; 00160 } 00161 00162