Merge branch 'improveCholmod' into 'master'

added update/downdate functionality



See merge request !14

Solver/CholmodSolver.cc

@@ -151,6 +151,124 @@ bool CholmodSolver::calc_system( const std::vector<int>&    _colptr,
 }
 
 
+//-----------------------------------------------------------------------------
+
+
+bool CholmodSolver::calc_system_prepare_pattern( const std::vector<int>&    _colptr,
+                                                 const std::vector<int>&    _rowind,
+                                                 const std::vector<double>& _values,
+                                                 const std::vector<int>&    _colptr2,
+                                                 const std::vector<int>&    _rowind2,
+                                                 const std::vector<double>& _values2 )
+{
+    if(show_timings_) sw_.start();
+
+    colptr_ = _colptr;
+    rowind_ = _rowind;
+    values_ = _values;
+
+    int n   = colptr_.size()-1;
+
+    // setup matrix matA
+    cholmod_sparse matA;
+
+    matA.nrow = n;
+    matA.ncol = n;
+    matA.nzmax = _values.size();
+
+    matA.p = &colptr_[0];
+    matA.i = &rowind_[0];
+    matA.x = &values_[0];
+    matA.nz = 0;
+    matA.z = 0;
+
+    //    matA.stype = -1;
+    matA.stype = 1;
+    matA.itype = CHOLMOD_INT;
+    matA.xtype = CHOLMOD_REAL;
+    matA.dtype = CHOLMOD_DOUBLE;
+    matA.sorted = 1;
+    matA.packed = 1;
+
+    // setup matrix matA_pattern
+    cholmod_sparse matA_pattern;
+
+    matA_pattern.nrow = n;
+    matA_pattern.ncol = n;
+    matA_pattern.nzmax = _values2.size();
+
+    matA_pattern.p = (int*)(&_colptr2[0]);
+    matA_pattern.i = (int*)(&_rowind2[0]);
+    matA_pattern.x = (double*)(&_values2[0]);
+    matA_pattern.nz = 0;
+    matA_pattern.z = 0;
+
+    //    matA_pattern.stype = -1;
+    matA_pattern.stype = 1;
+    matA_pattern.itype = CHOLMOD_INT;
+    matA_pattern.xtype = CHOLMOD_REAL;
+    matA_pattern.dtype = CHOLMOD_DOUBLE;
+    matA_pattern.sorted = 1;
+    matA_pattern.packed = 1;
+
+
+    // clean up
+    if( mp_L )
+    {
+        cholmod_free_factor( &mp_L, mp_cholmodCommon );
+        mp_L = 0;
+    }
+
+
+    // compute permutation based on full pattern
+    std::vector<int> perm(matA_pattern.nrow);
+//    cholmod_metis(&matA_pattern, 0, 0, 0, perm.data(), mp_cholmodCommon) ;
+    cholmod_amd(&matA_pattern, 0, 0, perm.data(), mp_cholmodCommon) ;
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod AMD ordering: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod Timing cleanup: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+    // fix permutation to given one
+    mp_cholmodCommon->nmethods = 1;
+    mp_cholmodCommon->method[0].ordering = CHOLMOD_GIVEN;
+
+    if( !(mp_L = cholmod_analyze_p( &matA, perm.data(), 0, 0, mp_cholmodCommon )) )
+    {
+        std::cout << "cholmod_analyze failed" << std::endl;
+        return false;
+    }
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod Timing analyze: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+    if( !cholmod_factorize( &matA, mp_L, mp_cholmodCommon ) )
+    {
+        std::cout << "cholmod_factorize failed" << std::endl;
+        return false;
+    }
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod Timing factorize: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+    return true;
+}
+
   //-----------------------------------------------------------------------------
 
     
@@ -197,16 +315,82 @@ bool CholmodSolver::update_system( const std::vector<int>& _colptr,
 }
 
 
+//-----------------------------------------------------------------------------
+
+
+bool CholmodSolver::update_downdate_factor( const std::vector<int>&    _colptr,
+                                            const std::vector<int>&    _rowind,
+                                            const std::vector<double>& _values,
+                                            const bool                 _upd)
+{
+    if(show_timings_) sw_.start();
+
+    colptr_ = _colptr;
+    rowind_ = _rowind;
+    values_ = _values;
+
+    int n   = colptr_.size()-1;
+
+    cholmod_sparse matA;
+
+    matA.nrow = mp_L->n;
+    matA.ncol = n;
+    matA.nzmax = _values.size();
+
+    matA.p = &colptr_[0];
+    matA.i = &rowind_[0];
+    matA.x = &values_[0];
+    matA.nz = 0;
+    matA.z = 0;
+
+    matA.stype = 0;
+    matA.itype = CHOLMOD_INT;
+    matA.xtype = CHOLMOD_REAL;
+    matA.dtype = CHOLMOD_DOUBLE;
+    matA.sorted = 1;
+    matA.packed = 1;
+
+    // get permuted matrix
+    cholmod_sparse* matAp = cholmod_submatrix ( &matA, (int*)mp_L->Perm, mp_L->n, 0, -1, true, true, mp_cholmodCommon);
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod conversion Timing: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+    if( !cholmod_updown( _upd, matAp, mp_L, mp_cholmodCommon))
+    {
+      std::cerr << "Warning: Cholmod update/downdate failed!" << std::endl;
+    }
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod update/downdate Timing: " << sw_.stop()/1000.0 << "s\n";
+      sw_.start();
+    }
+
+    // clean up permuted matrix
+    if(matAp)
+      cholmod_free_sparse(&matAp, mp_cholmodCommon);
+
+    if(show_timings_)
+    {
+      std::cerr << " Cholmod free sparse Timing: " << sw_.stop()/1000.0 << "s\n";
+    }
+
+    return true;
+}
+
+
 //-----------------------------------------------------------------------------
   
 
 bool CholmodSolver::solve( double * _x, double * _b)
 {
-    const unsigned int n = colptr_.size() - 1;
+    const unsigned int n = mp_L->n;
 
     cholmod_dense *x, b;
 
-
     b.nrow = n;
     b.ncol = 1;
     b.nzmax = n;

Solver/CholmodSolver.hh

@@ -46,6 +46,7 @@
 #include <vector>
 
 #include "cholmod.h"
+//#include "cholmod_matrixops.h"
 
 // typedef struct cholmod_common_struct cholmod_common;
 // typedef struct cholmod_factor_struct cholmod_factor;
@@ -68,6 +69,13 @@ public:
 		      const std::vector<int>&    _rowind, 
 		      const std::vector<double>& _values );
 
+    bool calc_system_prepare_pattern( const std::vector<int>&    _colptr,
+                                      const std::vector<int>&    _rowind,
+                                      const std::vector<double>& _values,
+                                      const std::vector<int>&    _colptr2,
+                                      const std::vector<int>&    _rowind2,
+                                      const std::vector<double>& _values2 );
+
 
     template< class GMM_MatrixT>
     bool calc_system_gmm( const GMM_MatrixT& _mat);
@@ -75,6 +83,10 @@ public:
     template< class Eigen_MatrixT>
     bool calc_system_eigen( const Eigen_MatrixT& _mat);
 
+    // factorize matrix _mat by optimizing the permutation for the pattern of _mat_pattern
+    template< class Eigen_MatrixT>
+    bool calc_system_eigen_prepare_pattern( const Eigen_MatrixT& _mat, const Eigen_MatrixT& _mat_pattern);
+
 
     bool update_system( const std::vector<int>&    _colptr, 
  			const std::vector<int>&    _rowind, 
@@ -88,6 +100,15 @@ public:
     bool update_system_eigen( const Eigen_MatrixT& _mat);
 
 
+    bool update_downdate_factor( const std::vector<int>&    _colptr,
+                                 const std::vector<int>&    _rowind,
+                                 const std::vector<double>& _values,
+                                 const bool                 _upd);
+
+    template< class Eigen_MatrixT>
+    bool update_downdate_factor_eigen( const Eigen_MatrixT& _mat, const bool _upd = true);
+
+
     bool solve ( double *             _x0, double *             _b);
 
     bool solve ( std::vector<double>& _x0, std::vector<double>& _b);

Solver/CholmodSolverT.cc

@@ -107,7 +107,46 @@ bool CholmodSolver::calc_system_eigen( const Eigen_MatrixT& _mat)
 
     return calc_system( colptr_, rowind_, values_);
 }
-  
+
+
+//-----------------------------------------------------------------------------
+
+
+template< class Eigen_MatrixT>
+bool CholmodSolver::calc_system_eigen_prepare_pattern( const Eigen_MatrixT& _mat, const Eigen_MatrixT& _mat_pattern)
+{
+    if(show_timings_) sw_.start();
+
+#if COMISO_EIGEN3_AVAILABLE
+    COMISO_EIGEN::get_ccs_symmetric_data( _mat,
+                                         'u',
+                                         values_,
+                                         rowind_,
+                                         colptr_ );
+#endif
+
+    std::vector<double> values2;
+    std::vector<int>    colptr2;
+    std::vector<int>    rowind2;
+
+#if COMISO_EIGEN3_AVAILABLE
+    COMISO_EIGEN::get_ccs_symmetric_data( _mat_pattern,
+                                         'u',
+                                         values2,
+                                         rowind2,
+                                         colptr2 );
+#endif
+
+    if(show_timings_)
+    {
+      std::cerr << "Cholmod Timing EIGEN convert: " << sw_.stop()/1000.0 << "s\n";
+      std::cerr << "#nnz: " << values_.size() << std::endl;
+    }
+
+    return calc_system_prepare_pattern( colptr_, rowind_, values_, colptr2, rowind2, values2);
+}
+
+
 //-----------------------------------------------------------------------------
 
 template< class Eigen_MatrixT>
@@ -123,6 +162,30 @@ bool CholmodSolver::update_system_eigen( const Eigen_MatrixT& _mat)
   return update_system( colptr_, rowind_, values_);
 }
 
+//-----------------------------------------------------------------------------
+
+template< class Eigen_MatrixT>
+bool CholmodSolver::update_downdate_factor_eigen( const Eigen_MatrixT& _mat, const bool _upd)
+{
+    if(show_timings_) sw_.start();
+
+#if COMISO_EIGEN3_AVAILABLE
+    COMISO_EIGEN::get_ccs_symmetric_data( _mat,
+                                         'c',
+                                         values_,
+                                         rowind_,
+                                         colptr_ );
+#endif
+
+    if(show_timings_)
+    {
+      std::cerr << "Cholmod Timing EIGEN convert: " << sw_.stop()/1000.0 << "s\n";
+      std::cerr << "#nnz: " << values_.size() << std::endl;
+    }
+
+    return update_downdate_factor( colptr_, rowind_, values_, _upd);
+}
+
 
 }
 

cmake/FindSPECTRA.cmake

+# - Try to find SPECTRA
+# Once done this will define
+#  SPECTRA_FOUND         - System has SPECTRA
+#  SPECTRA_INCLUDE_DIRS  - The SPECTRA include directories
+
+if (SPECTRA_INCLUDE_DIR)
+  # in cache already
+  set(SPECTRA_FOUND TRUE)
+  set(SPECTRA_INCLUDE_DIRS "${SPECTRA_INCLUDE_DIR}" )
+else (SPECTRA_INCLUDE_DIR)
+
+# Check if the base path is set
+if ( NOT CMAKE_WINDOWS_LIBS_DIR )
+  # This is the base directory for windows library search used in the finders we shipp.
+  set(CMAKE_WINDOWS_LIBS_DIR "c:/libs" CACHE STRING "Default Library search dir on windows." )
+endif()
+
+if ( CMAKE_GENERATOR MATCHES "^Visual Studio 11.*Win64" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2012/x64/")
+elseif ( CMAKE_GENERATOR MATCHES "^Visual Studio 11.*" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2012/x32/")
+elseif ( CMAKE_GENERATOR MATCHES "^Visual Studio 12.*Win64" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2013/x64/")
+elseif ( CMAKE_GENERATOR MATCHES "^Visual Studio 12.*" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2013/x32/")
+elseif ( CMAKE_GENERATOR MATCHES "^Visual Studio 14.*Win64" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2015/x64/")
+elseif ( CMAKE_GENERATOR MATCHES "^Visual Studio 14.*" )
+  SET(VS_SEARCH_PATH "${CMAKE_WINDOWS_LIBS_DIR}/vs2015/x32/")
+endif()
+
+
+find_path( SPECTRA_INCLUDE_DIR 
+	   NAMES SymEigsSolver.h 
+           PATHS $ENV{SPECTRA_DIR}
+                 /usr/include/spectra
+                 /usr/local/include
+                 /usr/local/include/spectra/
+                 /opt/local/include/spectra/
+                 "${CMAKE_WINDOWS_LIBS_DIR}/general/spectra"
+                 "${CMAKE_WINDOWS_LIBS_DIR}/spectra"
+                 "${CMAKE_WINDOWS_LIBS_DIR}/spectra/include"
+		 "${CMAKE_WINDOWS_LIBS_DIR}/eigen/include"
+		  ${PROJECT_SOURCE_DIR}/MacOS/Libs/SPECTRA/include
+                  ../../External/include
+                  ${module_file_path}/../../../External/include
+          )
+
+set(SPECTRA_INCLUDE_DIRS "${SPECTRA_INCLUDE_DIR}" )
+
+
+include(FindPackageHandleStandardArgs)
+# handle the QUIETLY and REQUIRED arguments and set SPECTRA_FOUND to TRUE
+# if all listed variables are TRUE
+find_package_handle_standard_args(SPECTRA  DEFAULT_MSG
+                                  SPECTRA_INCLUDE_DIR)
+
+mark_as_advanced(SPECTRA_INCLUDE_DIR)
+
+endif(SPECTRA_INCLUDE_DIR)