some updates in openNL linear solver

SCHNApps/CMakeLists.txt

@@ -17,6 +17,7 @@ find_package(GLEW REQUIRED)
 find_package(Qt4 REQUIRED)
 find_package(QGLViewer REQUIRED)
 find_package(PythonLibs 2.7 REQUIRED)
+find_package(SuiteSparse REQUIRED)
 
 SET( QT_USE_QTOPENGL TRUE )
 SET( QT_USE_QTXML TRUE )
@@ -66,6 +67,7 @@ SET (EXT_INCLUDES
 	${QT_INCLUDE_DIR}
 	${QGLVIEWER_INCLUDE_DIR}
 	${PYTHON_INCLUDE_DIRS}
+	${SUITESPARSE_INCLUDE_DIRS}
 )
 
 # define libs for external libs
@@ -81,6 +83,7 @@ SET (EXT_LIBS
 	${QT_LIBRARIES}
 	${QGLVIEWER_LIBRARIES}
 	${PYTHON_LIBRARIES}
+	${SUITESPARSE_LIBRARIES}
 )
 
 

SCHNApps/bin/init.py

@@ -5,7 +5,7 @@ differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
 subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
 surfaceDeformationPlugin = schnapps.loadPlugin("SurfaceDeformation");
 
-obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/midRes/camel_10k.off");
+obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
 
 v = schnapps.getView("view_0");
 

SCHNApps/src/main.cpp

-
 #include <QSplashScreen>
 #include "window.h"
 #include <QFileInfo>
 #include "PythonQt/PythonQt.h"
 #include "PythonQt/gui/PythonQtScriptingConsole.h"
 
-
 int main(int argc, char* argv[])
 {
 	QApplication app(argc, argv);

ThirdParty/include/OpenNL/cholesky_solver.h

@@ -32,7 +32,7 @@
 #define __CHOLESKY_SOLVER__
 
 #include <cassert>
-#include "cholmod.h"
+#include <cholmod.h>
 
 template <typename CoeffType>
 class Solver_CHOLESKY
@@ -58,43 +58,41 @@ public:
 //		assert(A_in.has_symmetric_storage()) ;
 
 		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);
+		int NNZ = A_in.nonZeros() ;
+		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 ;
-
-		cholmod_sparse A = Eigen::viewAsCholmod(A_in);
+		int count = 0 ;
+		for(int j = 0; j < N; j++) {
+			colptr[j] = count;
+			for(typename Eigen::SparseMatrix<CoeffType>::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) ;
+		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 Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& b_in, Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& x_out) {
     	assert(factorized_) ;
-        assert(L->n == b_in.dimension()) ;
-        assert(L->n == x_out.dimension()) ;
+		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) ;
@@ -111,8 +109,8 @@ public:
 			x_out[i] = cx[i] ;
 
 		// Cleanup
-		cholmod_free_dense(&x, &c) ;
 		cholmod_free_dense(&b, &c) ;
+		cholmod_free_dense(&x, &c) ;
     }
 
     void reset() {

ThirdParty/include/OpenNL/linear_solver.h

@@ -36,10 +36,8 @@
 #include <cassert>
 #include <cstdlib>
 
-//#include "cholesky_solver.h"
-
-#include <Eigen/Cholesky>
 #include <Eigen/Sparse>
+#include "cholesky_solver.h"
 
 
 template <class CoeffType>
@@ -83,7 +81,7 @@ public:
 		delete A_ ;
 		delete x_ ;
 		delete b_ ;
-//		delete direct_solver_ ;
+		delete direct_solver_ ;
 	}
 
 	// __________________ Parameters ________________________
@@ -191,13 +189,7 @@ private:
 	Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* x_ ;
 	Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* b_ ;
 
-//	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_;
+	Solver_CHOLESKY<CoeffType>* direct_solver_ ;
 } ;
 
 #include "OpenNL/linear_solver.hpp"

ThirdParty/include/OpenNL/linear_solver.hpp

@@ -33,9 +33,7 @@ LinearSolver<CoeffType>::LinearSolver(unsigned int nb_variables) {
 	A_ = NULL ;
 	x_ = NULL ;
 	b_ = NULL ;
-//	symmetric_solver_ = NULL;
-//	nonsymmetric_solver_ = NULL;
-	direct_solver_ = NULL;
+	direct_solver_ = new Solver_CHOLESKY<CoeffType>() ;
 }
 
 template <typename CoeffType>
@@ -138,7 +136,7 @@ void LinearSolver<CoeffType>::end_row() {
 		if(!matrix_already_set_) {
 			for(unsigned int i = 0; i < nf; i++) {
 				for(unsigned int j = 0; j < nf; j++) {
-					A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j];
+					A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j] ;
 				}
 			}
 		}
@@ -156,11 +154,11 @@ void LinearSolver<CoeffType>::end_row() {
 		// Construct the matrix coefficients of the current row
 		if(!matrix_already_set_) {
 			for(unsigned int i = 0; i < nf; i++) {
-				A_->coeffRef(current_row_, if_[i]) += af_[i];
+				A_->coeffRef(current_row_, if_[i]) += af_[i] ;
 			}
 		}
 		// Construct the right-hand side of the current row
-		(*b_)[current_row_] = bk_ ;
+		(*b_)[current_row_] = -bk_ ;
 		for(unsigned int i = 0; i < nl; i++) {
 			(*b_)[current_row_] -= al_[i] * xl_[i] ;
 		}
@@ -171,13 +169,10 @@ void LinearSolver<CoeffType>::end_row() {
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::end_system() {
-	if(least_squares_ && direct_ && !direct_solver_)
-		direct_solver_ = new Eigen::LDLT<Eigen::SparseMatrix<CoeffType> >(*A_);
+	if(least_squares_ && direct_ && !direct_solver_->factorized())
+		direct_solver_->factorize(*A_) ;
 
-//	if(least_squares_ && direct_ && !direct_solver_->factorized())
-//		direct_solver_->factorize(*A_);
-
-    transition(IN_SYSTEM, CONSTRUCTED) ;
+	transition(IN_SYSTEM, CONSTRUCTED) ;
 }
 
 template <typename CoeffType>
@@ -186,45 +181,48 @@ void LinearSolver<CoeffType>::solve() {
 	if(least_squares_)
 	{
 		if(direct_) {
-			*x_ = direct_solver_->solve(*b_) ;
+			direct_solver_->solve(*b_, *x_) ;
 		} else {
-			Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> > solver(*A_) ;
+			Eigen::ConjugateGradient<Eigen::MatrixXf > solver(*A_) ;
+//			int n = x_->rows() ;
+//			Eigen::Matrix<CoeffType, Eigen::Dynamic, 1> guess(n) ;
+//			for(int i = 0; i < n; ++i)
+//				guess[i] = (*x_)[i] ;
+//			*x_ = solver.solveWithGuess(*b_, guess) ;
 			*x_ = solver.solve(*b_) ;
 		}
 	} else {
-		Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> > solver(*A_) ;
+		Eigen::BiCGSTAB<Eigen::MatrixXf > solver(*A_) ;
 		*x_ = solver.solve(*b_) ;
 	}
+
 	vector_to_variables() ;
 	transition(CONSTRUCTED, SOLVED) ;
 }
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::reset(bool keep_matrix) {
-    if(keep_matrix) {
-    	x_->setZero();
-    	b_->setZero();
-    	matrix_already_set_ = true ;
-    } else {
-    	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(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() ;
-    	}
-    }
+	if(keep_matrix) {
+		x_->setZero() ;
+		b_->setZero() ;
+		matrix_already_set_ = true ;
+	} else {
+		delete A_ ; A_ = NULL ;
+		delete x_ ; x_ = NULL ;
+		delete b_ ; b_ = NULL ;
+		direct_solver_->reset();
+		matrix_already_set_ = false ;
+		for(unsigned int i = 0; i < nb_variables_; ++i) {
+			variable_[i].unlock() ;
+		}
+	}
 	state_ = INITIAL ;
 }
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::vector_to_variables() {
 	for(unsigned int i=0; i < nb_variables(); i++) {
-	Variable<CoeffType>& v = variable(i) ;
+		Variable<CoeffType>& v = variable(i) ;
 		if(!v.is_locked()) {
 			v.set_value((*x_)[v.index()]) ;
 		}
@@ -234,7 +232,7 @@ void LinearSolver<CoeffType>::vector_to_variables() {
 template <typename CoeffType>
 void LinearSolver<CoeffType>::variables_to_vector() {
 	for(unsigned int i=0; i < nb_variables(); i++) {
-	Variable<CoeffType>& v = variable(i) ;
+		Variable<CoeffType>& v = variable(i) ;
 		if(!v.is_locked()) {
 			(*x_)[v.index()] = v.value() ;
 		}