/******************************************************************************* * CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps * * version 0.1 * * Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg * * * * This library is free software; you can redistribute it and/or modify it * * under the terms of the GNU Lesser General Public License as published by the * * Free Software Foundation; either version 2.1 of the License, or (at your * * option) any later version. * * * * This library is distributed in the hope that it will be useful, but WITHOUT * * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * * for more details. * * * * You should have received a copy of the GNU Lesser General Public License * * along with this library; if not, write to the Free Software Foundation, * * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. * * * * Web site: http://cgogn.unistra.fr/ * * Contact information: cgogn@unistra.fr * * * *******************************************************************************/ /* * This interface to CHOLMOD is largely inspired from the one found in * Graphite developped by Bruno Levy, ALICE Project, LORIA, INRIA Lorraine */ #ifndef __CHOLESKY_SOLVER__ #define __CHOLESKY_SOLVER__ #include #include template class Solver_CHOLESKY { public: Solver_CHOLESKY() { N = 0 ; factorized_ = false ; L = NULL ; // Initialize CHOLMOD library cholmod_start(&c) ; } ~Solver_CHOLESKY() { cholmod_free_factor(&L, &c) ; cholmod_finish(&c) ; } bool factorized() { return factorized_ ; } void factorize(const Eigen::SparseMatrix& A_in) { assert(A_in.rows() == A_in.cols()) ; // assert(A_in.has_symmetric_storage()) ; N = A_in.rows() ; // Translate sparse matrix into cholmod format int NNZ = A_in.nonZeros() ; cholmod_sparse* A = cholmod_allocate_sparse(N, N, NNZ, false, true, -1, CHOLMOD_REAL, &c); int* colptr = static_cast(A->p) ; int* rowind = static_cast(A->i) ; double* a = static_cast(A->x) ; // Convert local Matrix into CHOLMOD Matrix int count = 0 ; for(int j = 0; j < N; j++) { colptr[j] = count; for(typename Eigen::SparseMatrix::InnerIterator it(A_in, j); it; ++it) { a[count] = it.value(); rowind[count] = it.row(); ++count; } } colptr[N] = NNZ ; // Factorize L = cholmod_analyze(A, &c) ; cholmod_factorize(A, L, &c) ; factorized_ = true ; // Cleanup cholmod_free_sparse(&A, &c) ; } void solve(const Eigen::Matrix& b_in, Eigen::Matrix& x_out) { assert(factorized_) ; assert(L->n == b_in.rows()) ; assert(L->n == x_out.rows()) ; // Translate right-hand side into cholmod format cholmod_dense* b = cholmod_allocate_dense(N, 1, N, CHOLMOD_REAL, &c) ; double* cb = static_cast(b->x) ; for(int i = 0; i < N; ++i) cb[i] = b_in[i] ; // Solve cholmod_dense* x = cholmod_solve(CHOLMOD_A, L, b, &c) ; // Get result double* cx = static_cast(x->x) ; for(int i = 0; i < N; ++i) x_out[i] = cx[i] ; // Cleanup cholmod_free_dense(&b, &c) ; cholmod_free_dense(&x, &c) ; } void reset() { if(factorized_) { cholmod_free_factor(&L, &c) ; factorized_ = false ; L = NULL ; } } private: int N ; bool factorized_ ; cholmod_common c ; cholmod_factor* L ; } ; #endif