$search
00001 00002 00003 00006 00007 00008 #define WANT_STREAM 00009 00010 #define WANT_MATH 00011 00012 #include "newmat.h" 00013 00014 #include "tmt.h" 00015 00016 #ifdef use_namespace 00017 using namespace NEWMAT; 00018 #endif 00019 00020 00021 00022 // test maxima and minima 00023 00024 Real TestMax(const GenericMatrix& GM) 00025 { 00026 Matrix M = GM; 00027 ColumnVector CV = GM.AsColumn(); 00028 00029 int c, i, j, k; int n = CV.Nrows(), nr = M.Nrows(), nc = M.Ncols(); 00030 Real x1, x2, x3; 00031 00032 MatrixBandWidth mbw = GM.BandWidth(); 00033 int u = mbw.Upper(); int l = mbw.Lower(); 00034 bool IsBandMatrix = u > 0 && l > 0 && !(u == 0 && l == 0); 00035 00036 x1 = GM.MaximumAbsoluteValue(); 00037 x2 = GM.MaximumAbsoluteValue1(i); 00038 if (fabs(CV(i)) != x2) return 1.1; 00039 x3 = GM.MaximumAbsoluteValue2(i,j); 00040 if (fabs(M(i,j)) != x3) return 1.2; 00041 if (x3 != x1) return 1.3; 00042 if (x2 != x1) return 1.4; 00043 c = 0; 00044 for (k = 1; k <= n; k++) 00045 { 00046 Real cvk = fabs(CV(k)); 00047 if (x1 <= cvk) { if (x1 < cvk) return 1.5; else c++; } 00048 } 00049 if (c == 0) return 1.6; 00050 00051 00052 x1 = GM.MinimumAbsoluteValue(); 00053 x2 = GM.MinimumAbsoluteValue1(i); 00054 if (fabs(CV(i)) != x2) return 2.1; 00055 x3 = GM.MinimumAbsoluteValue2(i,j); 00056 if (fabs(M(i,j)) != x3) return 2.2; 00057 if (x3 != x1) return 2.3; 00058 if (x2 != CV.MinimumAbsoluteValue()) return 2.4; 00059 c = 0; 00060 if (IsBandMatrix) 00061 { 00062 for (i = 1; i <= nr; i++) for (j = 1; j <= nc; j++) 00063 if (i - j <= l && j - i <= u) 00064 { 00065 Real mij = fabs(M(i,j)); 00066 if (x1 >= mij) { if (x1 > mij) return 2.51; else c++; } 00067 } 00068 } 00069 else 00070 { 00071 for (k = 1; k <= n; k++) 00072 { 00073 Real cvk = fabs(CV(k)); 00074 if (x1 >= cvk) { if (x1 > cvk) return 2.52; else c++; } 00075 } 00076 } 00077 if (c == 0) return 2.6; 00078 00079 x1 = GM.Maximum(); 00080 x2 = GM.Maximum1(i); 00081 if (CV(i) != x2) return 3.1; 00082 x3 = GM.Maximum2(i,j); 00083 if (M(i,j) != x3) return 3.2; 00084 if (x3 != x1) return 3.3; 00085 if (x2 != CV.Maximum()) return 3.4; 00086 c = 0; 00087 if (IsBandMatrix) 00088 { 00089 for (i = 1; i <= nr; i++) for (j = 1; j <= nc; j++) 00090 if (i - j <= l && j - i <= u) 00091 { 00092 Real mij = M(i,j); 00093 if (x1 <= mij) { if (x1 < mij) return 3.51; else c++; } 00094 } 00095 } 00096 else 00097 { 00098 for (k = 1; k <= n; k++) 00099 { 00100 Real cvk = CV(k); 00101 if (x1 <= cvk) { if (x1 < cvk) return 3.52; else c++; } 00102 } 00103 } 00104 if (c == 0) return 3.6; 00105 00106 x1 = GM.Minimum(); 00107 x2 = GM.Minimum1(i); 00108 if (CV(i) != x2) return 4.1; 00109 x3 = GM.Minimum2(i,j); 00110 if (M(i,j) != x3) return 4.2; 00111 if (x3 != x1) return 4.3; 00112 if (x2 != CV.Minimum()) return 4.4; 00113 c = 0; 00114 if (IsBandMatrix) 00115 { 00116 for (i = 1; i <= nr; i++) for (j = 1; j <= nc; j++) 00117 if (i - j <= l && j - i <= u) 00118 { 00119 Real mij = M(i,j); 00120 if (x1 >= mij) { if (x1 > mij) return 4.51; else c++; } 00121 } 00122 } 00123 else 00124 { 00125 for (k = 1; k <= n; k++) 00126 { 00127 Real cvk = CV(k); 00128 if (x1 >= cvk) { if (x1 > cvk) return 4.52; else c++; } 00129 } 00130 } 00131 if (c == 0) return 4.6; 00132 00133 return 0; 00134 00135 } 00136 00137 00138 void trymatl() 00139 { 00140 Tracer et("Twenty first test of Matrix package"); 00141 Tracer::PrintTrace(); 00142 00143 Matrix M(4, 4); 00144 M << 1.5 << 1.0 << -4.0 << 4.5 00145 << 3.5 << 7.0 << 2.0 << -1.0 00146 << -1.5 << 7.5 << -6.0 << 5.0 00147 << 0.5 << -8.0 << 5.5 << 4.0; 00148 UpperTriangularMatrix UM; UM << M; 00149 LowerTriangularMatrix LM; LM << -M; 00150 SymmetricMatrix SM; SM << (UM + UM.t()); 00151 BandMatrix BM(4, 1, 2); BM.Inject(M); 00152 DiagonalMatrix DM; DM << M; 00153 SymmetricBandMatrix SBM(4,1); SBM.Inject(M); 00154 RowVector Test(28); int k = 0; 00155 00156 // tests for different shapes 00157 Test(++k) = TestMax(GenericMatrix(M)); 00158 Test(++k) = TestMax(GenericMatrix(UM)); 00159 Test(++k) = TestMax(GenericMatrix(LM)); 00160 Test(++k) = TestMax(GenericMatrix(SM)); 00161 Test(++k) = TestMax(GenericMatrix(DM)); 00162 Test(++k) = TestMax(GenericMatrix(BM)); 00163 Test(++k) = TestMax(GenericMatrix(SBM)); 00164 00165 // tests for constant matrices - don't put non-zeros in corners of BM 00166 Matrix MX(4,4); 00167 MX = 1; Test(++k) = TestMax(GenericMatrix(MX)); 00168 BM.Inject(MX); Test(++k) = TestMax(GenericMatrix(BM)); 00169 MX = 0; Test(++k) = TestMax(GenericMatrix(MX)); 00170 BM.Inject(MX); Test(++k) = TestMax(GenericMatrix(BM)); 00171 MX = -1; Test(++k) = TestMax(GenericMatrix(MX)); 00172 BM.Inject(MX); Test(++k) = TestMax(GenericMatrix(BM)); 00173 00174 // test for non-square 00175 MX = M | (2 * M).t(); Test(++k) = TestMax(GenericMatrix(MX)); 00176 00177 // test on row and column vector 00178 RowVector RV(6); 00179 RV << 1 << 3 << -5 << 2 << -4 << 3; 00180 Test(++k) = TestMax(GenericMatrix(RV)); 00181 Test(++k) = TestMax(GenericMatrix(RV.t())); 00182 00183 // test for function form 00184 Test(++k) = (MaximumAbsoluteValue(RV) - 5); 00185 Test(++k) = (MinimumAbsoluteValue(RV) - 1); 00186 Test(++k) = (Maximum(RV) - 3); 00187 Test(++k) = (Minimum(RV) + 5); 00188 Test(++k) = (MaximumAbsoluteValue(-RV) - 5); 00189 Test(++k) = (MinimumAbsoluteValue(-RV) - 1); 00190 Test(++k) = (Maximum(-RV) - 5); 00191 Test(++k) = (Minimum(-RV) + 3); 00192 00193 // test formulae 00194 RowVector RV2(6); 00195 RV2 << 2 << 8 << 1 << 9 << 3 << -1; 00196 Test(++k) = (MaximumAbsoluteValue(RV+RV2) - 11); 00197 Test(++k) = (MinimumAbsoluteValue(RV+RV2) - 1); 00198 Test(++k) = (Maximum(RV+RV2) - 11); 00199 Test(++k) = (Minimum(RV+RV2) + 4); 00200 00201 00202 Print(Test); 00203 } 00204 00205