260 lines
12 KiB
C++
260 lines
12 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
// SPDX-License-Identifier: MPL-2.0
|
|
|
|
#include "sparse.h"
|
|
|
|
template <typename Scalar>
|
|
void initSPD(double density, Matrix<Scalar, Dynamic, Dynamic>& refMat, SparseMatrix<Scalar>& sparseMat) {
|
|
Matrix<Scalar, Dynamic, Dynamic> aux(refMat.rows(), refMat.cols());
|
|
initSparse(density, refMat, sparseMat);
|
|
refMat = refMat * refMat.adjoint();
|
|
for (int k = 0; k < 2; ++k) {
|
|
initSparse(density, aux, sparseMat, ForceNonZeroDiag);
|
|
refMat += aux * aux.adjoint();
|
|
}
|
|
sparseMat.setZero();
|
|
for (int j = 0; j < sparseMat.cols(); ++j)
|
|
for (int i = j; i < sparseMat.rows(); ++i)
|
|
if (refMat(i, j) != Scalar(0)) sparseMat.insert(i, j) = refMat(i, j);
|
|
sparseMat.finalize();
|
|
}
|
|
|
|
template <typename Scalar>
|
|
void sparse_solvers(int rows, int cols) {
|
|
double density = (std::max)(8. / (rows * cols), 0.01);
|
|
typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
|
|
typedef Matrix<Scalar, Dynamic, 1> DenseVector;
|
|
// Scalar eps = 1e-6;
|
|
|
|
DenseVector vec1 = DenseVector::Random(rows);
|
|
|
|
std::vector<Vector2i> zeroCoords;
|
|
std::vector<Vector2i> nonzeroCoords;
|
|
|
|
// test triangular solver
|
|
{
|
|
DenseVector vec2 = vec1, vec3 = vec1;
|
|
SparseMatrix<Scalar> m2(rows, cols);
|
|
DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
|
|
|
|
// lower - dense
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
|
|
VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
|
|
m2.template triangularView<Lower>().solve(vec3));
|
|
|
|
// upper - dense
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
|
|
VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
|
|
m2.template triangularView<Upper>().solve(vec3));
|
|
VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
|
|
m2.conjugate().template triangularView<Upper>().solve(vec3));
|
|
{
|
|
SparseMatrix<Scalar> cm2(m2);
|
|
// Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
|
|
Map<SparseMatrix<Scalar> > mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(),
|
|
cm2.valuePtr());
|
|
VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
|
|
mm2.conjugate().template triangularView<Upper>().solve(vec3));
|
|
}
|
|
|
|
// lower - transpose
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
|
|
VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
|
|
m2.transpose().template triangularView<Upper>().solve(vec3));
|
|
|
|
// upper - transpose
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
|
|
VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
|
|
m2.transpose().template triangularView<Lower>().solve(vec3));
|
|
|
|
SparseMatrix<Scalar> matB(rows, rows);
|
|
DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);
|
|
|
|
// lower - sparse
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular);
|
|
initSparse<Scalar>(density, refMatB, matB);
|
|
refMat2.template triangularView<Lower>().solveInPlace(refMatB);
|
|
m2.template triangularView<Lower>().solveInPlace(matB);
|
|
VERIFY_IS_APPROX(matB.toDense(), refMatB);
|
|
|
|
// upper - sparse
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular);
|
|
initSparse<Scalar>(density, refMatB, matB);
|
|
refMat2.template triangularView<Upper>().solveInPlace(refMatB);
|
|
m2.template triangularView<Upper>().solveInPlace(matB);
|
|
VERIFY_IS_APPROX(matB, refMatB);
|
|
|
|
// A triangularView is a view of the triangular PART of a possibly-general matrix,
|
|
// so the stored matrix need not be strictly triangular. Exercise a general lhs, a
|
|
// SparseVector rhs, a mismatched-StorageIndex rhs, an uncompressed lhs, an
|
|
// expression lhs, and a unit diagonal -- none of which the checks above cover.
|
|
{
|
|
SparseMatrix<Scalar> mg(rows, rows);
|
|
DenseMatrix refMatG = DenseMatrix::Zero(rows, rows);
|
|
initSparse<Scalar>(density, refMatG, mg, ForceNonZeroDiag); // GENERAL (both triangles stored)
|
|
initSparse<Scalar>(density, refMatB, matB);
|
|
|
|
// general matrix through a lower / upper / unit-upper view, sparse rhs
|
|
for (int mode = 0; mode < 3; ++mode) {
|
|
DenseMatrix rb = refMatB;
|
|
SparseMatrix<Scalar> mb = matB;
|
|
if (mode == 0) {
|
|
refMatG.template triangularView<Lower>().solveInPlace(rb);
|
|
mg.template triangularView<Lower>().solveInPlace(mb);
|
|
} else if (mode == 1) {
|
|
refMatG.template triangularView<Upper>().solveInPlace(rb);
|
|
mg.template triangularView<Upper>().solveInPlace(mb);
|
|
} else {
|
|
refMatG.template triangularView<UnitUpper>().solveInPlace(rb);
|
|
mg.template triangularView<UnitUpper>().solveInPlace(mb);
|
|
}
|
|
VERIFY_IS_APPROX(mb.toDense(), rb);
|
|
}
|
|
|
|
// expression lhs (no raw storage -> iterator path)
|
|
{
|
|
DenseMatrix rb = refMatB;
|
|
SparseMatrix<Scalar> mb = matB;
|
|
refMatG.template triangularView<Upper>().solveInPlace(rb);
|
|
(Scalar(1) * mg).template triangularView<Upper>().solveInPlace(mb);
|
|
VERIFY_IS_APPROX(mb.toDense(), rb);
|
|
}
|
|
|
|
// uncompressed lhs (innerNonZeroPtr != null)
|
|
{
|
|
DenseMatrix rb = refMatB;
|
|
SparseMatrix<Scalar> mb = matB, mu = mg;
|
|
mu.reserve(Matrix<int, Dynamic, 1>::Constant(mu.cols(), rows)); // -> uncompressed
|
|
refMatG.template triangularView<Lower>().solveInPlace(rb);
|
|
mu.template triangularView<Lower>().solveInPlace(mb);
|
|
VERIFY_IS_APPROX(mb.toDense(), rb);
|
|
}
|
|
|
|
// mismatched-StorageIndex rhs (-> InnerIterator fallback)
|
|
{
|
|
DenseMatrix rb = refMatB;
|
|
SparseMatrix<Scalar, ColMajor, long> mbl = matB;
|
|
refMatG.template triangularView<Lower>().solveInPlace(rb);
|
|
mg.template triangularView<Lower>().solveInPlace(mbl);
|
|
VERIFY_IS_APPROX(DenseMatrix(mbl), rb);
|
|
}
|
|
|
|
// SparseVector rhs (sets CompressedAccessBit but its outerIndexPtr() is null)
|
|
{
|
|
DenseVector rv = DenseVector::Zero(rows);
|
|
SparseVector<Scalar> vb(rows);
|
|
for (Index i = 0; i < rows; ++i)
|
|
if (internal::random<int>(0, 2) == 0) {
|
|
Scalar s = internal::random<Scalar>();
|
|
vb.coeffRef(i) = s;
|
|
rv(i) = s;
|
|
}
|
|
DenseVector rref = refMatG.template triangularView<Lower>().solve(rv);
|
|
SparseVector<Scalar> vx = vb;
|
|
mg.template triangularView<Lower>().solveInPlace(vx);
|
|
VERIFY_IS_APPROX(DenseVector(vx), rref);
|
|
}
|
|
|
|
// explicitly-stored zero rhs entries must not expand into stored zeros: the reach
|
|
// is a structural bound, so a zero rhs coefficient (or one that cancels to zero)
|
|
// is pruned at insertion rather than materialized across the whole reach.
|
|
{
|
|
SparseMatrix<Scalar> mb(rows, matB.cols());
|
|
DenseMatrix rb = DenseMatrix::Zero(rows, matB.cols());
|
|
for (Index c = 0; c < mb.cols(); ++c)
|
|
for (Index i = 0; i < rows; ++i)
|
|
if (internal::random<int>(0, 3) == 0) {
|
|
Scalar s = internal::random<int>(0, 2) == 0 ? Scalar(0) : internal::random<Scalar>(); // some explicit 0
|
|
mb.insert(i, c) = s;
|
|
rb(i, c) = s;
|
|
}
|
|
mb.makeCompressed();
|
|
refMatG.template triangularView<Lower>().solveInPlace(rb);
|
|
mg.template triangularView<Lower>().solveInPlace(mb);
|
|
VERIFY_IS_APPROX(mb.toDense(), rb);
|
|
for (Index c = 0; c < mb.cols(); ++c)
|
|
for (typename SparseMatrix<Scalar>::InnerIterator it(mb, c); it; ++it)
|
|
VERIFY(!numext::is_exactly_zero(it.value())); // no stored zeros
|
|
}
|
|
|
|
// A reached column with no stored diagonal (non-unit) is out of contract: it must
|
|
// assert in debug on every path (pointer/iterator x lower/upper), rather than the
|
|
// failure being silently keyed to has_compressed_access. In release these divide by
|
|
// zero -> inf/NaN with no out-of-bounds read (covered by the sanitizer drivers).
|
|
{
|
|
SparseMatrix<Scalar> us(3, 3);
|
|
us.insert(0, 1) = Scalar(1);
|
|
us.insert(1, 1) = Scalar(2);
|
|
us.insert(2, 2) = Scalar(3);
|
|
us.makeCompressed(); // column 0 empty -> reached from rhs(0) but no diagonal
|
|
SparseMatrix<Scalar> ub(3, 1);
|
|
ub.insert(0, 0) = Scalar(1);
|
|
ub.makeCompressed();
|
|
SparseMatrix<Scalar> up = ub, ui = ub;
|
|
VERIFY_RAISES_ASSERT(us.template triangularView<Upper>().solveInPlace(up)); // pointer upper
|
|
VERIFY_RAISES_ASSERT((Scalar(1) * us).template triangularView<Upper>().solveInPlace(ui)); // iterator upper
|
|
|
|
SparseMatrix<Scalar> ls(3, 3);
|
|
ls.insert(0, 0) = Scalar(2);
|
|
ls.insert(1, 1) = Scalar(3);
|
|
ls.makeCompressed(); // column 2 empty (last) -> reached from rhs(2) but no diagonal
|
|
SparseMatrix<Scalar> lb(3, 1);
|
|
lb.insert(2, 0) = Scalar(1);
|
|
lb.makeCompressed();
|
|
SparseMatrix<Scalar> lp = lb, li = lb;
|
|
VERIFY_RAISES_ASSERT(ls.template triangularView<Lower>().solveInPlace(lp)); // pointer lower
|
|
VERIFY_RAISES_ASSERT((Scalar(1) * ls).template triangularView<Lower>().solveInPlace(li)); // iterator lower
|
|
}
|
|
|
|
// mixed-scalar rhs: a real lhs applied to a rhs must accumulate in the rhs scalar.
|
|
// For real Scalar this is the ordinary path; for complex Scalar it is the
|
|
// real-factor / complex-data case that must both compile and be correct.
|
|
{
|
|
typedef typename NumTraits<Scalar>::Real Real;
|
|
SparseMatrix<Real> mr(rows, rows);
|
|
Matrix<Real, Dynamic, Dynamic> refMatR = Matrix<Real, Dynamic, Dynamic>::Zero(rows, rows);
|
|
initSparse<Real>(density, refMatR, mr, ForceNonZeroDiag);
|
|
DenseMatrix rb = refMatB;
|
|
SparseMatrix<Scalar> mb = matB;
|
|
refMatR.template cast<Scalar>().template triangularView<Lower>().solveInPlace(rb);
|
|
mr.template triangularView<Lower>().solveInPlace(mb); // SparseMatrix<Real> lhs, <Scalar> rhs
|
|
VERIFY_IS_APPROX(mb.toDense(), rb);
|
|
}
|
|
}
|
|
|
|
// test deprecated API
|
|
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
|
|
VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
|
|
m2.template triangularView<Lower>().solve(vec3));
|
|
|
|
// test empty triangular matrix
|
|
{
|
|
m2.resize(0, 0);
|
|
refMatB.resize(0, refMatB.cols());
|
|
DenseMatrix res = m2.template triangularView<Lower>().solve(refMatB);
|
|
VERIFY_IS_EQUAL(res.rows(), 0);
|
|
VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
|
|
res = refMatB;
|
|
m2.template triangularView<Lower>().solveInPlace(res);
|
|
VERIFY_IS_EQUAL(res.rows(), 0);
|
|
VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
|
|
}
|
|
}
|
|
}
|
|
|
|
EIGEN_DECLARE_TEST(sparse_solvers) {
|
|
for (int i = 0; i < g_repeat; i++) {
|
|
CALL_SUBTEST_1(sparse_solvers<double>(8, 8));
|
|
int s = internal::random<int>(1, 300);
|
|
CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s, s));
|
|
CALL_SUBTEST_1(sparse_solvers<double>(s, s));
|
|
}
|
|
}
|