1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #define EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths 11 #include "main.h" 12 13 template<typename MatrixType, typename Index, typename Scalar> 14 typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type 15 block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) { 16 // check cwise-Functions: 17 VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1)); 18 VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1)); 19 20 VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1)); 21 VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1)); 22 23 return Scalar(0); 24 } 25 26 template<typename MatrixType, typename Index, typename Scalar> 27 typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type 28 block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) { 29 return Scalar(0); 30 } 31 32 33 template<typename MatrixType> void block(const MatrixType& m) 34 { 35 typedef typename MatrixType::Index Index; 36 typedef typename MatrixType::Scalar Scalar; 37 typedef typename MatrixType::RealScalar RealScalar; 38 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 39 typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; 40 typedef Matrix<Scalar, Dynamic, Dynamic> DynamicMatrixType; 41 typedef Matrix<Scalar, Dynamic, 1> DynamicVectorType; 42 43 Index rows = m.rows(); 44 Index cols = m.cols(); 45 46 MatrixType m1 = MatrixType::Random(rows, cols), 47 m1_copy = m1, 48 m2 = MatrixType::Random(rows, cols), 49 m3(rows, cols), 50 ones = MatrixType::Ones(rows, cols); 51 VectorType v1 = VectorType::Random(rows); 52 53 Scalar s1 = internal::random<Scalar>(); 54 55 Index r1 = internal::random<Index>(0,rows-1); 56 Index r2 = internal::random<Index>(r1,rows-1); 57 Index c1 = internal::random<Index>(0,cols-1); 58 Index c2 = internal::random<Index>(c1,cols-1); 59 60 block_real_only(m1, r1, r2, c1, c1, s1); 61 62 //check row() and col() 63 VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1)); 64 //check operator(), both constant and non-constant, on row() and col() 65 m1 = m1_copy; 66 m1.row(r1) += s1 * m1_copy.row(r2); 67 VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + s1 * m1_copy.row(r2)); 68 // check nested block xpr on lhs 69 m1.row(r1).row(0) += s1 * m1_copy.row(r2); 70 VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + Scalar(2) * s1 * m1_copy.row(r2)); 71 m1 = m1_copy; 72 m1.col(c1) += s1 * m1_copy.col(c2); 73 VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2)); 74 m1.col(c1).col(0) += s1 * m1_copy.col(c2); 75 VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2)); 76 77 78 //check block() 79 Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1); 80 81 RowVectorType br1(m1.block(r1,0,1,cols)); 82 VectorType bc1(m1.block(0,c1,rows,1)); 83 VERIFY_IS_EQUAL(b1, m1.block(r1,c1,1,1)); 84 VERIFY_IS_EQUAL(m1.row(r1), br1); 85 VERIFY_IS_EQUAL(m1.col(c1), bc1); 86 //check operator(), both constant and non-constant, on block() 87 m1.block(r1,c1,r2-r1+1,c2-c1+1) = s1 * m2.block(0, 0, r2-r1+1,c2-c1+1); 88 m1.block(r1,c1,r2-r1+1,c2-c1+1)(r2-r1,c2-c1) = m2.block(0, 0, r2-r1+1,c2-c1+1)(0,0); 89 90 enum { 91 BlockRows = 2, 92 BlockCols = 5 93 }; 94 if (rows>=5 && cols>=8) 95 { 96 // test fixed block() as lvalue 97 m1.template block<BlockRows,BlockCols>(1,1) *= s1; 98 // test operator() on fixed block() both as constant and non-constant 99 m1.template block<BlockRows,BlockCols>(1,1)(0, 3) = m1.template block<2,5>(1,1)(1,2); 100 // check that fixed block() and block() agree 101 Matrix<Scalar,Dynamic,Dynamic> b = m1.template block<BlockRows,BlockCols>(3,3); 102 VERIFY_IS_EQUAL(b, m1.block(3,3,BlockRows,BlockCols)); 103 104 // same tests with mixed fixed/dynamic size 105 m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols) *= s1; 106 m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols)(0,3) = m1.template block<2,5>(1,1)(1,2); 107 Matrix<Scalar,Dynamic,Dynamic> b2 = m1.template block<Dynamic,BlockCols>(3,3,2,5); 108 VERIFY_IS_EQUAL(b2, m1.block(3,3,BlockRows,BlockCols)); 109 } 110 111 if (rows>2) 112 { 113 // test sub vectors 114 VERIFY_IS_EQUAL(v1.template head<2>(), v1.block(0,0,2,1)); 115 VERIFY_IS_EQUAL(v1.template head<2>(), v1.head(2)); 116 VERIFY_IS_EQUAL(v1.template head<2>(), v1.segment(0,2)); 117 VERIFY_IS_EQUAL(v1.template head<2>(), v1.template segment<2>(0)); 118 Index i = rows-2; 119 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.block(i,0,2,1)); 120 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.tail(2)); 121 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.segment(i,2)); 122 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.template segment<2>(i)); 123 i = internal::random<Index>(0,rows-2); 124 VERIFY_IS_EQUAL(v1.segment(i,2), v1.template segment<2>(i)); 125 } 126 127 // stress some basic stuffs with block matrices 128 VERIFY(numext::real(ones.col(c1).sum()) == RealScalar(rows)); 129 VERIFY(numext::real(ones.row(r1).sum()) == RealScalar(cols)); 130 131 VERIFY(numext::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows)); 132 VERIFY(numext::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); 133 134 // now test some block-inside-of-block. 135 136 // expressions with direct access 137 VERIFY_IS_EQUAL( (m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2)) , (m1.block(r2,c2,rows-r2,cols-c2)) ); 138 VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , (m1.row(r1).segment(c1,c2-c1+1)) ); 139 VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , (m1.col(c1).segment(r1,r2-r1+1)) ); 140 VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0)) , (m1.row(r1).segment(c1,c2-c1+1)).transpose() ); 141 VERIFY_IS_EQUAL( (m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0)) , (m1.row(r1).segment(c1,c2-c1+1)).transpose() ); 142 143 // expressions without direct access 144 VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2)) , ((m1+m2).block(r2,c2,rows-r2,cols-c2)) ); 145 VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)) ); 146 VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , ((m1+m2).col(c1).segment(r1,r2-r1+1)) ); 147 VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() ); 148 VERIFY_IS_EQUAL( ((m1+m2).transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() ); 149 150 // evaluation into plain matrices from expressions with direct access (stress MapBase) 151 DynamicMatrixType dm; 152 DynamicVectorType dv; 153 dm.setZero(); 154 dm = m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2); 155 VERIFY_IS_EQUAL(dm, (m1.block(r2,c2,rows-r2,cols-c2))); 156 dm.setZero(); 157 dv.setZero(); 158 dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0).transpose(); 159 dv = m1.row(r1).segment(c1,c2-c1+1); 160 VERIFY_IS_EQUAL(dv, dm); 161 dm.setZero(); 162 dv.setZero(); 163 dm = m1.col(c1).segment(r1,r2-r1+1); 164 dv = m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0); 165 VERIFY_IS_EQUAL(dv, dm); 166 dm.setZero(); 167 dv.setZero(); 168 dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0); 169 dv = m1.row(r1).segment(c1,c2-c1+1); 170 VERIFY_IS_EQUAL(dv, dm); 171 dm.setZero(); 172 dv.setZero(); 173 dm = m1.row(r1).segment(c1,c2-c1+1).transpose(); 174 dv = m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0); 175 VERIFY_IS_EQUAL(dv, dm); 176 } 177 178 179 template<typename MatrixType> 180 void compare_using_data_and_stride(const MatrixType& m) 181 { 182 typedef typename MatrixType::Index Index; 183 Index rows = m.rows(); 184 Index cols = m.cols(); 185 Index size = m.size(); 186 Index innerStride = m.innerStride(); 187 Index outerStride = m.outerStride(); 188 Index rowStride = m.rowStride(); 189 Index colStride = m.colStride(); 190 const typename MatrixType::Scalar* data = m.data(); 191 192 for(int j=0;j<cols;++j) 193 for(int i=0;i<rows;++i) 194 VERIFY(m.coeff(i,j) == data[i*rowStride + j*colStride]); 195 196 if(!MatrixType::IsVectorAtCompileTime) 197 { 198 for(int j=0;j<cols;++j) 199 for(int i=0;i<rows;++i) 200 VERIFY(m.coeff(i,j) == data[(MatrixType::Flags&RowMajorBit) 201 ? i*outerStride + j*innerStride 202 : j*outerStride + i*innerStride]); 203 } 204 205 if(MatrixType::IsVectorAtCompileTime) 206 { 207 VERIFY(innerStride == int((&m.coeff(1))-(&m.coeff(0)))); 208 for (int i=0;i<size;++i) 209 VERIFY(m.coeff(i) == data[i*innerStride]); 210 } 211 } 212 213 template<typename MatrixType> 214 void data_and_stride(const MatrixType& m) 215 { 216 typedef typename MatrixType::Index Index; 217 Index rows = m.rows(); 218 Index cols = m.cols(); 219 220 Index r1 = internal::random<Index>(0,rows-1); 221 Index r2 = internal::random<Index>(r1,rows-1); 222 Index c1 = internal::random<Index>(0,cols-1); 223 Index c2 = internal::random<Index>(c1,cols-1); 224 225 MatrixType m1 = MatrixType::Random(rows, cols); 226 compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1)); 227 compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1)); 228 compare_using_data_and_stride(m1.row(r1)); 229 compare_using_data_and_stride(m1.col(c1)); 230 compare_using_data_and_stride(m1.row(r1).transpose()); 231 compare_using_data_and_stride(m1.col(c1).transpose()); 232 } 233 234 void test_block() 235 { 236 for(int i = 0; i < g_repeat; i++) { 237 CALL_SUBTEST_1( block(Matrix<float, 1, 1>()) ); 238 CALL_SUBTEST_2( block(Matrix4d()) ); 239 CALL_SUBTEST_3( block(MatrixXcf(3, 3)) ); 240 CALL_SUBTEST_4( block(MatrixXi(8, 12)) ); 241 CALL_SUBTEST_5( block(MatrixXcd(20, 20)) ); 242 CALL_SUBTEST_6( block(MatrixXf(20, 20)) ); 243 244 CALL_SUBTEST_8( block(Matrix<float,Dynamic,4>(3, 4)) ); 245 246 #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR 247 CALL_SUBTEST_6( data_and_stride(MatrixXf(internal::random(5,50), internal::random(5,50))) ); 248 CALL_SUBTEST_7( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(internal::random(5,50), internal::random(5,50))) ); 249 #endif 250 } 251 } 252
