sparse_solvers.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-2010 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 
12 template<typename Scalar> void
15  SparseMatrix<Scalar>& sparseMat)
16 {
17  Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
18  initSparse(density,refMat,sparseMat);
19  refMat = refMat * refMat.adjoint();
20  for (int k=0; k<2; ++k)
21  {
22  initSparse(density,aux,sparseMat,ForceNonZeroDiag);
23  refMat += aux * aux.adjoint();
24  }
25  sparseMat.setZero();
26  for (int j=0 ; j<sparseMat.cols(); ++j)
27  for (int i=j ; i<sparseMat.rows(); ++i)
28  if (refMat(i,j)!=Scalar(0))
29  sparseMat.insert(i,j) = refMat(i,j);
30  sparseMat.finalize();
31 }
32 
33 template<typename Scalar> void sparse_solvers(int rows, int cols)
34 {
35  double density = (std::max)(8./(rows*cols), 0.01);
38  // Scalar eps = 1e-6;
39 
40  DenseVector vec1 = DenseVector::Random(rows);
41 
42  std::vector<Vector2i> zeroCoords;
43  std::vector<Vector2i> nonzeroCoords;
44 
45  // test triangular solver
46  {
47  DenseVector vec2 = vec1, vec3 = vec1;
48  SparseMatrix<Scalar> m2(rows, cols);
49  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
50 
51  // lower - dense
52  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
53  VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
54  m2.template triangularView<Lower>().solve(vec3));
55 
56  // upper - dense
57  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
58  VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
59  m2.template triangularView<Upper>().solve(vec3));
60  VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
61  m2.conjugate().template triangularView<Upper>().solve(vec3));
62  {
63  SparseMatrix<Scalar> cm2(m2);
64  //Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
65  MappedSparseMatrix<Scalar> mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr());
66  VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
67  mm2.conjugate().template triangularView<Upper>().solve(vec3));
68  }
69 
70  // lower - transpose
71  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
72  VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
73  m2.transpose().template triangularView<Upper>().solve(vec3));
74 
75  // upper - transpose
76  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
77  VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
78  m2.transpose().template triangularView<Lower>().solve(vec3));
79 
80  SparseMatrix<Scalar> matB(rows, rows);
81  DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);
82 
83  // lower - sparse
84  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular);
85  initSparse<Scalar>(density, refMatB, matB);
86  refMat2.template triangularView<Lower>().solveInPlace(refMatB);
87  m2.template triangularView<Lower>().solveInPlace(matB);
88  VERIFY_IS_APPROX(matB.toDense(), refMatB);
89 
90  // upper - sparse
91  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular);
92  initSparse<Scalar>(density, refMatB, matB);
93  refMat2.template triangularView<Upper>().solveInPlace(refMatB);
94  m2.template triangularView<Upper>().solveInPlace(matB);
95  VERIFY_IS_APPROX(matB, refMatB);
96 
97  // test deprecated API
98  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
99  VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
100  m2.template triangularView<Lower>().solve(vec3));
101 
102  // test empty triangular matrix
103  {
104  m2.resize(0,0);
105  refMatB.resize(0,refMatB.cols());
106  DenseMatrix res = m2.template triangularView<Lower>().solve(refMatB);
107  VERIFY_IS_EQUAL(res.rows(),0);
108  VERIFY_IS_EQUAL(res.cols(),refMatB.cols());
109  res = refMatB;
110  m2.template triangularView<Lower>().solveInPlace(res);
111  VERIFY_IS_EQUAL(res.rows(),0);
112  VERIFY_IS_EQUAL(res.cols(),refMatB.cols());
113  }
114  }
115 }
116 
118 {
119  for(int i = 0; i < g_repeat; i++) {
120  CALL_SUBTEST_1(sparse_solvers<double>(8, 8) );
121  int s = internal::random<int>(1,300);
122  CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s,s) );
123  CALL_SUBTEST_1(sparse_solvers<double>(s,s) );
124  }
125 }
EIGEN_DECLARE_TEST(sparse_solvers)
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Index cols() const
Definition: SparseMatrix.h:140
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:168
Index rows() const
Definition: SparseMatrix.h:138
MatrixType m2(n_dims)
MatrixXf matB(2, 2)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
void sparse_solvers(int rows, int cols)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Matrix< Scalar, Dynamic, 1 > DenseVector
#define VERIFY_IS_APPROX(a, b)
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
#define CALL_SUBTEST_1(FUNC)
static int g_repeat
Definition: main.h:169
Scalar & insert(Index row, Index col)
RealScalar s
void initSparse(double density, Matrix< Scalar, Dynamic, Dynamic, Opt1 > &refMat, SparseMatrix< Scalar, Opt2, StorageIndex > &sparseMat, int flags=0, std::vector< Matrix< StorageIndex, 2, 1 > > *zeroCoords=0, std::vector< Matrix< StorageIndex, 2, 1 > > *nonzeroCoords=0)
Definition: sparse.h:51
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:159
RowVectorXd vec1(3)
A triangularView< Lower >().adjoint().solveInPlace(B)
static Vector9 vec3(const Matrix3 &R)
Definition: SO3.cpp:342
#define CALL_SUBTEST_2(FUNC)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
void initSPD(double density, Matrix< Scalar, Dynamic, Dynamic > &refMat, SparseMatrix< Scalar > &sparseMat)
const Scalar * valuePtr() const
Definition: SparseMatrix.h:150
std::ptrdiff_t j


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:54