blob: 3a8873d430b59bd83b41722121052852a1ca0496 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2010 Gael Guennebaud // // 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.h" template void initSPD(double density, Matrix& refMat, SparseMatrix& sparseMat) { Matrix 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 void sparse_solvers(int rows, int cols) { double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector; // Scalar eps = 1e-6; DenseVector vec1 = DenseVector::Random(rows); std::vector zeroCoords; std::vector nonzeroCoords; // test triangular solver { DenseVector vec2 = vec1, vec3 = vec1; SparseMatrix m2(rows, cols); DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); // lower - dense initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), m2.template triangularView().solve(vec3)); // upper - dense initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), m2.template triangularView().solve(vec3)); VERIFY_IS_APPROX(refMat2.conjugate().template triangularView().solve(vec2), m2.conjugate().template triangularView().solve(vec3)); { SparseMatrix cm2(m2); //Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr MappedSparseMatrix mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr()); VERIFY_IS_APPROX(refMat2.conjugate().template triangularView().solve(vec2), mm2.conjugate().template triangularView().solve(vec3)); } // lower - transpose initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.transpose().template triangularView().solve(vec2), m2.transpose().template triangularView().solve(vec3)); // upper - transpose initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.transpose().template triangularView().solve(vec2), m2.transpose().template triangularView().solve(vec3)); SparseMatrix matB(rows, rows); DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); // lower - sparse initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); initSparse(density, refMatB, matB); refMat2.template triangularView().solveInPlace(refMatB); m2.template triangularView().solveInPlace(matB); VERIFY_IS_APPROX(matB.toDense(), refMatB); // upper - sparse initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); initSparse(density, refMatB, matB); refMat2.template triangularView().solveInPlace(refMatB); m2.template triangularView().solveInPlace(matB); VERIFY_IS_APPROX(matB, refMatB); // test deprecated API initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), m2.template triangularView().solve(vec3)); } } void test_sparse_solvers() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(sparse_solvers(8, 8) ); int s = internal::random(1,300); CALL_SUBTEST_2(sparse_solvers >(s,s) ); CALL_SUBTEST_1(sparse_solvers(s,s) ); } }