1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "sparse.h" 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> void 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathinitSPD(double density, 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,Dynamic,Dynamic>& refMat, 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar>& sparseMat) 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols()); 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse(density,refMat,sparseMat); 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath refMat = refMat * refMat.adjoint(); 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k<2; ++k) 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse(density,aux,sparseMat,ForceNonZeroDiag); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath refMat += aux * aux.adjoint(); 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sparseMat.setZero(); 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j=0 ; j<sparseMat.cols(); ++j) 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=j ; i<sparseMat.rows(); ++i) 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (refMat(i,j)!=Scalar(0)) 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sparseMat.insert(i,j) = refMat(i,j); 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sparseMat.finalize(); 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> void sparse_solvers(int rows, int cols) 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double density = (std::max)(8./(rows*cols), 0.01); 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> DenseVector; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Scalar eps = 1e-6; 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseVector vec1 = DenseVector::Random(rows); 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::vector<Vector2i> zeroCoords; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::vector<Vector2i> nonzeroCoords; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test triangular solver 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseVector vec2 = vec1, vec3 = vec1; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar> m2(rows, cols); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lower - dense 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.template triangularView<Lower>().solve(vec3)); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // upper - dense 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2), 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.template triangularView<Upper>().solve(vec3)); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2), 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.conjugate().template triangularView<Upper>().solve(vec3)); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar> cm2(m2); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MappedSparseMatrix<Scalar> mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr()); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2), 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mm2.conjugate().template triangularView<Upper>().solve(vec3)); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lower - transpose 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2), 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.transpose().template triangularView<Upper>().solve(vec3)); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // upper - transpose 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2), 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.transpose().template triangularView<Lower>().solve(vec3)); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar> matB(rows, rows); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lower - sparse 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMatB, matB); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath refMat2.template triangularView<Lower>().solveInPlace(refMatB); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.template triangularView<Lower>().solveInPlace(matB); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(matB.toDense(), refMatB); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // upper - sparse 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMatB, matB); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath refMat2.template triangularView<Upper>().solveInPlace(refMatB); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.template triangularView<Upper>().solveInPlace(matB); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(matB, refMatB); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test deprecated API 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m2.template triangularView<Lower>().solve(vec3)); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_sparse_solvers() 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(sparse_solvers<double>(8, 8) ); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int s = internal::random<int>(1,300); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s,s) ); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(sparse_solvers<double>(s,s) ); 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 113