diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-03-04 09:34:27 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-03-04 09:34:27 +0100 |
commit | 05274219a7c5fdf04bfda089dc3f9eb2923fcc7e (patch) | |
tree | 1f5faca59ebd457ccf1d4beb9dbbb69a448f9376 /test | |
parent | f8390995127f9f73f2376c43f93eaa27bbad3675 (diff) |
Add a CG-based solver for rectangular least-square problems (bug #975).
Diffstat (limited to 'test')
-rw-r--r-- | test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | test/lscg.cpp | 29 | ||||
-rw-r--r-- | test/sparse_solver.h | 58 |
3 files changed, 84 insertions, 4 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 168749634..1712b8718 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -234,6 +234,7 @@ ei_add_test(sparse_permutations) ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) ei_add_test(bicgstab) +ei_add_test(lscg) ei_add_test(sparselu) ei_add_test(sparseqr) ei_add_test(umeyama) diff --git a/test/lscg.cpp b/test/lscg.cpp new file mode 100644 index 000000000..599ed5619 --- /dev/null +++ b/test/lscg.cpp @@ -0,0 +1,29 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud <g.gael@free.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/. + +#include "sparse_solver.h" +#include <Eigen/IterativeLinearSolvers> + +template<typename T> void test_lscg_T() +{ + LSCG<SparseMatrix<T> > lscg_colmajor_diag; + LSCG<SparseMatrix<T>, IdentityPreconditioner> lscg_colmajor_I; + + CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_diag) ); + CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_I) ); + + CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_diag) ); + CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_I) ); +} + +void test_lscg() +{ + CALL_SUBTEST_1(test_lscg_T<double>()); + CALL_SUBTEST_2(test_lscg_T<std::complex<double> >()); +} diff --git a/test/sparse_solver.h b/test/sparse_solver.h index f0a4691e3..f266e2c9a 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -17,9 +17,9 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, typedef typename Mat::Scalar Scalar; typedef typename Mat::StorageIndex StorageIndex; - DenseRhs refX = dA.lu().solve(db); + DenseRhs refX = dA.householderQr().solve(db); { - Rhs x(b.rows(), b.cols()); + Rhs x(A.cols(), b.cols()); Rhs oldb = b; solver.compute(A); @@ -94,7 +94,7 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, // test dense Block as the result and rhs: { - DenseRhs x(db.rows(), db.cols()); + DenseRhs x(refX.rows(), refX.cols()); DenseRhs oldb(db); x.setZero(); x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols())); @@ -119,7 +119,7 @@ void check_sparse_solving_real_cases(Solver& solver, const typename Solver::Matr typedef typename Mat::Scalar Scalar; typedef typename Mat::RealScalar RealScalar; - Rhs x(b.rows(), b.cols()); + Rhs x(A.cols(), b.cols()); solver.compute(A); if (solver.info() != Success) @@ -410,3 +410,53 @@ template<typename Solver> void check_sparse_square_abs_determinant(Solver& solve } } +template<typename Solver, typename DenseMat> +void generate_sparse_leastsquare_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300, int options = ForceNonZeroDiag) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + + int rows = internal::random<int>(1,maxSize); + int cols = internal::random<int>(1,rows); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,cols); + dA.resize(rows,cols); + + initSparse<Scalar>(density, dA, A, options); +} + +template<typename Solver> void check_sparse_leastsquare_solving(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef SparseMatrix<Scalar,ColMajor> SpMat; + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + + int rhsCols = internal::random<int>(1,16); + + Mat A; + DenseMatrix dA; + for (int i = 0; i < g_repeat; i++) { + generate_sparse_leastsquare_problem(solver, A, dA); + + A.makeCompressed(); + DenseVector b = DenseVector::Random(A.rows()); + DenseMatrix dB(A.rows(),rhsCols); + SpMat B(A.rows(),rhsCols); + double density = (std::max)(8./(A.rows()*rhsCols), 0.1); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + B.makeCompressed(); + check_sparse_solving(solver, A, b, dA, b); + check_sparse_solving(solver, A, dB, dA, dB); + check_sparse_solving(solver, A, B, dA, dB); + + // check only once + if(i==0) + { + b = DenseVector::Zero(A.rows()); + check_sparse_solving(solver, A, b, dA, b); + } + } +} |