1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2015-2016 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 // #define EIGEN_DONT_VECTORIZE
10 // #define EIGEN_MAX_ALIGN_BYTES 0
11 #include "sparse_solver.h"
12 #include <Eigen/IterativeLinearSolvers>
13 #include <unsupported/Eigen/IterativeSolvers>
14
test_incomplete_cholesky_T()15 template<typename T, typename I> void test_incomplete_cholesky_T()
16 {
17 typedef SparseMatrix<T,0,I> SparseMatrixType;
18 ConjugateGradient<SparseMatrixType, Lower, IncompleteCholesky<T, Lower, AMDOrdering<I> > > cg_illt_lower_amd;
19 ConjugateGradient<SparseMatrixType, Lower, IncompleteCholesky<T, Lower, NaturalOrdering<I> > > cg_illt_lower_nat;
20 ConjugateGradient<SparseMatrixType, Upper, IncompleteCholesky<T, Upper, AMDOrdering<I> > > cg_illt_upper_amd;
21 ConjugateGradient<SparseMatrixType, Upper, IncompleteCholesky<T, Upper, NaturalOrdering<I> > > cg_illt_upper_nat;
22 ConjugateGradient<SparseMatrixType, Upper|Lower, IncompleteCholesky<T, Lower, AMDOrdering<I> > > cg_illt_uplo_amd;
23
24
25 CALL_SUBTEST( check_sparse_spd_solving(cg_illt_lower_amd) );
26 CALL_SUBTEST( check_sparse_spd_solving(cg_illt_lower_nat) );
27 CALL_SUBTEST( check_sparse_spd_solving(cg_illt_upper_amd) );
28 CALL_SUBTEST( check_sparse_spd_solving(cg_illt_upper_nat) );
29 CALL_SUBTEST( check_sparse_spd_solving(cg_illt_uplo_amd) );
30 }
31
test_incomplete_cholesky()32 void test_incomplete_cholesky()
33 {
34 CALL_SUBTEST_1(( test_incomplete_cholesky_T<double,int>() ));
35 CALL_SUBTEST_2(( test_incomplete_cholesky_T<std::complex<double>, int>() ));
36 CALL_SUBTEST_3(( test_incomplete_cholesky_T<double,long int>() ));
37
38 #ifdef EIGEN_TEST_PART_1
39 // regression for bug 1150
40 for(int N = 1; N<20; ++N)
41 {
42 Eigen::MatrixXd b( N, N );
43 b.setOnes();
44
45 Eigen::SparseMatrix<double> m( N, N );
46 m.reserve(Eigen::VectorXi::Constant(N,4));
47 for( int i = 0; i < N; ++i )
48 {
49 m.insert( i, i ) = 1;
50 m.coeffRef( i, i / 2 ) = 2;
51 m.coeffRef( i, i / 3 ) = 2;
52 m.coeffRef( i, i / 4 ) = 2;
53 }
54
55 Eigen::SparseMatrix<double> A;
56 A = m * m.transpose();
57
58 Eigen::ConjugateGradient<Eigen::SparseMatrix<double>,
59 Eigen::Lower | Eigen::Upper,
60 Eigen::IncompleteCholesky<double> > solver( A );
61 VERIFY(solver.preconditioner().info() == Eigen::Success);
62 VERIFY(solver.info() == Eigen::Success);
63 }
64 #endif
65 }
66