Commit b03f81ee authored by Pierre Kraemer's avatar Pierre Kraemer

try some changes in Eigen cholesky solver

parent 734a59a7
......@@ -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
{
......
......@@ -3,6 +3,8 @@
#include "Algo/Selection/raySelector.h"
#include "Algo/Selection/collector.h"
#include "Utils/drawer.h"
#include <QKeyEvent>
#include <QMouseEvent>
......
......@@ -32,17 +32,12 @@
#define __CHOLESKY_SOLVER__
#include <cassert>
#include "OpenNL/sparse_matrix.h"
#include "CHOLMOD/cholmod.h"
#include "cholmod.h"
template< class MATRIX, class VECTOR >
template <typename CoeffType>
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<CoeffType>& 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<int*>(A->p) ;
int* rowind = static_cast<int*>(A->i) ;
double* a = static_cast<double*>(A->x) ;
// int* colptr = static_cast<int*>(A->p) ;
// int* rowind = static_cast<int*>(A->i) ;
// double* a = static_cast<double*>(A->x) ;
// Convert local Matrix into CHOLMOD Matrix
int count = 0 ;
for(int j = 0; j < N; j++) {
const typename SparseMatrix<CoeffType>::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<CoeffType>::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()) ;
......
......@@ -36,6 +36,9 @@
#include <cassert>
#include <cstdlib>
//#include "cholesky_solver.h"
#include <Eigen/Cholesky>
#include <Eigen/Sparse>
......@@ -188,9 +191,13 @@ private:
Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* x_ ;
Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* b_ ;
Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >* symmetric_solver_;
Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >* nonsymmetric_solver_;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >* direct_solver_;
// Solver_CHOLESKY<CoeffType>* direct_solver_;
Eigen::LDLT<Eigen::SparseMatrix<CoeffType> >* direct_solver_;
// Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >* symmetric_solver_;
// Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >* nonsymmetric_solver_;
// Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >* direct_solver_;
} ;
#include "OpenNL/linear_solver.hpp"
......
......@@ -33,8 +33,8 @@ LinearSolver<CoeffType>::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<CoeffType>::end_row() {
template <typename CoeffType>
void LinearSolver<CoeffType>::end_system() {
if(least_squares_ && direct_ && !direct_solver_)
direct_solver_ = new Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >(*A_);
direct_solver_ = new Eigen::LDLT<Eigen::SparseMatrix<CoeffType> >(*A_);
// if(least_squares_ && direct_ && !direct_solver_->factorized())
// direct_solver_->factorize(*A_);
transition(IN_SYSTEM, CONSTRUCTED) ;
}
......@@ -185,12 +188,12 @@ void LinearSolver<CoeffType>::solve() {
if(direct_) {
*x_ = direct_solver_->solve(*b_) ;
} else {
symmetric_solver_ = new Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >(*A_) ;
*x_ = symmetric_solver_->solve(*b_) ;
Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> > solver(*A_) ;
*x_ = solver->solve(*b_) ;
}
} else {
nonsymmetric_solver_ = new Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >(*A_) ;
*x_ = nonsymmetric_solver_->solve(*b_) ;
Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> > solver(*A_) ;
*x_ = solver->solve(*b_) ;
}
vector_to_variables() ;
transition(CONSTRUCTED, SOLVED) ;
......@@ -206,9 +209,10 @@ void LinearSolver<CoeffType>::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() ;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment