blob: 0a06c828bd5cf82b88f6f275466156daf87c378e [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();
Christoph Hertzberg0833b822014-10-31 17:12:13 +010024 const Index inner = ref.innerSize();
25 const Index outer = ref.outerSize();
26
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 Guennebaud4aac8722014-07-22 12:54:03 +0200160 VERIFY_IS_APPROX(m1*s1, refM1*s1);
Gael Guennebaud4e602832012-11-16 09:02:50 +0100161 VERIFY_IS_APPROX(m1+m2, refM1+refM2);
162 VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
163 VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2));
164 VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
165
166 VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
167 VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
168
169 VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
170 VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
171
172 if(SparseMatrixType::IsRowMajor)
173 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0)));
174 else
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100175 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.col(0)), refM1.col(0).dot(refM2.col(0)));
Gael Guennebaud3573a102014-02-17 13:46:17 +0100176
177 DenseVector rv = DenseVector::Random(m1.cols());
178 DenseVector cv = DenseVector::Random(m1.rows());
179 Index r = internal::random<Index>(0,m1.rows()-2);
180 Index c = internal::random<Index>(0,m1.cols()-1);
181 VERIFY_IS_APPROX(( m1.template block<1,Dynamic>(r,0,1,m1.cols()).dot(rv)) , refM1.row(r).dot(rv));
182 VERIFY_IS_APPROX(m1.row(r).dot(rv), refM1.row(r).dot(rv));
183 VERIFY_IS_APPROX(m1.col(c).dot(cv), refM1.col(c).dot(cv));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100184
185 VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate());
186 VERIFY_IS_APPROX(m1.real(), refM1.real());
187
188 refM4.setRandom();
189 // sparse cwise* dense
190 VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4));
Gael Guennebaud90275082015-11-04 17:42:07 +0100191 // dense cwise* sparse
192 VERIFY_IS_APPROX(refM4.cwiseProduct(m3), refM4.cwiseProduct(refM3));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100193// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
194
Gael Guennebaud15084cf2016-01-29 22:09:45 +0100195 VERIFY_IS_APPROX(refM4 + m3, refM4 + refM3);
196 VERIFY_IS_APPROX(m3 + refM4, refM3 + refM4);
197 VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3);
198 VERIFY_IS_APPROX(m3 - refM4, refM3 - refM4);
199
Gael Guennebaud4e602832012-11-16 09:02:50 +0100200 // test aliasing
201 VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1));
202 VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval()));
203 VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval()));
204 VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1));
205 }
206
207 // test transpose
208 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100209 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
210 SparseMatrixType m2(rows, cols);
Gael Guennebaud4e602832012-11-16 09:02:50 +0100211 initSparse<Scalar>(density, refMat2, m2);
212 VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
213 VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
214
215 VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
Gael Guennebaudff46ec02014-09-22 23:33:28 +0200216
217 // check isApprox handles opposite storage order
218 typename Transpose<SparseMatrixType>::PlainObject m3(m2);
219 VERIFY(m2.isApprox(m3));
Gael Guennebaud4e602832012-11-16 09:02:50 +0100220 }
221
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000222 // test prune
223 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100224 SparseMatrixType m2(rows, cols);
225 DenseMatrix refM2(rows, cols);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000226 refM2.setZero();
227 int countFalseNonZero = 0;
228 int countTrueNonZero = 0;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200229 m2.reserve(VectorXi::Constant(m2.outerSize(), int(m2.innerSize())));
230 for (Index j=0; j<m2.cols(); ++j)
Gael Guennebaud28293142009-05-04 14:25:12 +0000231 {
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200232 for (Index i=0; i<m2.rows(); ++i)
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000233 {
Benoit Jacob47160402010-10-25 10:15:22 -0400234 float x = internal::random<float>(0,1);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000235 if (x<0.1)
236 {
237 // do nothing
238 }
239 else if (x<0.5)
240 {
241 countFalseNonZero++;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200242 m2.insert(i,j) = Scalar(0);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000243 }
244 else
245 {
246 countTrueNonZero++;
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200247 m2.insert(i,j) = Scalar(1);
248 refM2(i,j) = Scalar(1);
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000249 }
250 }
Gael Guennebaud28293142009-05-04 14:25:12 +0000251 }
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200252 if(internal::random<bool>())
253 m2.makeCompressed();
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000254 VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
Gael Guennebauda44d91a2015-10-13 10:53:38 +0200255 if(countTrueNonZero>0)
256 VERIFY_IS_APPROX(m2, refM2);
Hauke Heibeld204ec42010-11-02 14:33:33 +0100257 m2.prune(Scalar(1));
Gael Guennebaud52cf07d2009-01-21 18:46:04 +0000258 VERIFY(countTrueNonZero==m2.nonZeros());
259 VERIFY_IS_APPROX(m2, refM2);
260 }
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100261
Gael Guennebaud87138072012-01-28 11:13:59 +0100262 // test setFromTriplets
263 {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100264 typedef Triplet<Scalar,StorageIndex> TripletType;
Gael Guennebaud87138072012-01-28 11:13:59 +0100265 std::vector<TripletType> triplets;
Gael Guennebaudb47ef142014-07-08 16:47:11 +0200266 Index ntriplets = rows*cols;
Gael Guennebaud87138072012-01-28 11:13:59 +0100267 triplets.reserve(ntriplets);
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200268 DenseMatrix refMat_sum = DenseMatrix::Zero(rows,cols);
269 DenseMatrix refMat_prod = DenseMatrix::Zero(rows,cols);
270 DenseMatrix refMat_last = DenseMatrix::Zero(rows,cols);
271
Gael Guennebaudb47ef142014-07-08 16:47:11 +0200272 for(Index i=0;i<ntriplets;++i)
Gael Guennebaud87138072012-01-28 11:13:59 +0100273 {
Gael Guennebaudaa6c5162015-02-16 13:19:05 +0100274 StorageIndex r = internal::random<StorageIndex>(0,StorageIndex(rows-1));
275 StorageIndex c = internal::random<StorageIndex>(0,StorageIndex(cols-1));
Gael Guennebaud87138072012-01-28 11:13:59 +0100276 Scalar v = internal::random<Scalar>();
Gael Guennebaud18e3ac02012-01-31 09:14:01 +0100277 triplets.push_back(TripletType(r,c,v));
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200278 refMat_sum(r,c) += v;
279 if(std::abs(refMat_prod(r,c))==0)
280 refMat_prod(r,c) = v;
281 else
282 refMat_prod(r,c) *= v;
283 refMat_last(r,c) = v;
Gael Guennebaud87138072012-01-28 11:13:59 +0100284 }
285 SparseMatrixType m(rows,cols);
286 m.setFromTriplets(triplets.begin(), triplets.end());
Gael Guennebaudb4c79ee2015-10-13 11:30:41 +0200287 VERIFY_IS_APPROX(m, refMat_sum);
288
289 m.setFromTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>());
290 VERIFY_IS_APPROX(m, refMat_prod);
291#if (defined(__cplusplus) && __cplusplus >= 201103L)
292 m.setFromTriplets(triplets.begin(), triplets.end(), [] (Scalar,Scalar b) { return b; });
293 VERIFY_IS_APPROX(m, refMat_last);
294#endif
Gael Guennebaud87138072012-01-28 11:13:59 +0100295 }
Gael Guennebaud3af29ca2015-02-09 10:23:45 +0100296
297 // test Map
298 {
299 DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
300 SparseMatrixType m2(rows, cols), m3(rows, cols);
301 initSparse<Scalar>(density, refMat2, m2);
302 initSparse<Scalar>(density, refMat3, m3);
303 {
304 Map<SparseMatrixType> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
305 Map<SparseMatrixType> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr());
306 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
307 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
308 }
309 {
Gael Guennebaudfe513192015-02-13 10:03:53 +0100310 MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
311 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 +0100312 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
313 VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
314 }
315 }
Gael Guennebaud87138072012-01-28 11:13:59 +0100316
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100317 // test triangularView
318 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100319 DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
320 SparseMatrixType m2(rows, cols), m3(rows, cols);
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100321 initSparse<Scalar>(density, refMat2, m2);
322 refMat3 = refMat2.template triangularView<Lower>();
323 m3 = m2.template triangularView<Lower>();
324 VERIFY_IS_APPROX(m3, refMat3);
325
326 refMat3 = refMat2.template triangularView<Upper>();
327 m3 = m2.template triangularView<Upper>();
328 VERIFY_IS_APPROX(m3, refMat3);
329
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100330 if(inner>=outer) // FIXME this should be implemented for outer>inner as well
331 {
332 refMat3 = refMat2.template triangularView<UnitUpper>();
333 m3 = m2.template triangularView<UnitUpper>();
334 VERIFY_IS_APPROX(m3, refMat3);
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100335
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100336 refMat3 = refMat2.template triangularView<UnitLower>();
337 m3 = m2.template triangularView<UnitLower>();
338 VERIFY_IS_APPROX(m3, refMat3);
339 }
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200340
341 refMat3 = refMat2.template triangularView<StrictlyUpper>();
342 m3 = m2.template triangularView<StrictlyUpper>();
343 VERIFY_IS_APPROX(m3, refMat3);
344
345 refMat3 = refMat2.template triangularView<StrictlyLower>();
346 m3 = m2.template triangularView<StrictlyLower>();
347 VERIFY_IS_APPROX(m3, refMat3);
Gael Guennebaudf6b1dee2015-11-04 17:02:32 +0100348
349 // check sparse-traingular to dense
350 refMat3 = m2.template triangularView<StrictlyUpper>();
351 VERIFY_IS_APPROX(refMat3, DenseMatrix(refMat2.template triangularView<StrictlyUpper>()));
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100352 }
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200353
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100354 // test selfadjointView
Gael Guennebaud9353bba2011-12-04 14:39:24 +0100355 if(!SparseMatrixType::IsRowMajor)
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100356 {
357 DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
358 SparseMatrixType m2(rows, rows), m3(rows, rows);
359 initSparse<Scalar>(density, refMat2, m2);
360 refMat3 = refMat2.template selfadjointView<Lower>();
361 m3 = m2.template selfadjointView<Lower>();
362 VERIFY_IS_APPROX(m3, refMat3);
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100363
364 // selfadjointView only works for square matrices:
365 SparseMatrixType m4(rows, rows+1);
366 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Lower>());
367 VERIFY_RAISES_ASSERT(m4.template selfadjointView<Upper>());
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100368 }
369
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200370 // test sparseView
371 {
372 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
373 SparseMatrixType m2(rows, rows);
Gael Guennebaud9a3ec632010-11-15 14:14:05 +0100374 initSparse<Scalar>(density, refMat2, m2);
Gael Guennebaudfa6d36e2010-07-22 15:57:01 +0200375 VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
376 }
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100377
378 // test diagonal
379 {
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100380 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
381 SparseMatrixType m2(rows, cols);
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100382 initSparse<Scalar>(density, refMat2, m2);
383 VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval());
Gael Guennebaudb26e6972014-12-01 14:41:39 +0100384 VERIFY_IS_APPROX(const_cast<const SparseMatrixType&>(m2).diagonal(), refMat2.diagonal().eval());
385
386 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag);
387 m2.diagonal() += refMat2.diagonal();
388 refMat2.diagonal() += refMat2.diagonal();
389 VERIFY_IS_APPROX(m2, refMat2);
Gael Guennebaud82f9aa12011-12-04 21:49:21 +0100390 }
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200391
Gael Guennebaud62f21e22015-06-24 17:55:00 +0200392 // test diagonal to sparse
393 {
394 DenseVector d = DenseVector::Random(rows);
395 DenseMatrix refMat2 = d.asDiagonal();
396 SparseMatrixType m2(rows, rows);
397 m2 = d.asDiagonal();
398 VERIFY_IS_APPROX(m2, refMat2);
Gael Guennebaud4c8cd132015-06-24 18:11:06 +0200399 SparseMatrixType m3(d.asDiagonal());
400 VERIFY_IS_APPROX(m3, refMat2);
Gael Guennebaud62f21e22015-06-24 17:55:00 +0200401 refMat2 += d.asDiagonal();
402 m2 += d.asDiagonal();
403 VERIFY_IS_APPROX(m2, refMat2);
404 }
405
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200406 // test conservative resize
407 {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100408 std::vector< std::pair<StorageIndex,StorageIndex> > inc;
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100409 if(rows > 3 && cols > 2)
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100410 inc.push_back(std::pair<StorageIndex,StorageIndex>(-3,-2));
411 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,0));
412 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,2));
413 inc.push_back(std::pair<StorageIndex,StorageIndex>(3,0));
414 inc.push_back(std::pair<StorageIndex,StorageIndex>(0,3));
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200415
416 for(size_t i = 0; i< inc.size(); i++) {
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100417 StorageIndex incRows = inc[i].first;
418 StorageIndex incCols = inc[i].second;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200419 SparseMatrixType m1(rows, cols);
420 DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols);
421 initSparse<Scalar>(density, refMat1, m1);
422
423 m1.conservativeResize(rows+incRows, cols+incCols);
424 refMat1.conservativeResize(rows+incRows, cols+incCols);
425 if (incRows > 0) refMat1.bottomRows(incRows).setZero();
426 if (incCols > 0) refMat1.rightCols(incCols).setZero();
427
428 VERIFY_IS_APPROX(m1, refMat1);
429
430 // Insert new values
431 if (incRows > 0)
Gael Guennebaud7ee378d2013-07-12 16:40:02 +0200432 m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200433 if (incCols > 0)
Gael Guennebaud7ee378d2013-07-12 16:40:02 +0200434 m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1;
Benjamin Piwowarski6bf49ce2012-07-19 00:07:06 +0200435
436 VERIFY_IS_APPROX(m1, refMat1);
437
438
439 }
440 }
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200441
442 // test Identity matrix
443 {
444 DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows);
445 SparseMatrixType m1(rows, rows);
446 m1.setIdentity();
447 VERIFY_IS_APPROX(m1, refMat1);
Gael Guennebaud8a211bb2015-10-25 22:01:58 +0100448 for(int k=0; k<rows*rows/4; ++k)
449 {
450 Index i = internal::random<Index>(0,rows-1);
451 Index j = internal::random<Index>(0,rows-1);
Gael Guennebaud73f692d2015-10-27 11:01:37 +0100452 Scalar v = internal::random<Scalar>();
Gael Guennebaud8a211bb2015-10-25 22:01:58 +0100453 m1.coeffRef(i,j) = v;
454 refMat1.coeffRef(i,j) = v;
455 VERIFY_IS_APPROX(m1, refMat1);
456 if(internal::random<Index>(0,10)<2)
457 m1.makeCompressed();
458 }
459 m1.setIdentity();
460 refMat1.setIdentity();
461 VERIFY_IS_APPROX(m1, refMat1);
Desire NUENTSA4cd82452013-06-11 14:42:29 +0200462 }
Gael Guennebaudec469702016-02-01 15:04:33 +0100463
464 // test array/vector of InnerIterator
465 {
466 typedef typename SparseMatrixType::InnerIterator IteratorType;
467
468 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
469 SparseMatrixType m2(rows, cols);
470 initSparse<Scalar>(density, refMat2, m2);
471 IteratorType static_array[2];
472 static_array[0] = IteratorType(m2,0);
473 static_array[1] = IteratorType(m2,m2.outerSize()-1);
474 VERIFY( static_array[0] || m2.innerVector(static_array[0].outer()).nonZeros() == 0 );
475 VERIFY( static_array[1] || m2.innerVector(static_array[1].outer()).nonZeros() == 0 );
476 if(static_array[0] && static_array[1])
477 {
478 ++(static_array[1]);
479 static_array[1] = IteratorType(m2,0);
480 VERIFY( static_array[1] );
481 VERIFY( static_array[1].index() == static_array[0].index() );
482 VERIFY( static_array[1].outer() == static_array[0].outer() );
483 VERIFY( static_array[1].value() == static_array[0].value() );
484 }
485
486 std::vector<IteratorType> iters(2);
487 iters[0] = IteratorType(m2,0);
488 iters[1] = IteratorType(m2,m2.outerSize()-1);
489 }
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000490}
491
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100492
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100493template<typename SparseMatrixType>
Christoph Hertzberge8cdbed2014-12-04 22:48:53 +0100494void big_sparse_triplet(Index rows, Index cols, double density) {
495 typedef typename SparseMatrixType::StorageIndex StorageIndex;
496 typedef typename SparseMatrixType::Scalar Scalar;
497 typedef Triplet<Scalar,Index> TripletType;
498 std::vector<TripletType> triplets;
499 double nelements = density * rows*cols;
500 VERIFY(nelements>=0 && nelements < NumTraits<StorageIndex>::highest());
501 Index ntriplets = Index(nelements);
502 triplets.reserve(ntriplets);
503 Scalar sum = Scalar(0);
504 for(Index i=0;i<ntriplets;++i)
505 {
506 Index r = internal::random<Index>(0,rows-1);
507 Index c = internal::random<Index>(0,cols-1);
508 Scalar v = internal::random<Scalar>();
509 triplets.push_back(TripletType(r,c,v));
510 sum += v;
511 }
512 SparseMatrixType m(rows,cols);
513 m.setFromTriplets(triplets.begin(), triplets.end());
514 VERIFY(m.nonZeros() <= ntriplets);
515 VERIFY_IS_APPROX(sum, m.sum());
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100516}
517
518
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000519void test_sparse_basic()
520{
521 for(int i = 0; i < g_repeat; i++) {
Gael Guennebaudc43154b2015-03-04 10:16:46 +0100522 int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200);
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100523 if(Eigen::internal::random<int>(0,4) == 0) {
524 r = c; // check square matrices in 25% of tries
525 }
526 EIGEN_UNUSED_VARIABLE(r+c);
527 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(1, 1)) ));
Gael Guennebaud91fe1502011-06-07 11:28:16 +0200528 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) ));
Christoph Hertzberg0833b822014-10-31 17:12:13 +0100529 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(r, c)) ));
530 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(r, c)) ));
531 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(r, c)) ));
Gael Guennebaudd7698c12015-03-19 15:11:05 +0100532 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,ColMajor,long int>(r, c)) ));
533 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,RowMajor,long int>(r, c)) ));
Gael Guennebaud47e89022013-06-10 10:34:03 +0200534
Gael Guennebaudc43154b2015-03-04 10:16:46 +0100535 r = Eigen::internal::random<int>(1,100);
536 c = Eigen::internal::random<int>(1,100);
537 if(Eigen::internal::random<int>(0,4) == 0) {
538 r = c; // check square matrices in 25% of tries
539 }
540
Gael Guennebaudd7698c12015-03-19 15:11:05 +0100541 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) ));
542 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) ));
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000543 }
Christoph Hertzbergc5a37772014-10-31 17:19:05 +0100544
545 // Regression test for bug 900: (manually insert higher values here, if you have enough RAM):
546 CALL_SUBTEST_3((big_sparse_triplet<SparseMatrix<float, RowMajor, int> >(10000, 10000, 0.125)));
547 CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125)));
Gael Guennebaudbfd6ee62015-11-06 15:05:37 +0100548
549 // Regression test for bug 1105
550#ifdef EIGEN_TEST_PART_6
551 {
552 int n = Eigen::internal::random<int>(200,600);
553 SparseMatrix<std::complex<double>,0, long> mat(n, n);
554 std::complex<double> val;
555
556 for(int i=0; i<n; ++i)
557 {
558 mat.coeffRef(i, i%(n/10)) = val;
559 VERIFY(mat.data().allocatedSize()<20*n);
560 }
561 }
562#endif
Gael Guennebaud86ccd992008-11-05 13:47:55 +0000563}