sparse_block.cpp (12153B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr> 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 #include "sparse.h" 11 #include "AnnoyingScalar.h" 12 13 template<typename T> 14 typename Eigen::internal::enable_if<(T::Flags&RowMajorBit)==RowMajorBit, typename T::RowXpr>::type 15 innervec(T& A, Index i) 16 { 17 return A.row(i); 18 } 19 20 template<typename T> 21 typename Eigen::internal::enable_if<(T::Flags&RowMajorBit)==0, typename T::ColXpr>::type 22 innervec(T& A, Index i) 23 { 24 return A.col(i); 25 } 26 27 template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& ref) 28 { 29 const Index rows = ref.rows(); 30 const Index cols = ref.cols(); 31 const Index inner = ref.innerSize(); 32 const Index outer = ref.outerSize(); 33 34 typedef typename SparseMatrixType::Scalar Scalar; 35 typedef typename SparseMatrixType::RealScalar RealScalar; 36 typedef typename SparseMatrixType::StorageIndex StorageIndex; 37 38 double density = (std::max)(8./(rows*cols), 0.01); 39 typedef Matrix<Scalar,Dynamic,Dynamic,SparseMatrixType::IsRowMajor?RowMajor:ColMajor> DenseMatrix; 40 typedef Matrix<Scalar,Dynamic,1> DenseVector; 41 typedef Matrix<Scalar,1,Dynamic> RowDenseVector; 42 typedef SparseVector<Scalar> SparseVectorType; 43 44 Scalar s1 = internal::random<Scalar>(); 45 { 46 SparseMatrixType m(rows, cols); 47 DenseMatrix refMat = DenseMatrix::Zero(rows, cols); 48 initSparse<Scalar>(density, refMat, m); 49 50 VERIFY_IS_APPROX(m, refMat); 51 52 // test InnerIterators and Block expressions 53 for (int t=0; t<10; ++t) 54 { 55 Index j = internal::random<Index>(0,cols-2); 56 Index i = internal::random<Index>(0,rows-2); 57 Index w = internal::random<Index>(1,cols-j); 58 Index h = internal::random<Index>(1,rows-i); 59 60 VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); 61 for(Index c=0; c<w; c++) 62 { 63 VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); 64 for(Index r=0; r<h; r++) 65 { 66 VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); 67 VERIFY_IS_APPROX(m.block(i,j,h,w).coeff(r,c), refMat.block(i,j,h,w).coeff(r,c)); 68 } 69 } 70 for(Index r=0; r<h; r++) 71 { 72 VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); 73 for(Index c=0; c<w; c++) 74 { 75 VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); 76 VERIFY_IS_APPROX(m.block(i,j,h,w).coeff(r,c), refMat.block(i,j,h,w).coeff(r,c)); 77 } 78 } 79 80 VERIFY_IS_APPROX(m.middleCols(j,w), refMat.middleCols(j,w)); 81 VERIFY_IS_APPROX(m.middleRows(i,h), refMat.middleRows(i,h)); 82 for(Index r=0; r<h; r++) 83 { 84 VERIFY_IS_APPROX(m.middleCols(j,w).row(r), refMat.middleCols(j,w).row(r)); 85 VERIFY_IS_APPROX(m.middleRows(i,h).row(r), refMat.middleRows(i,h).row(r)); 86 for(Index c=0; c<w; c++) 87 { 88 VERIFY_IS_APPROX(m.col(c).coeff(r), refMat.col(c).coeff(r)); 89 VERIFY_IS_APPROX(m.row(r).coeff(c), refMat.row(r).coeff(c)); 90 91 VERIFY_IS_APPROX(m.middleCols(j,w).coeff(r,c), refMat.middleCols(j,w).coeff(r,c)); 92 VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c)); 93 if(m.middleCols(j,w).coeff(r,c) != Scalar(0)) 94 { 95 VERIFY_IS_APPROX(m.middleCols(j,w).coeffRef(r,c), refMat.middleCols(j,w).coeff(r,c)); 96 } 97 if(m.middleRows(i,h).coeff(r,c) != Scalar(0)) 98 { 99 VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c)); 100 } 101 } 102 } 103 for(Index c=0; c<w; c++) 104 { 105 VERIFY_IS_APPROX(m.middleCols(j,w).col(c), refMat.middleCols(j,w).col(c)); 106 VERIFY_IS_APPROX(m.middleRows(i,h).col(c), refMat.middleRows(i,h).col(c)); 107 } 108 } 109 110 for(Index c=0; c<cols; c++) 111 { 112 VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c)); 113 VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c)); 114 } 115 116 for(Index r=0; r<rows; r++) 117 { 118 VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r)); 119 VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r)); 120 } 121 } 122 123 // test innerVector() 124 { 125 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 126 SparseMatrixType m2(rows, cols); 127 initSparse<Scalar>(density, refMat2, m2); 128 Index j0 = internal::random<Index>(0,outer-1); 129 Index j1 = internal::random<Index>(0,outer-1); 130 Index r0 = internal::random<Index>(0,rows-1); 131 Index c0 = internal::random<Index>(0,cols-1); 132 133 VERIFY_IS_APPROX(m2.innerVector(j0), innervec(refMat2,j0)); 134 VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), innervec(refMat2,j0)+innervec(refMat2,j1)); 135 136 m2.innerVector(j0) *= Scalar(2); 137 innervec(refMat2,j0) *= Scalar(2); 138 VERIFY_IS_APPROX(m2, refMat2); 139 140 m2.row(r0) *= Scalar(3); 141 refMat2.row(r0) *= Scalar(3); 142 VERIFY_IS_APPROX(m2, refMat2); 143 144 m2.col(c0) *= Scalar(4); 145 refMat2.col(c0) *= Scalar(4); 146 VERIFY_IS_APPROX(m2, refMat2); 147 148 m2.row(r0) /= Scalar(3); 149 refMat2.row(r0) /= Scalar(3); 150 VERIFY_IS_APPROX(m2, refMat2); 151 152 m2.col(c0) /= Scalar(4); 153 refMat2.col(c0) /= Scalar(4); 154 VERIFY_IS_APPROX(m2, refMat2); 155 156 SparseVectorType v1; 157 VERIFY_IS_APPROX(v1 = m2.col(c0) * 4, refMat2.col(c0)*4); 158 VERIFY_IS_APPROX(v1 = m2.row(r0) * 4, refMat2.row(r0).transpose()*4); 159 160 SparseMatrixType m3(rows,cols); 161 m3.reserve(VectorXi::Constant(outer,int(inner/2))); 162 for(Index j=0; j<outer; ++j) 163 for(Index k=0; k<(std::min)(j,inner); ++k) 164 m3.insertByOuterInner(j,k) = internal::convert_index<StorageIndex>(k+1); 165 for(Index j=0; j<(std::min)(outer, inner); ++j) 166 { 167 VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); 168 if(j>0) 169 VERIFY(RealScalar(j)==numext::real(m3.innerVector(j).lastCoeff())); 170 } 171 m3.makeCompressed(); 172 for(Index j=0; j<(std::min)(outer, inner); ++j) 173 { 174 VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); 175 if(j>0) 176 VERIFY(RealScalar(j)==numext::real(m3.innerVector(j).lastCoeff())); 177 } 178 179 VERIFY(m3.innerVector(j0).nonZeros() == m3.transpose().innerVector(j0).nonZeros()); 180 181 // m2.innerVector(j0) = 2*m2.innerVector(j1); 182 // refMat2.col(j0) = 2*refMat2.col(j1); 183 // VERIFY_IS_APPROX(m2, refMat2); 184 } 185 186 // test innerVectors() 187 { 188 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 189 SparseMatrixType m2(rows, cols); 190 initSparse<Scalar>(density, refMat2, m2); 191 if(internal::random<float>(0,1)>0.5f) m2.makeCompressed(); 192 Index j0 = internal::random<Index>(0,outer-2); 193 Index j1 = internal::random<Index>(0,outer-2); 194 Index n0 = internal::random<Index>(1,outer-(std::max)(j0,j1)); 195 if(SparseMatrixType::IsRowMajor) 196 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); 197 else 198 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); 199 if(SparseMatrixType::IsRowMajor) 200 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), 201 refMat2.middleRows(j0,n0)+refMat2.middleRows(j1,n0)); 202 else 203 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), 204 refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); 205 206 VERIFY_IS_APPROX(m2, refMat2); 207 208 VERIFY(m2.innerVectors(j0,n0).nonZeros() == m2.transpose().innerVectors(j0,n0).nonZeros()); 209 210 m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); 211 if(SparseMatrixType::IsRowMajor) 212 refMat2.middleRows(j0,n0) = (refMat2.middleRows(j0,n0) + refMat2.middleRows(j1,n0)).eval(); 213 else 214 refMat2.middleCols(j0,n0) = (refMat2.middleCols(j0,n0) + refMat2.middleCols(j1,n0)).eval(); 215 216 VERIFY_IS_APPROX(m2, refMat2); 217 } 218 219 // test generic blocks 220 { 221 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 222 SparseMatrixType m2(rows, cols); 223 initSparse<Scalar>(density, refMat2, m2); 224 Index j0 = internal::random<Index>(0,outer-2); 225 Index j1 = internal::random<Index>(0,outer-2); 226 Index n0 = internal::random<Index>(1,outer-(std::max)(j0,j1)); 227 if(SparseMatrixType::IsRowMajor) 228 VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols)); 229 else 230 VERIFY_IS_APPROX(m2.block(0,j0,rows,n0), refMat2.block(0,j0,rows,n0)); 231 232 if(SparseMatrixType::IsRowMajor) 233 VERIFY_IS_APPROX(m2.block(j0,0,n0,cols)+m2.block(j1,0,n0,cols), 234 refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols)); 235 else 236 VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0), 237 refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); 238 239 Index i = internal::random<Index>(0,m2.outerSize()-1); 240 if(SparseMatrixType::IsRowMajor) { 241 m2.innerVector(i) = m2.innerVector(i) * s1; 242 refMat2.row(i) = refMat2.row(i) * s1; 243 VERIFY_IS_APPROX(m2,refMat2); 244 } else { 245 m2.innerVector(i) = m2.innerVector(i) * s1; 246 refMat2.col(i) = refMat2.col(i) * s1; 247 VERIFY_IS_APPROX(m2,refMat2); 248 } 249 250 Index r0 = internal::random<Index>(0,rows-2); 251 Index c0 = internal::random<Index>(0,cols-2); 252 Index r1 = internal::random<Index>(1,rows-r0); 253 Index c1 = internal::random<Index>(1,cols-c0); 254 255 VERIFY_IS_APPROX(DenseVector(m2.col(c0)), refMat2.col(c0)); 256 VERIFY_IS_APPROX(m2.col(c0), refMat2.col(c0)); 257 258 VERIFY_IS_APPROX(RowDenseVector(m2.row(r0)), refMat2.row(r0)); 259 VERIFY_IS_APPROX(m2.row(r0), refMat2.row(r0)); 260 261 VERIFY_IS_APPROX(m2.block(r0,c0,r1,c1), refMat2.block(r0,c0,r1,c1)); 262 VERIFY_IS_APPROX((2*m2).block(r0,c0,r1,c1), (2*refMat2).block(r0,c0,r1,c1)); 263 264 if(m2.nonZeros()>0) 265 { 266 VERIFY_IS_APPROX(m2, refMat2); 267 SparseMatrixType m3(rows, cols); 268 DenseMatrix refMat3(rows, cols); refMat3.setZero(); 269 Index n = internal::random<Index>(1,10); 270 for(Index k=0; k<n; ++k) 271 { 272 Index o1 = internal::random<Index>(0,outer-1); 273 Index o2 = internal::random<Index>(0,outer-1); 274 if(SparseMatrixType::IsRowMajor) 275 { 276 m3.innerVector(o1) = m2.row(o2); 277 refMat3.row(o1) = refMat2.row(o2); 278 } 279 else 280 { 281 m3.innerVector(o1) = m2.col(o2); 282 refMat3.col(o1) = refMat2.col(o2); 283 } 284 if(internal::random<bool>()) 285 m3.makeCompressed(); 286 } 287 if(m3.nonZeros()>0) 288 VERIFY_IS_APPROX(m3, refMat3); 289 } 290 } 291 } 292 293 EIGEN_DECLARE_TEST(sparse_block) 294 { 295 for(int i = 0; i < g_repeat; i++) { 296 int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200); 297 if(Eigen::internal::random<int>(0,4) == 0) { 298 r = c; // check square matrices in 25% of tries 299 } 300 EIGEN_UNUSED_VARIABLE(r+c); 301 CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(1, 1)) )); 302 CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(8, 8)) )); 303 CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(r, c)) )); 304 CALL_SUBTEST_2(( sparse_block(SparseMatrix<std::complex<double>, ColMajor>(r, c)) )); 305 CALL_SUBTEST_2(( sparse_block(SparseMatrix<std::complex<double>, RowMajor>(r, c)) )); 306 307 CALL_SUBTEST_3(( sparse_block(SparseMatrix<double,ColMajor,long int>(r, c)) )); 308 CALL_SUBTEST_3(( sparse_block(SparseMatrix<double,RowMajor,long int>(r, c)) )); 309 310 r = Eigen::internal::random<int>(1,100); 311 c = Eigen::internal::random<int>(1,100); 312 if(Eigen::internal::random<int>(0,4) == 0) { 313 r = c; // check square matrices in 25% of tries 314 } 315 316 CALL_SUBTEST_4(( sparse_block(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) )); 317 CALL_SUBTEST_4(( sparse_block(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) )); 318 #ifndef EIGEN_TEST_ANNOYING_SCALAR_DONT_THROW 319 AnnoyingScalar::dont_throw = true; 320 #endif 321 CALL_SUBTEST_5(( sparse_block(SparseMatrix<AnnoyingScalar>(r,c)) )); 322 } 323 }