1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro (at) gmail.com> 6 // Copyright (C) 2013 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 static long g_realloc_count = 0; 13 #define EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN g_realloc_count++; 14 15 #include "sparse.h" 16 17 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) 18 { 19 typedef typename SparseMatrixType::StorageIndex StorageIndex; 20 typedef Matrix<StorageIndex,2,1> Vector2; 21 22 const Index rows = ref.rows(); 23 const Index cols = ref.cols(); 24 //const Index inner = ref.innerSize(); 25 //const Index outer = ref.outerSize(); 26 27 typedef typename SparseMatrixType::Scalar Scalar; 28 typedef typename SparseMatrixType::RealScalar RealScalar; 29 enum { Flags = SparseMatrixType::Flags }; 30 31 double density = (std::max)(8./(rows*cols), 0.01); 32 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 33 typedef Matrix<Scalar,Dynamic,1> DenseVector; 34 Scalar eps = 1e-6; 35 36 Scalar s1 = internal::random<Scalar>(); 37 { 38 SparseMatrixType m(rows, cols); 39 DenseMatrix refMat = DenseMatrix::Zero(rows, cols); 40 DenseVector vec1 = DenseVector::Random(rows); 41 42 std::vector<Vector2> zeroCoords; 43 std::vector<Vector2> nonzeroCoords; 44 initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); 45 46 // test coeff and coeffRef 47 for (std::size_t i=0; i<zeroCoords.size(); ++i) 48 { 49 VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps ); 50 if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value) 51 VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[i].x(),zeroCoords[i].y()) = 5 ); 52 } 53 VERIFY_IS_APPROX(m, refMat); 54 55 if(!nonzeroCoords.empty()) { 56 m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); 57 refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); 58 } 59 60 VERIFY_IS_APPROX(m, refMat); 61 62 // test assertion 63 VERIFY_RAISES_ASSERT( m.coeffRef(-1,1) = 0 ); 64 VERIFY_RAISES_ASSERT( m.coeffRef(0,m.cols()) = 0 ); 65 } 66 67 // test insert (inner random) 68 { 69 DenseMatrix m1(rows,cols); 70 m1.setZero(); 71 SparseMatrixType m2(rows,cols); 72 bool call_reserve = internal::random<int>()%2; 73 Index nnz = internal::random<int>(1,int(rows)/2); 74 if(call_reserve) 75 { 76 if(internal::random<int>()%2) 77 m2.reserve(VectorXi::Constant(m2.outerSize(), int(nnz))); 78 else 79 m2.reserve(m2.outerSize() * nnz); 80 } 81 g_realloc_count = 0; 82 for (Index j=0; j<cols; ++j) 83 { 84 for (Index k=0; k<nnz; ++k) 85 { 86 Index i = internal::random<Index>(0,rows-1); 87 if (m1.coeff(i,j)==Scalar(0)) 88 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 89 } 90 } 91 92 if(call_reserve && !SparseMatrixType::IsRowMajor) 93 { 94 VERIFY(g_realloc_count==0); 95 } 96 97 m2.finalize(); 98 VERIFY_IS_APPROX(m2,m1); 99 } 100 101 // test insert (fully random) 102 { 103 DenseMatrix m1(rows,cols); 104 m1.setZero(); 105 SparseMatrixType m2(rows,cols); 106 if(internal::random<int>()%2) 107 m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); 108 for (int k=0; k<rows*cols; ++k) 109 { 110 Index i = internal::random<Index>(0,rows-1); 111 Index j = internal::random<Index>(0,cols-1); 112 if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2)) 113 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 114 else 115 { 116 Scalar v = internal::random<Scalar>(); 117 m2.coeffRef(i,j) += v; 118 m1(i,j) += v; 119 } 120 } 121 VERIFY_IS_APPROX(m2,m1); 122 } 123 124 // test insert (un-compressed) 125 for(int mode=0;mode<4;++mode) 126 { 127 DenseMatrix m1(rows,cols); 128 m1.setZero(); 129 SparseMatrixType m2(rows,cols); 130 VectorXi r(VectorXi::Constant(m2.outerSize(), ((mode%2)==0) ? int(m2.innerSize()) : std::max<int>(1,int(m2.innerSize())/8))); 131 m2.reserve(r); 132 for (Index k=0; k<rows*cols; ++k) 133 { 134 Index i = internal::random<Index>(0,rows-1); 135 Index j = internal::random<Index>(0,cols-1); 136 if (m1.coeff(i,j)==Scalar(0)) 137 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 138 if(mode==3) 139 m2.reserve(r); 140 } 141 if(internal::random<int>()%2) 142 m2.makeCompressed(); 143 VERIFY_IS_APPROX(m2,m1); 144 } 145 146 // test basic computations 147 { 148 DenseMatrix refM1 = DenseMatrix::Zero(rows, cols); 149 DenseMatrix refM2 = DenseMatrix::Zero(rows, cols); 150 DenseMatrix refM3 = DenseMatrix::Zero(rows, cols); 151 DenseMatrix refM4 = DenseMatrix::Zero(rows, cols); 152 SparseMatrixType m1(rows, cols); 153 SparseMatrixType m2(rows, cols); 154 SparseMatrixType m3(rows, cols); 155 SparseMatrixType m4(rows, cols); 156 initSparse<Scalar>(density, refM1, m1); 157 initSparse<Scalar>(density, refM2, m2); 158 initSparse<Scalar>(density, refM3, m3); 159 initSparse<Scalar>(density, refM4, m4); 160 161 if(internal::random<bool>()) 162 m1.makeCompressed(); 163 164 Index m1_nnz = m1.nonZeros(); 165 166 VERIFY_IS_APPROX(m1*s1, refM1*s1); 167 VERIFY_IS_APPROX(m1+m2, refM1+refM2); 168 VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); 169 VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2)); 170 VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); 171 VERIFY_IS_APPROX(m4=m1/s1, refM1/s1); 172 VERIFY_IS_EQUAL(m4.nonZeros(), m1_nnz); 173 174 if(SparseMatrixType::IsRowMajor) 175 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); 176 else 177 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.col(0)), refM1.col(0).dot(refM2.col(0))); 178 179 DenseVector rv = DenseVector::Random(m1.cols()); 180 DenseVector cv = DenseVector::Random(m1.rows()); 181 Index r = internal::random<Index>(0,m1.rows()-2); 182 Index c = internal::random<Index>(0,m1.cols()-1); 183 VERIFY_IS_APPROX(( m1.template block<1,Dynamic>(r,0,1,m1.cols()).dot(rv)) , refM1.row(r).dot(rv)); 184 VERIFY_IS_APPROX(m1.row(r).dot(rv), refM1.row(r).dot(rv)); 185 VERIFY_IS_APPROX(m1.col(c).dot(cv), refM1.col(c).dot(cv)); 186 187 VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate()); 188 VERIFY_IS_APPROX(m1.real(), refM1.real()); 189 190 refM4.setRandom(); 191 // sparse cwise* dense 192 VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4)); 193 // dense cwise* sparse 194 VERIFY_IS_APPROX(refM4.cwiseProduct(m3), refM4.cwiseProduct(refM3)); 195 // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); 196 197 VERIFY_IS_APPROX(refM4 + m3, refM4 + refM3); 198 VERIFY_IS_APPROX(m3 + refM4, refM3 + refM4); 199 VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3); 200 VERIFY_IS_APPROX(m3 - refM4, refM3 - refM4); 201 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3); 202 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3*RealScalar(0.5)).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3); 203 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3.cwiseProduct(m3)).eval(), RealScalar(0.5)*refM4 + refM3.cwiseProduct(refM3)); 204 205 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3); 206 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3*RealScalar(0.5)).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3); 207 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (m3+m3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3)); 208 VERIFY_IS_APPROX(((refM3+m3)+RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM3 + (refM3+refM3)); 209 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (refM3+m3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3)); 210 VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (m3+refM3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3)); 211 212 213 VERIFY_IS_APPROX(m1.sum(), refM1.sum()); 214 215 m4 = m1; refM4 = m4; 216 217 VERIFY_IS_APPROX(m1*=s1, refM1*=s1); 218 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 219 VERIFY_IS_APPROX(m1/=s1, refM1/=s1); 220 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 221 222 VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); 223 VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); 224 225 if (rows>=2 && cols>=2) 226 { 227 VERIFY_RAISES_ASSERT( m1 += m1.innerVector(0) ); 228 VERIFY_RAISES_ASSERT( m1 -= m1.innerVector(0) ); 229 VERIFY_RAISES_ASSERT( refM1 -= m1.innerVector(0) ); 230 VERIFY_RAISES_ASSERT( refM1 += m1.innerVector(0) ); 231 m1 = m4; refM1 = refM4; 232 } 233 234 // test aliasing 235 VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1)); 236 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 237 m1 = m4; refM1 = refM4; 238 VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval())); 239 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 240 m1 = m4; refM1 = refM4; 241 VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval())); 242 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 243 m1 = m4; refM1 = refM4; 244 VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1)); 245 VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz); 246 m1 = m4; refM1 = refM4; 247 248 if(m1.isCompressed()) 249 { 250 VERIFY_IS_APPROX(m1.coeffs().sum(), m1.sum()); 251 m1.coeffs() += s1; 252 for(Index j = 0; j<m1.outerSize(); ++j) 253 for(typename SparseMatrixType::InnerIterator it(m1,j); it; ++it) 254 refM1(it.row(), it.col()) += s1; 255 VERIFY_IS_APPROX(m1, refM1); 256 } 257 258 // and/or 259 { 260 typedef SparseMatrix<bool, SparseMatrixType::Options, typename SparseMatrixType::StorageIndex> SpBool; 261 SpBool mb1 = m1.real().template cast<bool>(); 262 SpBool mb2 = m2.real().template cast<bool>(); 263 VERIFY_IS_EQUAL(mb1.template cast<int>().sum(), refM1.real().template cast<bool>().count()); 264 VERIFY_IS_EQUAL((mb1 && mb2).template cast<int>().sum(), (refM1.real().template cast<bool>() && refM2.real().template cast<bool>()).count()); 265 VERIFY_IS_EQUAL((mb1 || mb2).template cast<int>().sum(), (refM1.real().template cast<bool>() || refM2.real().template cast<bool>()).count()); 266 SpBool mb3 = mb1 && mb2; 267 if(mb1.coeffs().all() && mb2.coeffs().all()) 268 { 269 VERIFY_IS_EQUAL(mb3.nonZeros(), (refM1.real().template cast<bool>() && refM2.real().template cast<bool>()).count()); 270 } 271 } 272 } 273 274 // test reverse iterators 275 { 276 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 277 SparseMatrixType m2(rows, cols); 278 initSparse<Scalar>(density, refMat2, m2); 279 std::vector<Scalar> ref_value(m2.innerSize()); 280 std::vector<Index> ref_index(m2.innerSize()); 281 if(internal::random<bool>()) 282 m2.makeCompressed(); 283 for(Index j = 0; j<m2.outerSize(); ++j) 284 { 285 Index count_forward = 0; 286 287 for(typename SparseMatrixType::InnerIterator it(m2,j); it; ++it) 288 { 289 ref_value[ref_value.size()-1-count_forward] = it.value(); 290 ref_index[ref_index.size()-1-count_forward] = it.index(); 291 count_forward++; 292 } 293 Index count_reverse = 0; 294 for(typename SparseMatrixType::ReverseInnerIterator it(m2,j); it; --it) 295 { 296 VERIFY_IS_APPROX( std::abs(ref_value[ref_value.size()-count_forward+count_reverse])+1, std::abs(it.value())+1); 297 VERIFY_IS_EQUAL( ref_index[ref_index.size()-count_forward+count_reverse] , it.index()); 298 count_reverse++; 299 } 300 VERIFY_IS_EQUAL(count_forward, count_reverse); 301 } 302 } 303 304 // test transpose 305 { 306 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 307 SparseMatrixType m2(rows, cols); 308 initSparse<Scalar>(density, refMat2, m2); 309 VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); 310 VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); 311 312 VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); 313 314 // check isApprox handles opposite storage order 315 typename Transpose<SparseMatrixType>::PlainObject m3(m2); 316 VERIFY(m2.isApprox(m3)); 317 } 318 319 // test prune 320 { 321 SparseMatrixType m2(rows, cols); 322 DenseMatrix refM2(rows, cols); 323 refM2.setZero(); 324 int countFalseNonZero = 0; 325 int countTrueNonZero = 0; 326 m2.reserve(VectorXi::Constant(m2.outerSize(), int(m2.innerSize()))); 327 for (Index j=0; j<m2.cols(); ++j) 328 { 329 for (Index i=0; i<m2.rows(); ++i) 330 { 331 float x = internal::random<float>(0,1); 332 if (x<0.1f) 333 { 334 // do nothing 335 } 336 else if (x<0.5f) 337 { 338 countFalseNonZero++; 339 m2.insert(i,j) = Scalar(0); 340 } 341 else 342 { 343 countTrueNonZero++; 344 m2.insert(i,j) = Scalar(1); 345 refM2(i,j) = Scalar(1); 346 } 347 } 348 } 349 if(internal::random<bool>()) 350 m2.makeCompressed(); 351 VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros()); 352 if(countTrueNonZero>0) 353 VERIFY_IS_APPROX(m2, refM2); 354 m2.prune(Scalar(1)); 355 VERIFY(countTrueNonZero==m2.nonZeros()); 356 VERIFY_IS_APPROX(m2, refM2); 357 } 358 359 // test setFromTriplets 360 { 361 typedef Triplet<Scalar,StorageIndex> TripletType; 362 std::vector<TripletType> triplets; 363 Index ntriplets = rows*cols; 364 triplets.reserve(ntriplets); 365 DenseMatrix refMat_sum = DenseMatrix::Zero(rows,cols); 366 DenseMatrix refMat_prod = DenseMatrix::Zero(rows,cols); 367 DenseMatrix refMat_last = DenseMatrix::Zero(rows,cols); 368 369 for(Index i=0;i<ntriplets;++i) 370 { 371 StorageIndex r = internal::random<StorageIndex>(0,StorageIndex(rows-1)); 372 StorageIndex c = internal::random<StorageIndex>(0,StorageIndex(cols-1)); 373 Scalar v = internal::random<Scalar>(); 374 triplets.push_back(TripletType(r,c,v)); 375 refMat_sum(r,c) += v; 376 if(std::abs(refMat_prod(r,c))==0) 377 refMat_prod(r,c) = v; 378 else 379 refMat_prod(r,c) *= v; 380 refMat_last(r,c) = v; 381 } 382 SparseMatrixType m(rows,cols); 383 m.setFromTriplets(triplets.begin(), triplets.end()); 384 VERIFY_IS_APPROX(m, refMat_sum); 385 386 m.setFromTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>()); 387 VERIFY_IS_APPROX(m, refMat_prod); 388 #if (defined(__cplusplus) && __cplusplus >= 201103L) 389 m.setFromTriplets(triplets.begin(), triplets.end(), [] (Scalar,Scalar b) { return b; }); 390 VERIFY_IS_APPROX(m, refMat_last); 391 #endif 392 } 393 394 // test Map 395 { 396 DenseMatrix refMat2(rows, cols), refMat3(rows, cols); 397 SparseMatrixType m2(rows, cols), m3(rows, cols); 398 initSparse<Scalar>(density, refMat2, m2); 399 initSparse<Scalar>(density, refMat3, m3); 400 { 401 Map<SparseMatrixType> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr()); 402 Map<SparseMatrixType> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr()); 403 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); 404 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); 405 } 406 { 407 MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr()); 408 MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr()); 409 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); 410 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); 411 } 412 413 Index i = internal::random<Index>(0,rows-1); 414 Index j = internal::random<Index>(0,cols-1); 415 m2.coeffRef(i,j) = 123; 416 if(internal::random<bool>()) 417 m2.makeCompressed(); 418 Map<SparseMatrixType> mapMat2(rows, cols, m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr()); 419 VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(123)); 420 VERIFY_IS_EQUAL(mapMat2.coeff(i,j),Scalar(123)); 421 mapMat2.coeffRef(i,j) = -123; 422 VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(-123)); 423 } 424 425 // test triangularView 426 { 427 DenseMatrix refMat2(rows, cols), refMat3(rows, cols); 428 SparseMatrixType m2(rows, cols), m3(rows, cols); 429 initSparse<Scalar>(density, refMat2, m2); 430 refMat3 = refMat2.template triangularView<Lower>(); 431 m3 = m2.template triangularView<Lower>(); 432 VERIFY_IS_APPROX(m3, refMat3); 433 434 refMat3 = refMat2.template triangularView<Upper>(); 435 m3 = m2.template triangularView<Upper>(); 436 VERIFY_IS_APPROX(m3, refMat3); 437 438 { 439 refMat3 = refMat2.template triangularView<UnitUpper>(); 440 m3 = m2.template triangularView<UnitUpper>(); 441 VERIFY_IS_APPROX(m3, refMat3); 442 443 refMat3 = refMat2.template triangularView<UnitLower>(); 444 m3 = m2.template triangularView<UnitLower>(); 445 VERIFY_IS_APPROX(m3, refMat3); 446 } 447 448 refMat3 = refMat2.template triangularView<StrictlyUpper>(); 449 m3 = m2.template triangularView<StrictlyUpper>(); 450 VERIFY_IS_APPROX(m3, refMat3); 451 452 refMat3 = refMat2.template triangularView<StrictlyLower>(); 453 m3 = m2.template triangularView<StrictlyLower>(); 454 VERIFY_IS_APPROX(m3, refMat3); 455 456 // check sparse-triangular to dense 457 refMat3 = m2.template triangularView<StrictlyUpper>(); 458 VERIFY_IS_APPROX(refMat3, DenseMatrix(refMat2.template triangularView<StrictlyUpper>())); 459 } 460 461 // test selfadjointView 462 if(!SparseMatrixType::IsRowMajor) 463 { 464 DenseMatrix refMat2(rows, rows), refMat3(rows, rows); 465 SparseMatrixType m2(rows, rows), m3(rows, rows); 466 initSparse<Scalar>(density, refMat2, m2); 467 refMat3 = refMat2.template selfadjointView<Lower>(); 468 m3 = m2.template selfadjointView<Lower>(); 469 VERIFY_IS_APPROX(m3, refMat3); 470 471 refMat3 += refMat2.template selfadjointView<Lower>(); 472 m3 += m2.template selfadjointView<Lower>(); 473 VERIFY_IS_APPROX(m3, refMat3); 474 475 refMat3 -= refMat2.template selfadjointView<Lower>(); 476 m3 -= m2.template selfadjointView<Lower>(); 477 VERIFY_IS_APPROX(m3, refMat3); 478 479 // selfadjointView only works for square matrices: 480 SparseMatrixType m4(rows, rows+1); 481 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Lower>()); 482 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Upper>()); 483 } 484 485 // test sparseView 486 { 487 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 488 SparseMatrixType m2(rows, rows); 489 initSparse<Scalar>(density, refMat2, m2); 490 VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval()); 491 492 // sparse view on expressions: 493 VERIFY_IS_APPROX((s1*m2).eval(), (s1*refMat2).sparseView().eval()); 494 VERIFY_IS_APPROX((m2+m2).eval(), (refMat2+refMat2).sparseView().eval()); 495 VERIFY_IS_APPROX((m2*m2).eval(), (refMat2.lazyProduct(refMat2)).sparseView().eval()); 496 VERIFY_IS_APPROX((m2*m2).eval(), (refMat2*refMat2).sparseView().eval()); 497 } 498 499 // test diagonal 500 { 501 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 502 SparseMatrixType m2(rows, cols); 503 initSparse<Scalar>(density, refMat2, m2); 504 VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval()); 505 DenseVector d = m2.diagonal(); 506 VERIFY_IS_APPROX(d, refMat2.diagonal().eval()); 507 d = m2.diagonal().array(); 508 VERIFY_IS_APPROX(d, refMat2.diagonal().eval()); 509 VERIFY_IS_APPROX(const_cast<const SparseMatrixType&>(m2).diagonal(), refMat2.diagonal().eval()); 510 511 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag); 512 m2.diagonal() += refMat2.diagonal(); 513 refMat2.diagonal() += refMat2.diagonal(); 514 VERIFY_IS_APPROX(m2, refMat2); 515 } 516 517 // test diagonal to sparse 518 { 519 DenseVector d = DenseVector::Random(rows); 520 DenseMatrix refMat2 = d.asDiagonal(); 521 SparseMatrixType m2(rows, rows); 522 m2 = d.asDiagonal(); 523 VERIFY_IS_APPROX(m2, refMat2); 524 SparseMatrixType m3(d.asDiagonal()); 525 VERIFY_IS_APPROX(m3, refMat2); 526 refMat2 += d.asDiagonal(); 527 m2 += d.asDiagonal(); 528 VERIFY_IS_APPROX(m2, refMat2); 529 } 530 531 // test conservative resize 532 { 533 std::vector< std::pair<StorageIndex,StorageIndex> > inc; 534 if(rows > 3 && cols > 2) 535 inc.push_back(std::pair<StorageIndex,StorageIndex>(-3,-2)); 536 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,0)); 537 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,2)); 538 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,0)); 539 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,3)); 540 541 for(size_t i = 0; i< inc.size(); i++) { 542 StorageIndex incRows = inc[i].first; 543 StorageIndex incCols = inc[i].second; 544 SparseMatrixType m1(rows, cols); 545 DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols); 546 initSparse<Scalar>(density, refMat1, m1); 547 548 m1.conservativeResize(rows+incRows, cols+incCols); 549 refMat1.conservativeResize(rows+incRows, cols+incCols); 550 if (incRows > 0) refMat1.bottomRows(incRows).setZero(); 551 if (incCols > 0) refMat1.rightCols(incCols).setZero(); 552 553 VERIFY_IS_APPROX(m1, refMat1); 554 555 // Insert new values 556 if (incRows > 0) 557 m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1; 558 if (incCols > 0) 559 m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1; 560 561 VERIFY_IS_APPROX(m1, refMat1); 562 563 564 } 565 } 566 567 // test Identity matrix 568 { 569 DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows); 570 SparseMatrixType m1(rows, rows); 571 m1.setIdentity(); 572 VERIFY_IS_APPROX(m1, refMat1); 573 for(int k=0; k<rows*rows/4; ++k) 574 { 575 Index i = internal::random<Index>(0,rows-1); 576 Index j = internal::random<Index>(0,rows-1); 577 Scalar v = internal::random<Scalar>(); 578 m1.coeffRef(i,j) = v; 579 refMat1.coeffRef(i,j) = v; 580 VERIFY_IS_APPROX(m1, refMat1); 581 if(internal::random<Index>(0,10)<2) 582 m1.makeCompressed(); 583 } 584 m1.setIdentity(); 585 refMat1.setIdentity(); 586 VERIFY_IS_APPROX(m1, refMat1); 587 } 588 589 // test array/vector of InnerIterator 590 { 591 typedef typename SparseMatrixType::InnerIterator IteratorType; 592 593 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 594 SparseMatrixType m2(rows, cols); 595 initSparse<Scalar>(density, refMat2, m2); 596 IteratorType static_array[2]; 597 static_array[0] = IteratorType(m2,0); 598 static_array[1] = IteratorType(m2,m2.outerSize()-1); 599 VERIFY( static_array[0] || m2.innerVector(static_array[0].outer()).nonZeros() == 0 ); 600 VERIFY( static_array[1] || m2.innerVector(static_array[1].outer()).nonZeros() == 0 ); 601 if(static_array[0] && static_array[1]) 602 { 603 ++(static_array[1]); 604 static_array[1] = IteratorType(m2,0); 605 VERIFY( static_array[1] ); 606 VERIFY( static_array[1].index() == static_array[0].index() ); 607 VERIFY( static_array[1].outer() == static_array[0].outer() ); 608 VERIFY( static_array[1].value() == static_array[0].value() ); 609 } 610 611 std::vector<IteratorType> iters(2); 612 iters[0] = IteratorType(m2,0); 613 iters[1] = IteratorType(m2,m2.outerSize()-1); 614 } 615 } 616 617 618 template<typename SparseMatrixType> 619 void big_sparse_triplet(Index rows, Index cols, double density) { 620 typedef typename SparseMatrixType::StorageIndex StorageIndex; 621 typedef typename SparseMatrixType::Scalar Scalar; 622 typedef Triplet<Scalar,Index> TripletType; 623 std::vector<TripletType> triplets; 624 double nelements = density * rows*cols; 625 VERIFY(nelements>=0 && nelements < NumTraits<StorageIndex>::highest()); 626 Index ntriplets = Index(nelements); 627 triplets.reserve(ntriplets); 628 Scalar sum = Scalar(0); 629 for(Index i=0;i<ntriplets;++i) 630 { 631 Index r = internal::random<Index>(0,rows-1); 632 Index c = internal::random<Index>(0,cols-1); 633 Scalar v = internal::random<Scalar>(); 634 triplets.push_back(TripletType(r,c,v)); 635 sum += v; 636 } 637 SparseMatrixType m(rows,cols); 638 m.setFromTriplets(triplets.begin(), triplets.end()); 639 VERIFY(m.nonZeros() <= ntriplets); 640 VERIFY_IS_APPROX(sum, m.sum()); 641 } 642 643 644 void test_sparse_basic() 645 { 646 for(int i = 0; i < g_repeat; i++) { 647 int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200); 648 if(Eigen::internal::random<int>(0,4) == 0) { 649 r = c; // check square matrices in 25% of tries 650 } 651 EIGEN_UNUSED_VARIABLE(r+c); 652 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(1, 1)) )); 653 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) )); 654 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(r, c)) )); 655 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(r, c)) )); 656 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(r, c)) )); 657 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,ColMajor,long int>(r, c)) )); 658 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,RowMajor,long int>(r, c)) )); 659 660 r = Eigen::internal::random<int>(1,100); 661 c = Eigen::internal::random<int>(1,100); 662 if(Eigen::internal::random<int>(0,4) == 0) { 663 r = c; // check square matrices in 25% of tries 664 } 665 666 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) )); 667 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) )); 668 } 669 670 // Regression test for bug 900: (manually insert higher values here, if you have enough RAM): 671 CALL_SUBTEST_3((big_sparse_triplet<SparseMatrix<float, RowMajor, int> >(10000, 10000, 0.125))); 672 CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125))); 673 674 // Regression test for bug 1105 675 #ifdef EIGEN_TEST_PART_7 676 { 677 int n = Eigen::internal::random<int>(200,600); 678 SparseMatrix<std::complex<double>,0, long> mat(n, n); 679 std::complex<double> val; 680 681 for(int i=0; i<n; ++i) 682 { 683 mat.coeffRef(i, i%(n/10)) = val; 684 VERIFY(mat.data().allocatedSize()<20*n); 685 } 686 } 687 #endif 688 } 689
