Skip to content
GitLab
Commit 19de238c authored by David Bommes's avatar David Bommes
Browse files

added update/downdate functionality

parent 3d7d18a9
Pipeline #2911 passed with stage
in 5 minutes and 49 seconds
Hide whitespace changes
Inline Side-by-side
Solver/CholmodSolver.cc
View file @ 19de238c
......@@ -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
View file @ 19de238c
......@@ -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
View file @ 19de238c
......@@ -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 0 → 100644
View file @ 19de238c
# - 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)
Supports Markdown
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