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 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #include "sparse.h" 12 13 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) 14 { 15 typedef typename SparseMatrixType::Index Index; 16 17 const Index rows = ref.rows(); 18 const Index cols = ref.cols(); 19 typedef typename SparseMatrixType::Scalar Scalar; 20 enum { Flags = SparseMatrixType::Flags }; 21 22 double density = (std::max)(8./(rows*cols), 0.01); 23 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 24 typedef Matrix<Scalar,Dynamic,1> DenseVector; 25 Scalar eps = 1e-6; 26 27 SparseMatrixType m(rows, cols); 28 DenseMatrix refMat = DenseMatrix::Zero(rows, cols); 29 DenseVector vec1 = DenseVector::Random(rows); 30 Scalar s1 = internal::random<Scalar>(); 31 32 std::vector<Vector2i> zeroCoords; 33 std::vector<Vector2i> nonzeroCoords; 34 initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); 35 36 if (zeroCoords.size()==0 || nonzeroCoords.size()==0) 37 return; 38 39 // test coeff and coeffRef 40 for (int i=0; i<(int)zeroCoords.size(); ++i) 41 { 42 VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps ); 43 if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value) 44 VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 ); 45 } 46 VERIFY_IS_APPROX(m, refMat); 47 48 m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); 49 refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); 50 51 VERIFY_IS_APPROX(m, refMat); 52 /* 53 // test InnerIterators and Block expressions 54 for (int t=0; t<10; ++t) 55 { 56 int j = internal::random<int>(0,cols-1); 57 int i = internal::random<int>(0,rows-1); 58 int w = internal::random<int>(1,cols-j-1); 59 int h = internal::random<int>(1,rows-i-1); 60 61 // VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); 62 for(int c=0; c<w; c++) 63 { 64 VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); 65 for(int r=0; r<h; r++) 66 { 67 // VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); 68 } 69 } 70 // for(int 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(int 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 // } 77 // } 78 } 79 80 for(int c=0; c<cols; c++) 81 { 82 VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c)); 83 VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c)); 84 } 85 86 for(int r=0; r<rows; r++) 87 { 88 VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r)); 89 VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r)); 90 } 91 */ 92 93 // test insert (inner random) 94 { 95 DenseMatrix m1(rows,cols); 96 m1.setZero(); 97 SparseMatrixType m2(rows,cols); 98 if(internal::random<int>()%2) 99 m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); 100 for (int j=0; j<cols; ++j) 101 { 102 for (int k=0; k<rows/2; ++k) 103 { 104 int i = internal::random<int>(0,rows-1); 105 if (m1.coeff(i,j)==Scalar(0)) 106 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 107 } 108 } 109 m2.finalize(); 110 VERIFY_IS_APPROX(m2,m1); 111 } 112 113 // test insert (fully random) 114 { 115 DenseMatrix m1(rows,cols); 116 m1.setZero(); 117 SparseMatrixType m2(rows,cols); 118 if(internal::random<int>()%2) 119 m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); 120 for (int k=0; k<rows*cols; ++k) 121 { 122 int i = internal::random<int>(0,rows-1); 123 int j = internal::random<int>(0,cols-1); 124 if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2)) 125 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 126 else 127 { 128 Scalar v = internal::random<Scalar>(); 129 m2.coeffRef(i,j) += v; 130 m1(i,j) += v; 131 } 132 } 133 VERIFY_IS_APPROX(m2,m1); 134 } 135 136 // test insert (un-compressed) 137 for(int mode=0;mode<4;++mode) 138 { 139 DenseMatrix m1(rows,cols); 140 m1.setZero(); 141 SparseMatrixType m2(rows,cols); 142 VectorXi r(VectorXi::Constant(m2.outerSize(), ((mode%2)==0) ? m2.innerSize() : std::max<int>(1,m2.innerSize()/8))); 143 m2.reserve(r); 144 for (int k=0; k<rows*cols; ++k) 145 { 146 int i = internal::random<int>(0,rows-1); 147 int j = internal::random<int>(0,cols-1); 148 if (m1.coeff(i,j)==Scalar(0)) 149 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); 150 if(mode==3) 151 m2.reserve(r); 152 } 153 if(internal::random<int>()%2) 154 m2.makeCompressed(); 155 VERIFY_IS_APPROX(m2,m1); 156 } 157 158 // test basic computations 159 { 160 DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); 161 DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); 162 DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); 163 DenseMatrix refM4 = DenseMatrix::Zero(rows, rows); 164 SparseMatrixType m1(rows, rows); 165 SparseMatrixType m2(rows, rows); 166 SparseMatrixType m3(rows, rows); 167 SparseMatrixType m4(rows, rows); 168 initSparse<Scalar>(density, refM1, m1); 169 initSparse<Scalar>(density, refM2, m2); 170 initSparse<Scalar>(density, refM3, m3); 171 initSparse<Scalar>(density, refM4, m4); 172 173 VERIFY_IS_APPROX(m1+m2, refM1+refM2); 174 VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); 175 VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2)); 176 VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); 177 178 VERIFY_IS_APPROX(m1*=s1, refM1*=s1); 179 VERIFY_IS_APPROX(m1/=s1, refM1/=s1); 180 181 VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); 182 VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); 183 184 if(SparseMatrixType::IsRowMajor) 185 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); 186 else 187 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); 188 189 VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate()); 190 VERIFY_IS_APPROX(m1.real(), refM1.real()); 191 192 refM4.setRandom(); 193 // sparse cwise* dense 194 VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4)); 195 // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); 196 } 197 198 // test transpose 199 { 200 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 201 SparseMatrixType m2(rows, rows); 202 initSparse<Scalar>(density, refMat2, m2); 203 VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); 204 VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); 205 206 VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); 207 } 208 209 // test innerVector() 210 { 211 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 212 SparseMatrixType m2(rows, rows); 213 initSparse<Scalar>(density, refMat2, m2); 214 int j0 = internal::random<int>(0,rows-1); 215 int j1 = internal::random<int>(0,rows-1); 216 if(SparseMatrixType::IsRowMajor) 217 VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0)); 218 else 219 VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); 220 221 if(SparseMatrixType::IsRowMajor) 222 VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.row(j0)+refMat2.row(j1)); 223 else 224 VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); 225 226 SparseMatrixType m3(rows,rows); 227 m3.reserve(VectorXi::Constant(rows,rows/2)); 228 for(int j=0; j<rows; ++j) 229 for(int k=0; k<j; ++k) 230 m3.insertByOuterInner(j,k) = k+1; 231 for(int j=0; j<rows; ++j) 232 { 233 VERIFY(j==internal::real(m3.innerVector(j).nonZeros())); 234 if(j>0) 235 VERIFY(j==internal::real(m3.innerVector(j).lastCoeff())); 236 } 237 m3.makeCompressed(); 238 for(int j=0; j<rows; ++j) 239 { 240 VERIFY(j==internal::real(m3.innerVector(j).nonZeros())); 241 if(j>0) 242 VERIFY(j==internal::real(m3.innerVector(j).lastCoeff())); 243 } 244 245 //m2.innerVector(j0) = 2*m2.innerVector(j1); 246 //refMat2.col(j0) = 2*refMat2.col(j1); 247 //VERIFY_IS_APPROX(m2, refMat2); 248 } 249 250 // test innerVectors() 251 { 252 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 253 SparseMatrixType m2(rows, rows); 254 initSparse<Scalar>(density, refMat2, m2); 255 int j0 = internal::random<int>(0,rows-2); 256 int j1 = internal::random<int>(0,rows-2); 257 int n0 = internal::random<int>(1,rows-(std::max)(j0,j1)); 258 if(SparseMatrixType::IsRowMajor) 259 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); 260 else 261 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); 262 if(SparseMatrixType::IsRowMajor) 263 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), 264 refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols)); 265 else 266 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), 267 refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); 268 //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); 269 //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0); 270 } 271 272 // test prune 273 { 274 SparseMatrixType m2(rows, rows); 275 DenseMatrix refM2(rows, rows); 276 refM2.setZero(); 277 int countFalseNonZero = 0; 278 int countTrueNonZero = 0; 279 for (int j=0; j<m2.outerSize(); ++j) 280 { 281 m2.startVec(j); 282 for (int i=0; i<m2.innerSize(); ++i) 283 { 284 float x = internal::random<float>(0,1); 285 if (x<0.1) 286 { 287 // do nothing 288 } 289 else if (x<0.5) 290 { 291 countFalseNonZero++; 292 m2.insertBackByOuterInner(j,i) = Scalar(0); 293 } 294 else 295 { 296 countTrueNonZero++; 297 m2.insertBackByOuterInner(j,i) = Scalar(1); 298 if(SparseMatrixType::IsRowMajor) 299 refM2(j,i) = Scalar(1); 300 else 301 refM2(i,j) = Scalar(1); 302 } 303 } 304 } 305 m2.finalize(); 306 VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros()); 307 VERIFY_IS_APPROX(m2, refM2); 308 m2.prune(Scalar(1)); 309 VERIFY(countTrueNonZero==m2.nonZeros()); 310 VERIFY_IS_APPROX(m2, refM2); 311 } 312 313 // test setFromTriplets 314 { 315 typedef Triplet<Scalar,Index> TripletType; 316 std::vector<TripletType> triplets; 317 int ntriplets = rows*cols; 318 triplets.reserve(ntriplets); 319 DenseMatrix refMat(rows,cols); 320 refMat.setZero(); 321 for(int i=0;i<ntriplets;++i) 322 { 323 int r = internal::random<int>(0,rows-1); 324 int c = internal::random<int>(0,cols-1); 325 Scalar v = internal::random<Scalar>(); 326 triplets.push_back(TripletType(r,c,v)); 327 refMat(r,c) += v; 328 } 329 SparseMatrixType m(rows,cols); 330 m.setFromTriplets(triplets.begin(), triplets.end()); 331 VERIFY_IS_APPROX(m, refMat); 332 } 333 334 // test triangularView 335 { 336 DenseMatrix refMat2(rows, rows), refMat3(rows, rows); 337 SparseMatrixType m2(rows, rows), m3(rows, rows); 338 initSparse<Scalar>(density, refMat2, m2); 339 refMat3 = refMat2.template triangularView<Lower>(); 340 m3 = m2.template triangularView<Lower>(); 341 VERIFY_IS_APPROX(m3, refMat3); 342 343 refMat3 = refMat2.template triangularView<Upper>(); 344 m3 = m2.template triangularView<Upper>(); 345 VERIFY_IS_APPROX(m3, refMat3); 346 347 refMat3 = refMat2.template triangularView<UnitUpper>(); 348 m3 = m2.template triangularView<UnitUpper>(); 349 VERIFY_IS_APPROX(m3, refMat3); 350 351 refMat3 = refMat2.template triangularView<UnitLower>(); 352 m3 = m2.template triangularView<UnitLower>(); 353 VERIFY_IS_APPROX(m3, refMat3); 354 } 355 356 // test selfadjointView 357 if(!SparseMatrixType::IsRowMajor) 358 { 359 DenseMatrix refMat2(rows, rows), refMat3(rows, rows); 360 SparseMatrixType m2(rows, rows), m3(rows, rows); 361 initSparse<Scalar>(density, refMat2, m2); 362 refMat3 = refMat2.template selfadjointView<Lower>(); 363 m3 = m2.template selfadjointView<Lower>(); 364 VERIFY_IS_APPROX(m3, refMat3); 365 } 366 367 // test sparseView 368 { 369 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 370 SparseMatrixType m2(rows, rows); 371 initSparse<Scalar>(density, refMat2, m2); 372 VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval()); 373 } 374 375 // test diagonal 376 { 377 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 378 SparseMatrixType m2(rows, rows); 379 initSparse<Scalar>(density, refMat2, m2); 380 VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval()); 381 } 382 } 383 384 void test_sparse_basic() 385 { 386 for(int i = 0; i < g_repeat; i++) { 387 int s = Eigen::internal::random<int>(1,50); 388 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) )); 389 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(s, s)) )); 390 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(s, s)) )); 391 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(s, s)) )); 392 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,long int>(s, s)) )); 393 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,long int>(s, s)) )); 394 } 395 } 396
