From b03f81ee1502558cb14c6d35f0302ebf3862eed0 Mon Sep 17 00:00:00 2001 From: Pierre Kraemer Date: Wed, 6 Feb 2013 13:49:41 +0100 Subject: [PATCH] try some changes in Eigen cholesky solver --- .../include/surfaceDeformation.h | 4 +- .../src/surfaceDeformation.cpp | 2 + ThirdParty/include/OpenNL/cholesky_solver.h | 59 +++++++++---------- ThirdParty/include/OpenNL/linear_solver.h | 13 +++- ThirdParty/include/OpenNL/linear_solver.hpp | 22 ++++--- 5 files changed, 56 insertions(+), 44 deletions(-) diff --git a/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h b/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h index cba93292..567270f9 100644 --- a/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h +++ b/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h @@ -6,12 +6,14 @@ #include "ui_surfaceDeformation.h" -#include "Utils/drawer.h" +//#include "Utils/drawer.h" using namespace CGoGN; using namespace SCHNApps; +namespace CGoGN { namespace Utils { class Drawer; } } + enum SelectionMode { diff --git a/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp b/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp index cfa97fa9..d31a3ac0 100644 --- a/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp +++ b/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp @@ -3,6 +3,8 @@ #include "Algo/Selection/raySelector.h" #include "Algo/Selection/collector.h" +#include "Utils/drawer.h" + #include #include diff --git a/ThirdParty/include/OpenNL/cholesky_solver.h b/ThirdParty/include/OpenNL/cholesky_solver.h index 3856caaa..3ba0c8bd 100644 --- a/ThirdParty/include/OpenNL/cholesky_solver.h +++ b/ThirdParty/include/OpenNL/cholesky_solver.h @@ -32,17 +32,12 @@ #define __CHOLESKY_SOLVER__ #include -#include "OpenNL/sparse_matrix.h" -#include "CHOLMOD/cholmod.h" +#include "cholmod.h" -template< class MATRIX, class VECTOR > +template class Solver_CHOLESKY { public: - typedef MATRIX Matrix ; - typedef VECTOR Vector ; - typedef typename Vector::CoeffType CoeffType ; - Solver_CHOLESKY() { N = 0 ; factorized_ = false ; @@ -58,43 +53,45 @@ public: bool factorized() { return factorized_ ; } - void factorize(const MATRIX &A_in) { - assert(A_in.n() == A_in.m()) ; - assert(A_in.has_symmetric_storage()) ; + void factorize(const Eigen::SparseMatrix& A_in) { + assert(A_in.rows() == A_in.cols()) ; +// assert(A_in.has_symmetric_storage()) ; - N = A_in.n() ; - int NNZ = A_in.nnz() ; + N = A_in.rows() ; +// int NNZ = A_in.nonZeros() ; // Translate sparse matrix into cholmod format - cholmod_sparse* A = cholmod_allocate_sparse(N, N, NNZ, false, true, -1, CHOLMOD_REAL, &c); +// 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) ; +// 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++) { - const typename SparseMatrix::Column& Cj = A_in.column(j) ; - colptr[j] = count ; - for(unsigned int ii = 0; ii < Cj.nb_coeffs(); ii++) { - a[count] = Cj.coeff(ii).a ; - rowind[count] = Cj.coeff(ii).index ; - count++ ; - } - } - colptr[N] = NNZ ; +// int count = 0 ; +// for(int j = 0; j < N; j++) { +// const typename SparseMatrix::Column& Cj = A_in.column(j) ; +// colptr[j] = count ; +// for(unsigned int ii = 0; ii < Cj.nb_coeffs(); ii++) { +// a[count] = Cj.coeff(ii).a ; +// rowind[count] = Cj.coeff(ii).index ; +// count++ ; +// } +// } +// colptr[N] = NNZ ; + + cholmod_sparse A = Eigen::viewAsCholmod(A_in); // Factorize - L = cholmod_analyze(A, &c) ; - cholmod_factorize(A, L, &c) ; + L = cholmod_analyze(&A, &c) ; + cholmod_factorize(&A, L, &c) ; factorized_ = true ; // Cleanup - cholmod_free_sparse(&A, &c) ; +// cholmod_free_sparse(&A, &c) ; } - void solve(const VECTOR& b_in, VECTOR& x_out) { + void solve(const VECTOR& b_in, VECTOR& x_out) { assert(factorized_) ; assert(L->n == b_in.dimension()) ; assert(L->n == x_out.dimension()) ; diff --git a/ThirdParty/include/OpenNL/linear_solver.h b/ThirdParty/include/OpenNL/linear_solver.h index 972dd2c7..1cf12db3 100644 --- a/ThirdParty/include/OpenNL/linear_solver.h +++ b/ThirdParty/include/OpenNL/linear_solver.h @@ -36,6 +36,9 @@ #include #include +//#include "cholesky_solver.h" + +#include #include @@ -188,9 +191,13 @@ private: Eigen::Matrix* x_ ; Eigen::Matrix* b_ ; - Eigen::ConjugateGradient >* symmetric_solver_; - Eigen::BiCGSTAB >* nonsymmetric_solver_; - Eigen::SimplicialLDLT >* direct_solver_; +// Solver_CHOLESKY* direct_solver_; + + Eigen::LDLT >* direct_solver_; + +// Eigen::ConjugateGradient >* symmetric_solver_; +// Eigen::BiCGSTAB >* nonsymmetric_solver_; +// Eigen::SimplicialLDLT >* direct_solver_; } ; #include "OpenNL/linear_solver.hpp" diff --git a/ThirdParty/include/OpenNL/linear_solver.hpp b/ThirdParty/include/OpenNL/linear_solver.hpp index bdc792b7..c5b54d61 100644 --- a/ThirdParty/include/OpenNL/linear_solver.hpp +++ b/ThirdParty/include/OpenNL/linear_solver.hpp @@ -33,8 +33,8 @@ LinearSolver::LinearSolver(unsigned int nb_variables) { A_ = NULL ; x_ = NULL ; b_ = NULL ; - symmetric_solver_ = NULL; - nonsymmetric_solver_ = NULL; +// symmetric_solver_ = NULL; +// nonsymmetric_solver_ = NULL; direct_solver_ = NULL; } @@ -172,7 +172,10 @@ void LinearSolver::end_row() { template void LinearSolver::end_system() { if(least_squares_ && direct_ && !direct_solver_) - direct_solver_ = new Eigen::SimplicialLDLT >(*A_); + direct_solver_ = new Eigen::LDLT >(*A_); + +// if(least_squares_ && direct_ && !direct_solver_->factorized()) +// direct_solver_->factorize(*A_); transition(IN_SYSTEM, CONSTRUCTED) ; } @@ -185,12 +188,12 @@ void LinearSolver::solve() { if(direct_) { *x_ = direct_solver_->solve(*b_) ; } else { - symmetric_solver_ = new Eigen::ConjugateGradient >(*A_) ; - *x_ = symmetric_solver_->solve(*b_) ; + Eigen::ConjugateGradient > solver(*A_) ; + *x_ = solver->solve(*b_) ; } } else { - nonsymmetric_solver_ = new Eigen::BiCGSTAB >(*A_) ; - *x_ = nonsymmetric_solver_->solve(*b_) ; + Eigen::BiCGSTAB > solver(*A_) ; + *x_ = solver->solve(*b_) ; } vector_to_variables() ; transition(CONSTRUCTED, SOLVED) ; @@ -206,9 +209,10 @@ void LinearSolver::reset(bool keep_matrix) { delete A_ ; A_ = NULL ; delete x_ ; x_ = NULL ; delete b_ ; b_ = NULL ; - if(symmetric_solver_) delete symmetric_solver_ ; symmetric_solver_ = NULL ; - if(nonsymmetric_solver_) delete nonsymmetric_solver_ ; nonsymmetric_solver_ = NULL ; +// if(symmetric_solver_) delete symmetric_solver_ ; symmetric_solver_ = NULL ; +// if(nonsymmetric_solver_) delete nonsymmetric_solver_ ; nonsymmetric_solver_ = NULL ; if(direct_solver_) delete direct_solver_ ; direct_solver_ = NULL ; +// direct_solver_->reset(); matrix_already_set_ = false ; for(unsigned int i = 0; i < nb_variables_; ++i) { variable_[i].unlock() ; -- GitLab