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> void block(const MatrixType& m) 14 { 15 typedef typename MatrixType::Index Index; 16 typedef typename MatrixType::Scalar Scalar; 17 typedef typename MatrixType::RealScalar RealScalar; 18 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 19 typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; 20 typedef Matrix<Scalar, Dynamic, Dynamic> DynamicMatrixType; 21 typedef Matrix<Scalar, Dynamic, 1> DynamicVectorType; 22 23 Index rows = m.rows(); 24 Index cols = m.cols(); 25 26 MatrixType m1 = MatrixType::Random(rows, cols), 27 m1_copy = m1, 28 m2 = MatrixType::Random(rows, cols), 29 m3(rows, cols), 30 ones = MatrixType::Ones(rows, cols); 31 VectorType v1 = VectorType::Random(rows); 32 33 Scalar s1 = internal::random<Scalar>(); 34 35 Index r1 = internal::random<Index>(0,rows-1); 36 Index r2 = internal::random<Index>(r1,rows-1); 37 Index c1 = internal::random<Index>(0,cols-1); 38 Index c2 = internal::random<Index>(c1,cols-1); 39 40 //check row() and col() 41 VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1)); 42 //check operator(), both constant and non-constant, on row() and col() 43 m1 = m1_copy; 44 m1.row(r1) += s1 * m1_copy.row(r2); 45 VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + s1 * m1_copy.row(r2)); 46 // check nested block xpr on lhs 47 m1.row(r1).row(0) += s1 * m1_copy.row(r2); 48 VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + Scalar(2) * s1 * m1_copy.row(r2)); 49 m1 = m1_copy; 50 m1.col(c1) += s1 * m1_copy.col(c2); 51 VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2)); 52 m1.col(c1).col(0) += s1 * m1_copy.col(c2); 53 VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2)); 54 55 //check block() 56 Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1); 57 58 RowVectorType br1(m1.block(r1,0,1,cols)); 59 VectorType bc1(m1.block(0,c1,rows,1)); 60 VERIFY_IS_EQUAL(b1, m1.block(r1,c1,1,1)); 61 VERIFY_IS_EQUAL(m1.row(r1), br1); 62 VERIFY_IS_EQUAL(m1.col(c1), bc1); 63 //check operator(), both constant and non-constant, on block() 64 m1.block(r1,c1,r2-r1+1,c2-c1+1) = s1 * m2.block(0, 0, r2-r1+1,c2-c1+1); 65 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); 66 67 enum { 68 BlockRows = 2, 69 BlockCols = 5 70 }; 71 if (rows>=5 && cols>=8) 72 { 73 // test fixed block() as lvalue 74 m1.template block<BlockRows,BlockCols>(1,1) *= s1; 75 // test operator() on fixed block() both as constant and non-constant 76 m1.template block<BlockRows,BlockCols>(1,1)(0, 3) = m1.template block<2,5>(1,1)(1,2); 77 // check that fixed block() and block() agree 78 Matrix<Scalar,Dynamic,Dynamic> b = m1.template block<BlockRows,BlockCols>(3,3); 79 VERIFY_IS_EQUAL(b, m1.block(3,3,BlockRows,BlockCols)); 80 } 81 82 if (rows>2) 83 { 84 // test sub vectors 85 VERIFY_IS_EQUAL(v1.template head<2>(), v1.block(0,0,2,1)); 86 VERIFY_IS_EQUAL(v1.template head<2>(), v1.head(2)); 87 VERIFY_IS_EQUAL(v1.template head<2>(), v1.segment(0,2)); 88 VERIFY_IS_EQUAL(v1.template head<2>(), v1.template segment<2>(0)); 89 Index i = rows-2; 90 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.block(i,0,2,1)); 91 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.tail(2)); 92 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.segment(i,2)); 93 VERIFY_IS_EQUAL(v1.template tail<2>(), v1.template segment<2>(i)); 94 i = internal::random<Index>(0,rows-2); 95 VERIFY_IS_EQUAL(v1.segment(i,2), v1.template segment<2>(i)); 96 } 97 98 // stress some basic stuffs with block matrices 99 VERIFY(internal::real(ones.col(c1).sum()) == RealScalar(rows)); 100 VERIFY(internal::real(ones.row(r1).sum()) == RealScalar(cols)); 101 102 VERIFY(internal::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows)); 103 VERIFY(internal::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); 104 105 // now test some block-inside-of-block. 106 107 // expressions with direct access 108 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)) ); 109 VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , (m1.row(r1).segment(c1,c2-c1+1)) ); 110 VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , (m1.col(c1).segment(r1,r2-r1+1)) ); 111 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() ); 112 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() ); 113 114 // expressions without direct access 115 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)) ); 116 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)) ); 117 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)) ); 118 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() ); 119 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() ); 120 121 // evaluation into plain matrices from expressions with direct access (stress MapBase) 122 DynamicMatrixType dm; 123 DynamicVectorType dv; 124 dm.setZero(); 125 dm = m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2); 126 VERIFY_IS_EQUAL(dm, (m1.block(r2,c2,rows-r2,cols-c2))); 127 dm.setZero(); 128 dv.setZero(); 129 dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0).transpose(); 130 dv = m1.row(r1).segment(c1,c2-c1+1); 131 VERIFY_IS_EQUAL(dv, dm); 132 dm.setZero(); 133 dv.setZero(); 134 dm = m1.col(c1).segment(r1,r2-r1+1); 135 dv = m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0); 136 VERIFY_IS_EQUAL(dv, dm); 137 dm.setZero(); 138 dv.setZero(); 139 dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0); 140 dv = m1.row(r1).segment(c1,c2-c1+1); 141 VERIFY_IS_EQUAL(dv, dm); 142 dm.setZero(); 143 dv.setZero(); 144 dm = m1.row(r1).segment(c1,c2-c1+1).transpose(); 145 dv = m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0); 146 VERIFY_IS_EQUAL(dv, dm); 147 } 148 149 150 template<typename MatrixType> 151 void compare_using_data_and_stride(const MatrixType& m) 152 { 153 typedef typename MatrixType::Index Index; 154 Index rows = m.rows(); 155 Index cols = m.cols(); 156 Index size = m.size(); 157 Index innerStride = m.innerStride(); 158 Index outerStride = m.outerStride(); 159 Index rowStride = m.rowStride(); 160 Index colStride = m.colStride(); 161 const typename MatrixType::Scalar* data = m.data(); 162 163 for(int j=0;j<cols;++j) 164 for(int i=0;i<rows;++i) 165 VERIFY(m.coeff(i,j) == data[i*rowStride + j*colStride]); 166 167 if(!MatrixType::IsVectorAtCompileTime) 168 { 169 for(int j=0;j<cols;++j) 170 for(int i=0;i<rows;++i) 171 VERIFY(m.coeff(i,j) == data[(MatrixType::Flags&RowMajorBit) 172 ? i*outerStride + j*innerStride 173 : j*outerStride + i*innerStride]); 174 } 175 176 if(MatrixType::IsVectorAtCompileTime) 177 { 178 VERIFY(innerStride == int((&m.coeff(1))-(&m.coeff(0)))); 179 for (int i=0;i<size;++i) 180 VERIFY(m.coeff(i) == data[i*innerStride]); 181 } 182 } 183 184 template<typename MatrixType> 185 void data_and_stride(const MatrixType& m) 186 { 187 typedef typename MatrixType::Index Index; 188 Index rows = m.rows(); 189 Index cols = m.cols(); 190 191 Index r1 = internal::random<Index>(0,rows-1); 192 Index r2 = internal::random<Index>(r1,rows-1); 193 Index c1 = internal::random<Index>(0,cols-1); 194 Index c2 = internal::random<Index>(c1,cols-1); 195 196 MatrixType m1 = MatrixType::Random(rows, cols); 197 compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1)); 198 compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1)); 199 compare_using_data_and_stride(m1.row(r1)); 200 compare_using_data_and_stride(m1.col(c1)); 201 compare_using_data_and_stride(m1.row(r1).transpose()); 202 compare_using_data_and_stride(m1.col(c1).transpose()); 203 } 204 205 void test_block() 206 { 207 for(int i = 0; i < g_repeat; i++) { 208 CALL_SUBTEST_1( block(Matrix<float, 1, 1>()) ); 209 CALL_SUBTEST_2( block(Matrix4d()) ); 210 CALL_SUBTEST_3( block(MatrixXcf(3, 3)) ); 211 CALL_SUBTEST_4( block(MatrixXi(8, 12)) ); 212 CALL_SUBTEST_5( block(MatrixXcd(20, 20)) ); 213 CALL_SUBTEST_6( block(MatrixXf(20, 20)) ); 214 215 CALL_SUBTEST_8( block(Matrix<float,Dynamic,4>(3, 4)) ); 216 217 #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR 218 CALL_SUBTEST_6( data_and_stride(MatrixXf(internal::random(5,50), internal::random(5,50))) ); 219 CALL_SUBTEST_7( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(internal::random(5,50), internal::random(5,50))) ); 220 #endif 221 } 222 } 223
