- added Taucs Solver

- changed cmake-system of examples

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@83 1355f012-dd97-4b2f-ae87-10fa9f823a57

CoMISo/CMakeLists.txt

@@ -111,6 +111,15 @@ else ()
   set (COMISO_Eigen3_CONFIG_FILE_SETTINGS "#define COMISO_Eigen3_AVAILABLE 0" )
 endif ()
 
+find_package (Taucs)
+if (TAUCS_FOUND )
+  set (COMISO_TAUCS_CONFIG_FILE_SETTINGS "#define COMISO_TAUCS_AVAILABLE 1" )
+  list( APPEND COMISO_INCLUDE_DIRECTORIES ${TAUCS_INCLUDE_DIR} )
+else ()
+  message (STATUS "TAUCS not found!")
+  set (COMISO_TAUCS_CONFIG_FILE_SETTINGS "#define COMISO_TAUCS_AVAILABLE 0" )
+endif ()
+
 include_directories (
   ..
   ${CMAKE_CURRENT_SOURCE_DIR}

CoMISo/Config/config.hh.in

@@ -8,5 +8,6 @@
 @COMISO_TAO_CONFIG_FILE_SETTINGS@
 @COMISO_IPOPT_CONFIG_FILE_SETTINGS@
 @COMISO_MUMPS_CONFIG_FILE_SETTINGS@
+@COMISO_TAUCS_CONFIG_FILE_SETTINGS@
 
 

CoMISo/Examples/factored_solver/CMakeLists.txt

 include (ACGCommon)
-find_package(CoMISo)
-find_package(SUITESPARSE)
-find_package(MPI)
-find_package (MUMPS)
-find_package(IPOPT)
-
-if (MUMPS_FOUND )
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${MUMPS_LIBRARY} )
-endif ()
-
-find_package (IPOPT)
-if (IPOPT_FOUND AND MUMPS_FOUND)
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_DIRECTORIES ${IPOPT_LIBRARY_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${IPOPT_LIBRARY} )
-endif ()
-
-include_directories (
-  ..
-  ${CMAKE_SOURCE_DIR}
-  ${CMAKE_CURRENT_SOURCE_DIR}
-  ${CMAKE_CURRENT_BINARY_DIR}
-  ${COMISO_INCLUDE_DIR}
-  ${SUITESPARSE_INCLUDE_DIRS}
-  ${COMISO_INCLUDE_DIRECTORIES}
-)
-
-link_directories (
-  ${SUITESPARSE_LIBRARY_DIRS}
-  ${TAO_LIBRARY_DIR}
-  ${PETSC_LIBRARY_DIR}
-  ${COMISO_LINK_DIRECTORIES}
-)
+include (CoMISoExample)
 
 # source code directories
 set (directories 
@@ -60,10 +27,6 @@ set_target_properties(factored_solver PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
 
 target_link_libraries (factored_solver
   CoMISo
-  ${SUITESPARSE_LIBRARIES}
-  ${MPI_LIBRARY}
-  ${PETSC_LIBRARY}
-  ${TAO_LIBRARY}
   ${COMISO_LINK_LIBRARIES}
 )
 

CoMISo/Examples/quadratic_solver/CMakeLists.txt

 include (ACGCommon)
-find_package(CoMISo)
-find_package(SUITESPARSE)
-find_package(MPI)
-find_package (MUMPS)
-find_package(IPOPT)
-
-if (MUMPS_FOUND )
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${MUMPS_LIBRARY} )
-endif ()
-
-find_package (IPOPT)
-if (IPOPT_FOUND AND MUMPS_FOUND)
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_DIRECTORIES ${IPOPT_LIBRARY_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${IPOPT_LIBRARY} )
-endif ()
-
-include_directories (
-  ..
-  ${CMAKE_SOURCE_DIR}
-  ${CMAKE_CURRENT_SOURCE_DIR}
-  ${CMAKE_CURRENT_BINARY_DIR}
-  ${COMISO_INCLUDE_DIR}
-  ${SUITESPARSE_INCLUDE_DIRS}
-  ${COMISO_INCLUDE_DIRECTORIES}
-)
-
-link_directories (
-  ${SUITESPARSE_LIBRARY_DIRS}
-  ${TAO_LIBRARY_DIR}
-  ${PETSC_LIBRARY_DIR}
-  ${COMISO_LINK_DIRECTORIES}
-)
+include (CoMISoExample)
 
 # source code directories
 set (directories 
@@ -60,10 +27,6 @@ set_target_properties(quadratic_solver PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
 
 target_link_libraries (quadratic_solver
   CoMISo
-  ${SUITESPARSE_LIBRARIES}
-  ${MPI_LIBRARY}
-  ${PETSC_LIBRARY}
-  ${TAO_LIBRARY}
   ${COMISO_LINK_LIBRARIES}
 )
 

CoMISo/Examples/small_factored_example/CMakeLists.txt

 include (ACGCommon)
-find_package(CoMISo)
-find_package(SUITESPARSE)
-find_package(MPI)
-find_package (MUMPS)
-find_package(IPOPT)
-
-if (MUMPS_FOUND )
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${MUMPS_LIBRARY} )
-endif ()
-
-find_package (IPOPT)
-if (IPOPT_FOUND AND MUMPS_FOUND)
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_DIRECTORIES ${IPOPT_LIBRARY_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${IPOPT_LIBRARY} )
-endif ()
-
-include_directories (
-  ..
-  ${CMAKE_SOURCE_DIR}
-  ${CMAKE_CURRENT_SOURCE_DIR}
-  ${CMAKE_CURRENT_BINARY_DIR}
-  ${COMISO_INCLUDE_DIR}
-  ${SUITESPARSE_INCLUDE_DIRS}
-  ${COMISO_INCLUDE_DIRECTORIES}
-)
-
-link_directories (
-  ${SUITESPARSE_LIBRARY_DIRS}
-  ${TAO_LIBRARY_DIR}
-  ${PETSC_LIBRARY_DIR}
-  ${COMISO_LINK_DIRECTORIES}
-)
+include (CoMISoExample)
 
 # source code directories
 set (directories 
@@ -60,10 +27,6 @@ set_target_properties(small_factored_solver PROPERTIES INSTALL_RPATH_USE_LINK_PA
 
 target_link_libraries (small_factored_solver
   CoMISo
-  ${SUITESPARSE_LIBRARIES}
-  ${MPI_LIBRARY}
-  ${PETSC_LIBRARY}
-  ${TAO_LIBRARY}
   ${COMISO_LINK_LIBRARIES}
 )
 

CoMISo/Examples/small_quadratic_example/CMakeLists.txt

 include (ACGCommon)
-find_package(CoMISo)
-find_package(SUITESPARSE)
-find_package(MPI)
-find_package (MUMPS)
-find_package(IPOPT)
-
-if (MUMPS_FOUND )
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${MUMPS_LIBRARY} )
-endif ()
-
-find_package (IPOPT)
-if (IPOPT_FOUND AND MUMPS_FOUND)
-  list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} )
-  list( APPEND COMISO_LINK_DIRECTORIES ${IPOPT_LIBRARY_DIR} )
-  list( APPEND COMISO_LINK_LIBRARIES ${IPOPT_LIBRARY} )
-endif ()
-
-include_directories (
-  ..
-  ${CMAKE_SOURCE_DIR}
-  ${CMAKE_CURRENT_SOURCE_DIR}
-  ${CMAKE_CURRENT_BINARY_DIR}
-  ${COMISO_INCLUDE_DIR}
-  ${SUITESPARSE_INCLUDE_DIRS}
-  ${COMISO_INCLUDE_DIRECTORIES}
-)
-
-link_directories (
-  ${SUITESPARSE_LIBRARY_DIRS}
-  ${TAO_LIBRARY_DIR}
-  ${PETSC_LIBRARY_DIR}
-  ${COMISO_LINK_DIRECTORIES}
-)
+include (CoMISoExample)
 
 # source code directories
 set (directories 
@@ -60,10 +27,6 @@ set_target_properties(small_quadratic_solver PROPERTIES INSTALL_RPATH_USE_LINK_P
 
 target_link_libraries (small_quadratic_solver
   CoMISo
-  ${SUITESPARSE_LIBRARIES}
-  ${MPI_LIBRARY}
-  ${PETSC_LIBRARY}
-  ${TAO_LIBRARY}
   ${COMISO_LINK_LIBRARIES}
 )
 

CoMISo/Solver/TaucsSolver.cc

+#include "TaucsSolver.hh"
+
+
+namespace COMISO {
+
+TaucsSolver::TaucsSolver(int _size, bool _write_log)
+    : PAP_(0),
+      LSN_(0),
+      perm_(0),
+      invperm_(0),
+      px_(0),
+      pb_(0)
+  {
+    if (_write_log)
+       taucs_logfile((char*) "taucs.log");
+    resize( _size);
+  }
+
+
+  //-----------------------------------------------------------------------------
+  
+
+TaucsSolver::~TaucsSolver()
+  {
+    // free matrices
+    if(LSN_)
+    {
+      taucs_supernodal_factor_free_numeric(LSN_);
+      taucs_supernodal_factor_free( LSN_);
+    }
+
+    if(perm_)     delete [] perm_;
+    if(invperm_)  delete [] invperm_;
+    if(px_)       delete [] px_;
+    if(pb_)       delete [] pb_;
+
+  }
+  
+
+//-----------------------------------------------------------------------------
+
+
+void TaucsSolver::calc_system( const std::vector<int>&    _colptr, 
+		               const std::vector<int>&    _rowind, 
+			       const std::vector<double>& _values)
+  {
+    colptr_ = _colptr;
+    rowind_ = _rowind;
+    values_ = _values;
+
+    calc_system();
+  }
+
+//-----------------------------------------------------------------------------
+
+
+void TaucsSolver::calc_system()
+  {
+    int n   = colptr_.size()-1;
+
+    // free matrices
+    //    if(PAP_) taucs_ccs_free(PAP_);
+    if(LSN_)
+    {
+      taucs_supernodal_factor_free_numeric(LSN_);
+      taucs_supernodal_factor_free( LSN_);
+    }
+
+    mat_.n = n;
+    mat_.m = mat_.n;
+    mat_.flags = TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER;
+
+    mat_.colptr   = &colptr_[0];
+    mat_.rowind   = &rowind_[0];
+    mat_.values.d = &values_[0];
+
+    // calc permutation
+    taucs_ccs_order(&mat_, &perm_, &invperm_, (char*)"metis");
+    
+    // permute Matrix
+    PAP_ = taucs_ccs_permute_symmetrically(&mat_, perm_, invperm_);
+
+    // factor symbolic
+    LSN_ = taucs_ccs_factor_llt_symbolic( PAP_);
+
+    // factor numeric
+    taucs_ccs_factor_llt_numeric(PAP_,LSN_);
+
+    // free permuted matrix
+    taucs_ccs_free(PAP_);
+  }
+
+
+  //-----------------------------------------------------------------------------
+
+
+void TaucsSolver::update_system( const std::vector<int>&    _colptr, 
+				 const std::vector<int>&    _rowind, 
+				 const std::vector<double>& _values)
+  {
+    colptr_ = _colptr;
+    rowind_ = _rowind;
+    values_ = _values;
+    
+    update_system();
+  }
+
+
+  //-----------------------------------------------------------------------------
+
+
+void TaucsSolver::update_system()
+  {
+    int n   = colptr_.size()-1;
+
+    mat_.n = n;
+    mat_.m = mat_.n;
+    mat_.flags = TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER;
+
+    mat_.colptr   = &colptr_[0];
+    mat_.rowind   = &rowind_[0];
+    mat_.values.d = &values_[0];
+
+    // free numeric part of matrix
+    if(LSN_){taucs_supernodal_factor_free_numeric(LSN_);}
+
+    // permute Matrix
+    PAP_ = taucs_ccs_permute_symmetrically(&mat_, perm_, invperm_);
+
+    // factor numeric
+    taucs_ccs_factor_llt_numeric(PAP_,LSN_);
+
+    // free permuted matrix
+    taucs_ccs_free(PAP_);
+  }
+
+
+  //-----------------------------------------------------------------------------
+  
+
+  void TaucsSolver::solve( double * _x, double * _b)
+  {
+    // solve x coords
+    taucs_vec_permute(mat_.n, mat_.flags, _b, pb_, perm_);
+    taucs_supernodal_solve_llt(LSN_, px_, pb_);
+    taucs_vec_ipermute(mat_.n, mat_.flags, px_, _x, perm_);
+  }
+
+
+  //-----------------------------------------------------------------------------
+  
+
+  void TaucsSolver::solve_cg( double* _x0, double* _b, int _maxiter, double _max_error, bool /*_precond*/)
+  {
+//     // preconditioner
+//     taucs_ccs_matrix *M;
+//     M = taucs_ccs_factor_llt(&mat_, 1E-1, 1);
+
+//     // solve system
+//     taucs_conjugate_gradients( &mat_, taucs_ccs_solve_llt, M, _x0, _b, _maxiter, _max_error );
+
+    // solve system
+    taucs_conjugate_gradients( &mat_, 0, 0, _x0, _b, _maxiter, _max_error );
+  }
+
+
+  //-----------------------------------------------------------------------------
+
+
+  void TaucsSolver::eliminate_var( int _i, double _xi, double* _x, double* _rhs)
+  {
+    int n = colptr_.size()-1;
+
+    int iv_old(0);
+    int iv_new(0);
+
+    // for all columns
+    for(int j=0; j<(int)colptr_.size()-1; ++j)
+    {
+      // update x and rhs
+      if( j > _i)
+      {
+	_rhs[j-1] = _rhs[j];
+	_x  [j-1] = _x  [j];
+      }
+
+      if( j == _i)
+      {
+	// update rhs
+	for(int i=colptr_[j]; i<colptr_[j+1]; ++i)
+	{
+	  _rhs[rowind_[iv_old]] -= _xi*values_[iv_old];
+	  ++iv_old;
+	}
+      }
+      else
+      {
+	// store index to update colptr
+	int iv_new_save = iv_new;
+
+	for(int i=colptr_[j]; i<colptr_[j+1]; ++i)
+	{
+	  if( rowind_[iv_old] < _i)
+	  {
+	    rowind_[iv_new] = rowind_[iv_old];
+	    values_[iv_new] = values_[iv_old];
+	    ++iv_new;
+	  }
+	  else if( rowind_[iv_old] > _i)
+	  {
+	    rowind_[iv_new] = rowind_[iv_old]-1;
+	    values_[iv_new] = values_[iv_old];
+	    ++iv_new;
+	  }
+	  else
+	  {
+	    if( j< _i)
+	      _rhs[j] -= _xi*values_[iv_old];
+	    else
+	      _rhs[j-1] -= _xi*values_[iv_old];
+
+	  }
+	  ++iv_old;
+	}
+	
+	if( j<_i)
+	  colptr_[j] = iv_new_save;
+	else
+	  if( j>_i)
+	    colptr_[j-1] = iv_new_save;
+      }
+    }
+
+    // store index to end
+    colptr_[colptr_.size()-2] = iv_new;
+    // resize vectors
+    colptr_.resize( colptr_.size()-1);
+    values_.resize( iv_new);
+    rowind_.resize( iv_new);
+
+    mat_.n = n-1;
+    mat_.m = n-1;
+  }
+
+
+  //-----------------------------------------------------------------------------
+
+
+  void TaucsSolver::resize( int _size)
+  {
+    if(perm_)     delete [] perm_;
+    if(invperm_)  delete [] invperm_;
+    if(px_)       delete [] px_;
+    if(pb_)       delete [] pb_;
+
+    perm_    = new int[_size];
+    invperm_ = new int[_size];
+    px_      = new double[_size];
+    pb_      = new double[_size];
+  }
+
+  
+}

CoMISo/Solver/TaucsSolver.hh

+//=============================================================================
+//
+//  CLASS TaucsSolver
+//
+//  Author:  David Bommes <bommes@cs.rwth-aachen.de>
+//
+//  Version: $Revision: 1$
+//  Date:    $Date: 02-04-200X$
+//
+//=============================================================================
+
+#ifndef COMISO_TAUCS_SOLVER_HH
+#define COMISO_TAUCS_SOLVER_HH
+
+
+//== INCLUDES =================================================================
+#include <CoMISo/Config/CoMISoDefines.hh>
+
+
+#ifdef WIN32
+        extern "C"
+        {
+                #include <taucs.h>
+        }
+        #undef min
+        #undef max
+#else
+        #include <taucs.h>
+#endif
+
+
+#include "GMM_Tools.hh"
+#include <iostream>
+#include <vector>
+
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+class COMISODLLEXPORT TaucsSolver
+  {
+  public:
+
+    // _size is maximal size this instance can handle (smaller problems are possible!!!)
+    TaucsSolver(int _size, bool _write_log = false);
+    ~TaucsSolver();
+    
+    void calc_system( const std::vector<int>&    _colptr, 
+		      const std::vector<int>&    _rowind, 
+		      const std::vector<double>& _values );
+
+
+    template< class GMM_MatrixT>
+    void calc_system_gmm( const GMM_MatrixT& _mat);
+
+
+    void update_system( const std::vector<int>&    _colptr, 
+			const std::vector<int>&    _rowind, 
+			const std::vector<double>& _values );
+
+
+    template< class GMM_MatrixT>
+    void update_system_gmm( const GMM_MatrixT& _mat);
+
+
+    void solve( std::vector<double>& _x0, std::vector<double>& _b) { solve( &(_x0[0]), &(_b[0])); }
+    void solve( double* _x0, double* _b);
+
+    void resize       ( int _size);
+
+    // CG support
+
+    template< class GMM_MatrixT>
+    void init_system_gmm( const GMM_MatrixT& _mat);
+
+    void solve_cg( std::vector<double>& _x0, std::vector<double>& _b, int _maxiter = 50, double _max_error=1e-6, bool _precond=true) { solve_cg( &(_x0[0]), &(_b[0]), _maxiter, _max_error, _precond); }
+    void solve_cg( double* _x0, double* _b, int _maxiter = 50, double _max_error=1e-6, bool _precond=true);
+    
+    void eliminate_var( int _i, double _xi, std::vector<double>& _x, std::vector<double>& _rhs) {eliminate_var(_i, _xi, &(_x[0]), &(_rhs[0])); _rhs.resize(_rhs.size()-1); _x.resize(_x.size()-1);}
+    void eliminate_var( int _i, double _xi, double* _x, double* _rhs);
+
+    template< class GMM_MatrixT>
+    void get_matrix_gmm( GMM_MatrixT& _mat);
+
+
+    void calc_system();
+    void update_system();
+
+  private:
+
+    taucs_ccs_matrix   mat_;
+    taucs_ccs_matrix  *PAP_;
+    
+    void              *LSN_;
+
+    int               *perm_;
+    int               *invperm_;
+
+    double            *px_;
+    double            *pb_;
+
+
+    std::vector<double> values_;
+    std::vector<int>    rowind_;
+    std::vector<int>    colptr_;
+  };
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_TAUCS_SOLVER_TEMPLATES_C)
+#define COMISO_TAUCS_SOLVER_TEMPLATES
+#include "TaucsSolverT.cc"
+#endif
+//=============================================================================
+#endif // COMISO_TAUCS_SOLVER_HH defined
+//=============================================================================

CoMISo/Solver/TaucsSolverT.cc

+#define COMISO_TAUCS_SOLVER_TEMPLATES_C
+
+#include "TaucsSolver.hh"
+