Functions
sparse_solvers.cpp File Reference
#include "sparse.h"

Functions

template<typename Scalar >
void initSPD (double density, Matrix< Scalar, Dynamic, Dynamic > &refMat, SparseMatrix< Scalar > &sparseMat)
 
template<typename Scalar >
void sparse_solvers (int rows, int cols)
 
 EIGEN_DECLARE_TEST (sparse_solvers)
 

Function Documentation

◆ EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST ( sparse_solvers  )
113  {
114  for (int i = 0; i < g_repeat; i++) {
115  CALL_SUBTEST_1(sparse_solvers<double>(8, 8));
116  int s = internal::random<int>(1, 300);
117  CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s, s));
118  CALL_SUBTEST_1(sparse_solvers<double>(s, s));
119  }
120 }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
s
RealScalar s
Definition: level1_cplx_impl.h:130
Eigen::g_repeat
static int g_repeat
Definition: main.h:191
sparse_solvers
void sparse_solvers(int rows, int cols)
Definition: sparse_solvers.cpp:29
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10

References CALL_SUBTEST_1, CALL_SUBTEST_2, Eigen::g_repeat, i, s, and sparse_solvers().

◆ initSPD()

template<typename Scalar >
void initSPD ( double  density,
Matrix< Scalar, Dynamic, Dynamic > &  refMat,
SparseMatrix< Scalar > &  sparseMat 
)
13  {
14  Matrix<Scalar, Dynamic, Dynamic> aux(refMat.rows(), refMat.cols());
15  initSparse(density, refMat, sparseMat);
16  refMat = refMat * refMat.adjoint();
17  for (int k = 0; k < 2; ++k) {
18  initSparse(density, aux, sparseMat, ForceNonZeroDiag);
19  refMat += aux * aux.adjoint();
20  }
21  sparseMat.setZero();
22  for (int j = 0; j < sparseMat.cols(); ++j)
23  for (int i = j; i < sparseMat.rows(); ++i)
24  if (refMat(i, j) != Scalar(0)) sparseMat.insert(i, j) = refMat(i, j);
25  sparseMat.finalize();
26 }
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
Eigen::PlainObjectBase::cols
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: PlainObjectBase.h:192
Eigen::PlainObjectBase::rows
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: PlainObjectBase.h:191
Eigen::SparseMatrix::setZero
void setZero()
Definition: SparseMatrix.h:303
Eigen::SparseMatrix::cols
Index cols() const
Definition: SparseMatrix.h:161
Eigen::SparseMatrix::finalize
void finalize()
Definition: SparseMatrix.h:461
Eigen::SparseMatrix::rows
Index rows() const
Definition: SparseMatrix.h:159
Eigen::SparseMatrix::insert
Scalar & insert(Index row, Index col)
Definition: SparseMatrix.h:1586
k
char char char int int * k
Definition: level2_impl.h:374
UniformPSDSelfTest.density
density
Definition: UniformPSDSelfTest.py:19
ForceNonZeroDiag
@ ForceNonZeroDiag
Definition: sparse.h:32
initSparse
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:42
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::cols(), Eigen::PlainObjectBase< Derived >::cols(), UniformPSDSelfTest::density, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::finalize(), ForceNonZeroDiag, i, initSparse(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::insert(), j, k, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::rows(), Eigen::PlainObjectBase< Derived >::rows(), and Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::setZero().

◆ sparse_solvers()

template<typename Scalar >
void sparse_solvers ( int  rows,
int  cols 
)
29  {
30  double density = (std::max)(8. / (rows * cols), 0.01);
31  typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
32  typedef Matrix<Scalar, Dynamic, 1> DenseVector;
33  // Scalar eps = 1e-6;
34 
35  DenseVector vec1 = DenseVector::Random(rows);
36 
37  std::vector<Vector2i> zeroCoords;
38  std::vector<Vector2i> nonzeroCoords;
39 
40  // test triangular solver
41  {
42  DenseVector vec2 = vec1, vec3 = vec1;
43  SparseMatrix<Scalar> m2(rows, cols);
44  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
45 
46  // lower - dense
47  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
48  VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
49  m2.template triangularView<Lower>().solve(vec3));
50 
51  // upper - dense
52  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
53  VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
54  m2.template triangularView<Upper>().solve(vec3));
55  VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
56  m2.conjugate().template triangularView<Upper>().solve(vec3));
57  {
58  SparseMatrix<Scalar> cm2(m2);
59  // Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
60  Map<SparseMatrix<Scalar> > mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(),
61  cm2.valuePtr());
62  VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
63  mm2.conjugate().template triangularView<Upper>().solve(vec3));
64  }
65 
66  // lower - transpose
67  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
68  VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
69  m2.transpose().template triangularView<Upper>().solve(vec3));
70 
71  // upper - transpose
72  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
73  VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
74  m2.transpose().template triangularView<Lower>().solve(vec3));
75 
76  SparseMatrix<Scalar> matB(rows, rows);
77  DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);
78 
79  // lower - sparse
80  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular);
81  initSparse<Scalar>(density, refMatB, matB);
82  refMat2.template triangularView<Lower>().solveInPlace(refMatB);
83  m2.template triangularView<Lower>().solveInPlace(matB);
84  VERIFY_IS_APPROX(matB.toDense(), refMatB);
85 
86  // upper - sparse
87  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular);
88  initSparse<Scalar>(density, refMatB, matB);
89  refMat2.template triangularView<Upper>().solveInPlace(refMatB);
90  m2.template triangularView<Upper>().solveInPlace(matB);
91  VERIFY_IS_APPROX(matB, refMatB);
92 
93  // test deprecated API
94  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
95  VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
96  m2.template triangularView<Lower>().solve(vec3));
97 
98  // test empty triangular matrix
99  {
100  m2.resize(0, 0);
101  refMatB.resize(0, refMatB.cols());
102  DenseMatrix res = m2.template triangularView<Lower>().solve(refMatB);
103  VERIFY_IS_EQUAL(res.rows(), 0);
104  VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
105  res = refMatB;
106  m2.template triangularView<Lower>().solveInPlace(res);
107  VERIFY_IS_EQUAL(res.rows(), 0);
108  VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
109  }
110  }
111 }
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
DenseVector
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: BenchSparseUtil.h:24
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
matB
MatrixXf matB(2, 2)
vec1
RowVectorXd vec1(3)
m2
MatrixType m2(n_dims)
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:294
max
#define max(a, b)
Definition: datatypes.h:23
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:13
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:367
oomph::PseudoSolidHelper::Zero
double Zero
Definition: pseudosolid_node_update_elements.cc:35
MakeUpperTriangular
@ MakeUpperTriangular
Definition: sparse.h:32
MakeLowerTriangular
@ MakeLowerTriangular
Definition: sparse.h:32

References cols, Eigen::PlainObjectBase< Derived >::cols(), UniformPSDSelfTest::density, ForceNonZeroDiag, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::innerIndexPtr(), m2(), MakeLowerTriangular, MakeUpperTriangular, matB(), max, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::nonZeros(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::outerIndexPtr(), res, Eigen::PlainObjectBase< Derived >::resize(), rows, Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::valuePtr(), vec1(), VERIFY_IS_APPROX, VERIFY_IS_EQUAL, and oomph::PseudoSolidHelper::Zero.

Referenced by EIGEN_DECLARE_TEST().