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

gtsam
Author(s):
autogenerated on Wed May 28 2025 03:03:23