blob: 7b5f3eb385d6faa665305233b6868aa5e00ce288 [file] [log] [blame]
Gael Guennebaud86ccd992008-11-05 13:47:55 +00001// This file is part of Eigen, a lightweight C++ template library
Benoit Jacob6347b1d2009-05-22 20:25:33 +02002// for linear algebra.
Gael Guennebaud86ccd992008-11-05 13:47:55 +00003//
Gael Guennebaud22c76092011-03-22 11:58:22 +01004// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
Gael Guennebaud86ccd992008-11-05 13:47:55 +00005// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
Desire NUENTSA4cd82452013-06-11 14:42:29 +02006// Copyright (C) 2013 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
Gael Guennebaud86ccd992008-11-05 13:47:55 +00007//
Benoit Jacob69124cf2012-07-13 14:42:47 -04008// 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/.
Gael Guennebaud86ccd992008-11-05 13:47:55 +000011
Gael Guennebaudc43154b2015-03-04 10:16:46 +010012static long g_realloc_count = 0;
13#define EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN g_realloc_count++;
14
Gael Guennebaud86ccd992008-11-05 13:47:55 +000015#include "sparse.h"
16
Gael Guennebaud178858f2009-01-19 15:20:45 +000017template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
18{
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +010019 typedef typename SparseMatrixType::StorageIndex StorageIndex;
20 typedef Matrix<StorageIndex,2,1> Vector2;
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +020021
Gael Guennebaudfc202ba2015-02-13 18:57:41 +010022 const Index rows = ref.rows();
23 const Index cols = ref.cols();
Gael Guennebaud0a537cb2016-02-12 15:58:31 +010024 //const Index inner = ref.innerSize();
25 //const Index outer = ref.outerSize();
Christoph Hertzberg0833b822014-10-31 17:12:13 +010026
Gael Guennebaud178858f2009-01-19 15:20:45 +000027 typedef typename SparseMatrixType::Scalar Scalar;
28 enum { Flags = SparseMatrixType::Flags };
Gael Guennebaud9f795582009-12-16 19:18:40 +010029
Gael Guennebaud42e25782011-08-19 14:18:05 +020030 double density = (std::max)(8./(rows*cols), 0.01);
Gael Guennebaud86ccd992008-11-05 13:47:55 +000031 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
32 typedef Matrix<Scalar,Dynamic,1> DenseVector;
33 Scalar eps = 1e-6;
34
Benoit Jacob47160402010-10-25 10:15:22 -040035 Scalar s1 = internal::random<Scalar>();
Gael Guennebaud86ccd992008-11-05 13:47:55 +000036 {
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020037 SparseMatrixType m(rows, cols);
38 DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
39 DenseVector vec1 = DenseVector::Random(rows);
Gael Guennebaud86ccd992008-11-05 13:47:55 +000040
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +020041 std::vector<Vector2> zeroCoords;
42 std::vector<Vector2> nonzeroCoords;
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020043 initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
Gael Guennebaud86ccd992008-11-05 13:47:55 +000044
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020045 // test coeff and coeffRef
Christoph Hertzberg0833b822014-10-31 17:12:13 +010046 for (std::size_t i=0; i<zeroCoords.size(); ++i)
Gael Guennebaud86ccd992008-11-05 13:47:55 +000047 {
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020048 VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
49 if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value)
Christoph Hertzberg0833b822014-10-31 17:12:13 +010050 VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[i].x(),zeroCoords[i].y()) = 5 );
Gael Guennebaud86ccd992008-11-05 13:47:55 +000051 }
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020052 VERIFY_IS_APPROX(m, refMat);
Gael Guennebaud86ccd992008-11-05 13:47:55 +000053
Christoph Hertzberg0833b822014-10-31 17:12:13 +010054 if(!nonzeroCoords.empty()) {
55 m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
56 refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
57 }
Gael Guennebaud86ccd992008-11-05 13:47:55 +000058
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020059 VERIFY_IS_APPROX(m, refMat);
Christoph Hertzberg0833b822014-10-31 17:12:13 +010060
Gael Guennebauda915f022013-06-28 16:16:02 +020061 // test assertion
62 VERIFY_RAISES_ASSERT( m.coeffRef(-1,1) = 0 );
63 VERIFY_RAISES_ASSERT( m.coeffRef(0,m.cols()) = 0 );
Gael Guennebaudd1d7a1a2013-06-23 19:11:32 +020064 }
Gael Guennebaud86ccd992008-11-05 13:47:55 +000065
Gael Guennebaud28293142009-05-04 14:25:12 +000066 // test insert (inner random)
Gael Guennebaud5015e482008-12-11 18:26:24 +000067 {
68 DenseMatrix m1(rows,cols);
69 m1.setZero();
Gael Guennebaud178858f2009-01-19 15:20:45 +000070 SparseMatrixType m2(rows,cols);
Gael Guennebaudc43154b2015-03-04 10:16:46 +010071 bool call_reserve = internal::random<int>()%2;
72 Index nnz = internal::random<int>(1,int(rows)/2);
73 if(call_reserve)
74 {
75 if(internal::random<int>()%2)
76 m2.reserve(VectorXi::Constant(m2.outerSize(), int(nnz)));
77 else
78 m2.reserve(m2.outerSize() * nnz);
79 }
80 g_realloc_count = 0;
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +020081 for (Index j=0; j<cols; ++j)
Gael Guennebaud5015e482008-12-11 18:26:24 +000082 {
Gael Guennebaudc43154b2015-03-04 10:16:46 +010083 for (Index k=0; k<nnz; ++k)
Gael Guennebaud5015e482008-12-11 18:26:24 +000084 {
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +020085 Index i = internal::random<Index>(0,rows-1);
Gael Guennebaud5015e482008-12-11 18:26:24 +000086 if (m1.coeff(i,j)==Scalar(0))
Benoit Jacob47160402010-10-25 10:15:22 -040087 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
Gael Guennebaud5015e482008-12-11 18:26:24 +000088 }
89 }
Gael Guennebaudc43154b2015-03-04 10:16:46 +010090
91 if(call_reserve && !SparseMatrixType::IsRowMajor)
92 {
93 VERIFY(g_realloc_count==0);
94 }
95
Gael Guennebaud28293142009-05-04 14:25:12 +000096 m2.finalize();
97 VERIFY_IS_APPROX(m2,m1);
98 }
Gael Guennebaud9f795582009-12-16 19:18:40 +010099
Gael Guennebaud28293142009-05-04 14:25:12 +0000100 // test insert (fully random)
101 {
102 DenseMatrix m1(rows,cols);
103 m1.setZero();
104 SparseMatrixType m2(rows,cols);
Gael Guennebaud4ca89f32011-12-02 19:00:16 +0100105 if(internal::random<int>()%2)
106 m2.reserve(VectorXi::Constant(m2.outerSize(), 2));
Gael Guennebaud28293142009-05-04 14:25:12 +0000107 for (int k=0; k<rows*cols; ++k)
108 {
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +0200109 Index i = internal::random<Index>(0,rows-1);
110 Index j = internal::random<Index>(0,cols-1);
Gael Guennebaud4ca89f32011-12-02 19:00:16 +0100111 if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2))
Benoit Jacob47160402010-10-25 10:15:22 -0400112 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
Gael Guennebaud4ca89f32011-12-02 19:00:16 +0100113 else
114 {
115 Scalar v = internal::random<Scalar>();
116 m2.coeffRef(i,j) += v;
117 m1(i,j) += v;
118 }
Gael Guennebaud28293142009-05-04 14:25:12 +0000119 }
Gael Guennebaudc4c70662009-01-14 14:24:10 +0000120 VERIFY_IS_APPROX(m2,m1);
Gael Guennebaud5015e482008-12-11 18:26:24 +0000121 }
Gael Guennebaud7706baf2011-09-08 13:42:54 +0200122
123 // test insert (un-compressed)
124 for(int mode=0;mode<4;++mode)
125 {
126 DenseMatrix m1(rows,cols);
127 m1.setZero();
128 SparseMatrixType m2(rows,cols);
Gael Guennebaudb47ef142014-07-08 16:47:11 +0200129 VectorXi r(VectorXi::Constant(m2.outerSize(), ((mode%2)==0) ? int(m2.innerSize()) : std::max<int>(1,int(m2.innerSize())/8)));
Gael Guennebaud7706baf2011-09-08 13:42:54 +0200130 m2.reserve(r);
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100131 for (Index k=0; k<rows*cols; ++k)
Gael Guennebaud7706baf2011-09-08 13:42:54 +0200132 {
Gael Guennebaud6d1f5db2013-07-10 23:48:26 +0200133 Index i = internal::random<Index>(0,rows-1);
134 Index j = internal::random<Index>(0,cols-1);
Gael Guennebaud7706baf2011-09-08 13:42:54 +0200135 if (m1.coeff(i,j)==Scalar(0))
136 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
137 if(mode==3)
138 m2.reserve(r);
139 }
Gael Guennebaud4ca89f32011-12-02 19:00:16 +0100140 if(internal::random<int>()%2)
141 m2.makeCompressed();
Gael Guennebaud7706baf2011-09-08 13:42:54 +0200142 VERIFY_IS_APPROX(m2,m1);
143 }
Gael Guennebaud9f795582009-12-16 19:18:40 +0100144
Gael Guennebaud4e602832012-11-16 09:02:50 +0100145 // test basic computations
146 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100147 DenseMatrix refM1 = DenseMatrix::Zero(rows, cols);
148 DenseMatrix refM2 = DenseMatrix::Zero(rows, cols);
149 DenseMatrix refM3 = DenseMatrix::Zero(rows, cols);
150 DenseMatrix refM4 = DenseMatrix::Zero(rows, cols);
151 SparseMatrixType m1(rows, cols);
152 SparseMatrixType m2(rows, cols);
153 SparseMatrixType m3(rows, cols);
154 SparseMatrixType m4(rows, cols);
Gael Guennebaud4e602832012-11-16 09:02:50 +0100155 initSparse<Scalar>(density, refM1, m1);
156 initSparse<Scalar>(density, refM2, m2);
157 initSparse<Scalar>(density, refM3, m3);
158 initSparse<Scalar>(density, refM4, m4);
159
Gael Guennebaud2c1b56f2016-05-31 10:56:53 +0200160 if(internal::random<bool>())
161 m1.makeCompressed();
162
Gael Guennebaud4aac8722014-07-22 12:54:03 +0200163 VERIFY_IS_APPROX(m1*s1, refM1*s1);
Gael Guennebaud4e602832012-11-16 09:02:50 +0100164 VERIFY_IS_APPROX(m1+m2, refM1+refM2);
165 VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
166 VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2));
167 VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
168
Gael Guennebaud4e602832012-11-16 09:02:50 +0100169 if(SparseMatrixType::IsRowMajor)
170 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0)));
171 else
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100172 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.col(0)), refM1.col(0).dot(refM2.col(0)));
Gael Guennebaud3573a102014-02-17 13:46:17 +0100173
174 DenseVector rv = DenseVector::Random(m1.cols());
175 DenseVector cv = DenseVector::Random(m1.rows());
176 Index r = internal::random<Index>(0,m1.rows()-2);
177 Index c = internal::random<Index>(0,m1.cols()-1);
178 VERIFY_IS_APPROX(( m1.template block<1,Dynamic>(r,0,1,m1.cols()).dot(rv)) , refM1.row(r).dot(rv));
179 VERIFY_IS_APPROX(m1.row(r).dot(rv), refM1.row(r).dot(rv));
180 VERIFY_IS_APPROX(m1.col(c).dot(cv), refM1.col(c).dot(cv));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100181
182 VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate());
183 VERIFY_IS_APPROX(m1.real(), refM1.real());
184
185 refM4.setRandom();
186 // sparse cwise* dense
187 VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4));
Gael Guennebaud90275082015-11-04 17:42:07 +0100188 // dense cwise* sparse
189 VERIFY_IS_APPROX(refM4.cwiseProduct(m3), refM4.cwiseProduct(refM3));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100190// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
191
Gael Guennebaud15084cf2016-01-29 22:09:45 +0100192 VERIFY_IS_APPROX(refM4 + m3, refM4 + refM3);
193 VERIFY_IS_APPROX(m3 + refM4, refM3 + refM4);
194 VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3);
195 VERIFY_IS_APPROX(m3 - refM4, refM3 - refM4);
196
Gael Guennebaud2c1b56f2016-05-31 10:56:53 +0200197 VERIFY_IS_APPROX(m1.sum(), refM1.sum());
198
199 VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
200 VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
201
202 VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
203 VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
204
Gael Guennebaud4e602832012-11-16 09:02:50 +0100205 // test aliasing
206 VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1));
207 VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval()));
208 VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval()));
209 VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1));
Gael Guennebaud7e029d12016-08-29 12:06:37 +0200210
211 if(m1.isCompressed())
212 {
213 VERIFY_IS_APPROX(m1.coeffs().sum(), m1.sum());
214 m1.coeffs() += s1;
215 for(Index j = 0; j<m1.outerSize(); ++j)
216 for(typename SparseMatrixType::InnerIterator it(m1,j); it; ++it)
217 refM1(it.row(), it.col()) += s1;
218 VERIFY_IS_APPROX(m1, refM1);
219 }
Gael Guennebaud4e602832012-11-16 09:02:50 +0100220 }
221
222 // test transpose
223 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100224 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
225 SparseMatrixType m2(rows, cols);
Gael Guennebaud4e602832012-11-16 09:02:50 +0100226 initSparse<Scalar>(density, refMat2, m2);
227 VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
228 VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
229
230 VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
Gael Guennebaudff46ec02014-09-22 23:33:28 +0200231
232 // check isApprox handles opposite storage order
233 typename Transpose<SparseMatrixType>::PlainObject m3(m2);
234 VERIFY(m2.isApprox(m3));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100235 }
236
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000237 // test prune
238 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100239 SparseMatrixType m2(rows, cols);
240 DenseMatrix refM2(rows, cols);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000241 refM2.setZero();
242 int countFalseNonZero = 0;
243 int countTrueNonZero = 0;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200244 m2.reserve(VectorXi::Constant(m2.outerSize(), int(m2.innerSize())));
245 for (Index j=0; j<m2.cols(); ++j)
Gael Guennebaud28293142009-05-04 14:25:12 +0000246 {
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200247 for (Index i=0; i<m2.rows(); ++i)
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000248 {
Benoit Jacob47160402010-10-25 10:15:22 -0400249 float x = internal::random<float>(0,1);
Christoph Hertzbergdacb4692016-05-05 13:35:45 +0200250 if (x<0.1f)
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000251 {
252 // do nothing
253 }
Christoph Hertzbergdacb4692016-05-05 13:35:45 +0200254 else if (x<0.5f)
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000255 {
256 countFalseNonZero++;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200257 m2.insert(i,j) = Scalar(0);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000258 }
259 else
260 {
261 countTrueNonZero++;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200262 m2.insert(i,j) = Scalar(1);
263 refM2(i,j) = Scalar(1);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000264 }
265 }
Gael Guennebaud28293142009-05-04 14:25:12 +0000266 }
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200267 if(internal::random<bool>())
268 m2.makeCompressed();
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000269 VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200270 if(countTrueNonZero>0)
271 VERIFY_IS_APPROX(m2, refM2);
Hauke Heibeld204ec42010-11-02 14:33:33 +0100272 m2.prune(Scalar(1));
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000273 VERIFY(countTrueNonZero==m2.nonZeros());
274 VERIFY_IS_APPROX(m2, refM2);
275 }
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100276
Gael Guennebaud87138072012-01-28 11:13:59 +0100277 // test setFromTriplets
278 {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100279 typedef Triplet<Scalar,StorageIndex> TripletType;
Gael Guennebaud87138072012-01-28 11:13:59 +0100280 std::vector<TripletType> triplets;
Gael Guennebaudb47ef142014-07-08 16:47:11 +0200281 Index ntriplets = rows*cols;
Gael Guennebaud87138072012-01-28 11:13:59 +0100282 triplets.reserve(ntriplets);
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200283 DenseMatrix refMat_sum = DenseMatrix::Zero(rows,cols);
284 DenseMatrix refMat_prod = DenseMatrix::Zero(rows,cols);
285 DenseMatrix refMat_last = DenseMatrix::Zero(rows,cols);
286
Gael Guennebaudb47ef142014-07-08 16:47:11 +0200287 for(Index i=0;i<ntriplets;++i)
Gael Guennebaud87138072012-01-28 11:13:59 +0100288 {
Gael Guennebaudaa6c5162015-02-16 13:19:05 +0100289 StorageIndex r = internal::random<StorageIndex>(0,StorageIndex(rows-1));
290 StorageIndex c = internal::random<StorageIndex>(0,StorageIndex(cols-1));
Gael Guennebaud87138072012-01-28 11:13:59 +0100291 Scalar v = internal::random<Scalar>();
Gael Guennebaud18e3ac02012-01-31 09:14:01 +0100292 triplets.push_back(TripletType(r,c,v));
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200293 refMat_sum(r,c) += v;
294 if(std::abs(refMat_prod(r,c))==0)
295 refMat_prod(r,c) = v;
296 else
297 refMat_prod(r,c) *= v;
298 refMat_last(r,c) = v;
Gael Guennebaud87138072012-01-28 11:13:59 +0100299 }
300 SparseMatrixType m(rows,cols);
301 m.setFromTriplets(triplets.begin(), triplets.end());
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200302 VERIFY_IS_APPROX(m, refMat_sum);
303
304 m.setFromTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>());
305 VERIFY_IS_APPROX(m, refMat_prod);
306#if (defined(__cplusplus) && __cplusplus >= 201103L)
307 m.setFromTriplets(triplets.begin(), triplets.end(), [] (Scalar,Scalar b) { return b; });
308 VERIFY_IS_APPROX(m, refMat_last);
309#endif
Gael Guennebaud87138072012-01-28 11:13:59 +0100310 }
Gael Guennebaud3af29ca2015-02-09 10:23:45 +0100311
312 // test Map
313 {
314 DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
315 SparseMatrixType m2(rows, cols), m3(rows, cols);
316 initSparse<Scalar>(density, refMat2, m2);
317 initSparse<Scalar>(density, refMat3, m3);
318 {
319 Map<SparseMatrixType> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
320 Map<SparseMatrixType> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr());
321 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
322 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
323 }
324 {
Gael Guennebaudfe513192015-02-13 10:03:53 +0100325 MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
326 MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr());
Gael Guennebaud3af29ca2015-02-09 10:23:45 +0100327 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
328 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
329 }
Gael Guennebaud757971e2016-07-26 09:40:19 +0200330
331 Index i = internal::random<Index>(0,rows-1);
332 Index j = internal::random<Index>(0,cols-1);
333 m2.coeffRef(i,j) = 123;
334 if(internal::random<bool>())
335 m2.makeCompressed();
336 Map<SparseMatrixType> mapMat2(rows, cols, m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
337 VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(123));
338 VERIFY_IS_EQUAL(mapMat2.coeff(i,j),Scalar(123));
339 mapMat2.coeffRef(i,j) = -123;
340 VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(-123));
Gael Guennebaud3af29ca2015-02-09 10:23:45 +0100341 }
Gael Guennebaud87138072012-01-28 11:13:59 +0100342
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100343 // test triangularView
344 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100345 DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
346 SparseMatrixType m2(rows, cols), m3(rows, cols);
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100347 initSparse<Scalar>(density, refMat2, m2);
348 refMat3 = refMat2.template triangularView<Lower>();
349 m3 = m2.template triangularView<Lower>();
350 VERIFY_IS_APPROX(m3, refMat3);
351
352 refMat3 = refMat2.template triangularView<Upper>();
353 m3 = m2.template triangularView<Upper>();
354 VERIFY_IS_APPROX(m3, refMat3);
355
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100356 {
357 refMat3 = refMat2.template triangularView<UnitUpper>();
358 m3 = m2.template triangularView<UnitUpper>();
359 VERIFY_IS_APPROX(m3, refMat3);
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100360
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100361 refMat3 = refMat2.template triangularView<UnitLower>();
362 m3 = m2.template triangularView<UnitLower>();
363 VERIFY_IS_APPROX(m3, refMat3);
364 }
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200365
366 refMat3 = refMat2.template triangularView<StrictlyUpper>();
367 m3 = m2.template triangularView<StrictlyUpper>();
368 VERIFY_IS_APPROX(m3, refMat3);
369
370 refMat3 = refMat2.template triangularView<StrictlyLower>();
371 m3 = m2.template triangularView<StrictlyLower>();
372 VERIFY_IS_APPROX(m3, refMat3);
Gael Guennebaudf6b1dee2015-11-04 17:02:32 +0100373
374 // check sparse-traingular to dense
375 refMat3 = m2.template triangularView<StrictlyUpper>();
376 VERIFY_IS_APPROX(refMat3, DenseMatrix(refMat2.template triangularView<StrictlyUpper>()));
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100377 }
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200378
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100379 // test selfadjointView
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100380 if(!SparseMatrixType::IsRowMajor)
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100381 {
382 DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
383 SparseMatrixType m2(rows, rows), m3(rows, rows);
384 initSparse<Scalar>(density, refMat2, m2);
385 refMat3 = refMat2.template selfadjointView<Lower>();
386 m3 = m2.template selfadjointView<Lower>();
387 VERIFY_IS_APPROX(m3, refMat3);
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100388
389 // selfadjointView only works for square matrices:
390 SparseMatrixType m4(rows, rows+1);
391 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Lower>());
392 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Upper>());
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100393 }
394
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200395 // test sparseView
396 {
397 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
398 SparseMatrixType m2(rows, rows);
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100399 initSparse<Scalar>(density, refMat2, m2);
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200400 VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
Gael Guennebaud8456bbb2016-05-18 16:53:28 +0200401
402 // sparse view on expressions:
403 VERIFY_IS_APPROX((s1*m2).eval(), (s1*refMat2).sparseView().eval());
404 VERIFY_IS_APPROX((m2+m2).eval(), (refMat2+refMat2).sparseView().eval());
405 VERIFY_IS_APPROX((m2*m2).eval(), (refMat2.lazyProduct(refMat2)).sparseView().eval());
406 VERIFY_IS_APPROX((m2*m2).eval(), (refMat2*refMat2).sparseView().eval());
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200407 }
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100408
409 // test diagonal
410 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100411 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
412 SparseMatrixType m2(rows, cols);
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100413 initSparse<Scalar>(density, refMat2, m2);
414 VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval());
Gael Guennebaudb26e6972014-12-01 14:41:39 +0100415 VERIFY_IS_APPROX(const_cast<const SparseMatrixType&>(m2).diagonal(), refMat2.diagonal().eval());
416
417 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag);
418 m2.diagonal() += refMat2.diagonal();
419 refMat2.diagonal() += refMat2.diagonal();
420 VERIFY_IS_APPROX(m2, refMat2);
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100421 }
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200422
Gael Guennebaud62f21e22015-06-24 17:55:00 +0200423 // test diagonal to sparse
424 {
425 DenseVector d = DenseVector::Random(rows);
426 DenseMatrix refMat2 = d.asDiagonal();
427 SparseMatrixType m2(rows, rows);
428 m2 = d.asDiagonal();
429 VERIFY_IS_APPROX(m2, refMat2);
Gael Guennebaud4c8cd132015-06-24 18:11:06 +0200430 SparseMatrixType m3(d.asDiagonal());
431 VERIFY_IS_APPROX(m3, refMat2);
Gael Guennebaud62f21e22015-06-24 17:55:00 +0200432 refMat2 += d.asDiagonal();
433 m2 += d.asDiagonal();
434 VERIFY_IS_APPROX(m2, refMat2);
435 }
436
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200437 // test conservative resize
438 {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100439 std::vector< std::pair<StorageIndex,StorageIndex> > inc;
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100440 if(rows > 3 && cols > 2)
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100441 inc.push_back(std::pair<StorageIndex,StorageIndex>(-3,-2));
442 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,0));
443 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,2));
444 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,0));
445 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,3));
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200446
447 for(size_t i = 0; i< inc.size(); i++) {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100448 StorageIndex incRows = inc[i].first;
449 StorageIndex incCols = inc[i].second;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200450 SparseMatrixType m1(rows, cols);
451 DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols);
452 initSparse<Scalar>(density, refMat1, m1);
453
454 m1.conservativeResize(rows+incRows, cols+incCols);
455 refMat1.conservativeResize(rows+incRows, cols+incCols);
456 if (incRows > 0) refMat1.bottomRows(incRows).setZero();
457 if (incCols > 0) refMat1.rightCols(incCols).setZero();
458
459 VERIFY_IS_APPROX(m1, refMat1);
460
461 // Insert new values
462 if (incRows > 0)
Gael Guennebaud7ee378d2013-07-12 16:40:02 +0200463 m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200464 if (incCols > 0)
Gael Guennebaud7ee378d2013-07-12 16:40:02 +0200465 m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200466
467 VERIFY_IS_APPROX(m1, refMat1);
468
469
470 }
471 }
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200472
473 // test Identity matrix
474 {
475 DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows);
476 SparseMatrixType m1(rows, rows);
477 m1.setIdentity();
478 VERIFY_IS_APPROX(m1, refMat1);
Gael Guennebaud8a211bb2015-10-25 22:01:58 +0100479 for(int k=0; k<rows*rows/4; ++k)
480 {
481 Index i = internal::random<Index>(0,rows-1);
482 Index j = internal::random<Index>(0,rows-1);
Gael Guennebaud73f692d2015-10-27 11:01:37 +0100483 Scalar v = internal::random<Scalar>();
Gael Guennebaud8a211bb2015-10-25 22:01:58 +0100484 m1.coeffRef(i,j) = v;
485 refMat1.coeffRef(i,j) = v;
486 VERIFY_IS_APPROX(m1, refMat1);
487 if(internal::random<Index>(0,10)<2)
488 m1.makeCompressed();
489 }
490 m1.setIdentity();
491 refMat1.setIdentity();
492 VERIFY_IS_APPROX(m1, refMat1);
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200493 }
Gael Guennebaudec469702016-02-01 15:04:33 +0100494
495 // test array/vector of InnerIterator
496 {
497 typedef typename SparseMatrixType::InnerIterator IteratorType;
498
499 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
500 SparseMatrixType m2(rows, cols);
501 initSparse<Scalar>(density, refMat2, m2);
502 IteratorType static_array[2];
503 static_array[0] = IteratorType(m2,0);
504 static_array[1] = IteratorType(m2,m2.outerSize()-1);
505 VERIFY( static_array[0] || m2.innerVector(static_array[0].outer()).nonZeros() == 0 );
506 VERIFY( static_array[1] || m2.innerVector(static_array[1].outer()).nonZeros() == 0 );
507 if(static_array[0] && static_array[1])
508 {
509 ++(static_array[1]);
510 static_array[1] = IteratorType(m2,0);
511 VERIFY( static_array[1] );
512 VERIFY( static_array[1].index() == static_array[0].index() );
513 VERIFY( static_array[1].outer() == static_array[0].outer() );
514 VERIFY( static_array[1].value() == static_array[0].value() );
515 }
516
517 std::vector<IteratorType> iters(2);
518 iters[0] = IteratorType(m2,0);
519 iters[1] = IteratorType(m2,m2.outerSize()-1);
520 }
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000521}
522
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100523
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100524template<typename SparseMatrixType>
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100525void big_sparse_triplet(Index rows, Index cols, double density) {
526 typedef typename SparseMatrixType::StorageIndex StorageIndex;
527 typedef typename SparseMatrixType::Scalar Scalar;
528 typedef Triplet<Scalar,Index> TripletType;
529 std::vector<TripletType> triplets;
530 double nelements = density * rows*cols;
531 VERIFY(nelements>=0 && nelements < NumTraits<StorageIndex>::highest());
532 Index ntriplets = Index(nelements);
533 triplets.reserve(ntriplets);
534 Scalar sum = Scalar(0);
535 for(Index i=0;i<ntriplets;++i)
536 {
537 Index r = internal::random<Index>(0,rows-1);
538 Index c = internal::random<Index>(0,cols-1);
539 Scalar v = internal::random<Scalar>();
540 triplets.push_back(TripletType(r,c,v));
541 sum += v;
542 }
543 SparseMatrixType m(rows,cols);
544 m.setFromTriplets(triplets.begin(), triplets.end());
545 VERIFY(m.nonZeros() <= ntriplets);
546 VERIFY_IS_APPROX(sum, m.sum());
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100547}
548
549
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000550void test_sparse_basic()
551{
552 for(int i = 0; i < g_repeat; i++) {
Gael Guennebaudc43154b2015-03-04 10:16:46 +0100553 int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200);
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100554 if(Eigen::internal::random<int>(0,4) == 0) {
555 r = c; // check square matrices in 25% of tries
556 }
557 EIGEN_UNUSED_VARIABLE(r+c);
558 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(1, 1)) ));
Gael Guennebaud91fe1502011-06-07 11:28:16 +0200559 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) ));
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100560 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(r, c)) ));
561 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(r, c)) ));
562 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(r, c)) ));
Gael Guennebaudd7698c12015-03-19 15:11:05 +0100563 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,ColMajor,long int>(r, c)) ));
564 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,RowMajor,long int>(r, c)) ));
Gael Guennebaud47e89022013-06-10 10:34:03 +0200565
Gael Guennebaudc43154b2015-03-04 10:16:46 +0100566 r = Eigen::internal::random<int>(1,100);
567 c = Eigen::internal::random<int>(1,100);
568 if(Eigen::internal::random<int>(0,4) == 0) {
569 r = c; // check square matrices in 25% of tries
570 }
571
Gael Guennebaudd7698c12015-03-19 15:11:05 +0100572 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) ));
573 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) ));
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000574 }
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100575
576 // Regression test for bug 900: (manually insert higher values here, if you have enough RAM):
577 CALL_SUBTEST_3((big_sparse_triplet<SparseMatrix<float, RowMajor, int> >(10000, 10000, 0.125)));
578 CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125)));
Gael Guennebaudbfd6ee62015-11-06 15:05:37 +0100579
580 // Regression test for bug 1105
Christoph Hertzberg7268b102016-05-11 19:41:53 +0200581#ifdef EIGEN_TEST_PART_7
Gael Guennebaudbfd6ee62015-11-06 15:05:37 +0100582 {
583 int n = Eigen::internal::random<int>(200,600);
584 SparseMatrix<std::complex<double>,0, long> mat(n, n);
585 std::complex<double> val;
586
587 for(int i=0; i<n; ++i)
588 {
589 mat.coeffRef(i, i%(n/10)) = val;
590 VERIFY(mat.data().allocatedSize()<20*n);
591 }
592 }
593#endif
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000594}