Commit 01fc65b2 authored by David Bommes's avatar David Bommes
Browse files

added UMFPACKSolver and SparseQRSolver (both from SuiteSparse)

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@51 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 2cf18d25
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
//=============================================================================
//
// CLASS CholmodSolver
//
//=============================================================================
#ifndef COMISO_SPARSE_QR_SOLVER_HH
#define COMISO_SPARSE_QR_SOLVER_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "GMM_Tools.hh"
#include <CoMISo/Utils/StopWatch.hh>
#include <iostream>
#include <vector>
#include "SuiteSparseQR.hpp"
#include "cholmod.h"
//typedef struct cholmod_common_struct cholmod_common;
//== NAMESPACES ===============================================================
namespace COMISO {
typedef struct SuiteSparseQR_factorization<double> sparseqr_factor;
//== CLASS DEFINITION =========================================================
class COMISODLLEXPORT SparseQRSolver
{
public:
typedef UF_long Int;
SparseQRSolver();
~SparseQRSolver();
bool calc_system( const std::vector<Int>& _colptr,
const std::vector<Int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool calc_system_gmm( const GMM_MatrixT& _mat);
bool update_system( const std::vector<Int>& _colptr,
const std::vector<Int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool update_system_gmm( const GMM_MatrixT& _mat);
bool solve ( double * _x0, double * _b);
bool solve ( std::vector<double>& _x0, const std::vector<double>& _b)
{return solve( &(_x0[0]), (double*)&(_b[0]));}
double tolerance() { return tolerance_; }
int ordering() { return ordering_; }
bool& show_timings() { return show_timings_;}
// factorize _A*P = _Q*_R and return rank
template< class GMM_MatrixT, class GMM_MatrixT2, class GMM_MatrixT3, class IntT>
int factorize_system_gmm( const GMM_MatrixT& _A, GMM_MatrixT2& _Q, GMM_MatrixT3& _R, std::vector<IntT>& _P);
private:
cholmod_common * mp_cholmodCommon;
sparseqr_factor * mp_F;
double tolerance_;
// specify Permutation ordering
// 1. user-provided permutation (only for cholmod analyze p).
// 2. AMD with default parameters.
// 3. METIS with default parameters.
// 4. NESDIS with default parameters
// 5. natural ordering (with weighted postorder).
// 6. NESDIS, nd small = 20000, prune dense = 10.
// 7. NESDIS, nd small = 4, prune dense = 10, no constrained minimum degree.
// 8. NESDIS, nd small = 200, prune dense = 0.
// 9. COLAMD for A*A^T or AMD for A
int ordering_;
// dimension of the mxn matrix
int m_;
int n_;
std::vector<double> values_;
std::vector<Int> colptr_;
std::vector<Int> rowind_;
bool show_timings_;
StopWatch sw_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_SPARSE_QR_SOLVER_TEMPLATES_C)
#define COMISO_SPARSE_QR_SOLVER_TEMPLATES
#include "SparseQRSolverT.cc"
#endif
//=============================================================================
#endif // COMISO_SPARSE_QR_SOLVER_HH defined
//=============================================================================
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#define COMISO_SPARSE_QR_SOLVER_TEMPLATES_C
#include "SparseQRSolver.hh"
namespace COMISO {
template< class GMM_MatrixT>
bool SparseQRSolver::calc_system_gmm( const GMM_MatrixT& _mat)
{
// std::vector<int> colptr;
// std::vector<int> rowind;
// std::vector<double> values;
if(show_timings_) sw_.start();
COMISO_GMM::get_ccs_symmetric_data( _mat,
'c',
values_,
rowind_,
colptr_ );
if(show_timings_)
{
std::cerr << "Cholmod Timing GMM convert: " << sw_.stop()/1000.0 << "s\n";
}
return calc_system( colptr_, rowind_, values_);
}
//-----------------------------------------------------------------------------
template< class GMM_MatrixT>
bool SparseQRSolver::update_system_gmm( const GMM_MatrixT& _mat)
{
// std::vector<int> colptr;
// std::vector<int> rowind;
// std::vector<double> values;
COMISO_GMM::get_ccs_symmetric_data( _mat,
'c',
values_,
rowind_,
colptr_ );
return update_system( colptr_, rowind_, values_);
}
//-----------------------------------------------------------------------------
template< class GMM_MatrixT, class GMM_MatrixT2, class GMM_MatrixT3, class IntT>
int
SparseQRSolver::
factorize_system_gmm( const GMM_MatrixT& _A, GMM_MatrixT2& _Q, GMM_MatrixT3& _R, std::vector<IntT>& _P)
{
std::cerr << "factorize_system_gmm" << std::endl;
// get dimensions
int m = gmm::mat_nrows(_A);
int n = gmm::mat_ncols(_A);
// 1. _A -> cholmod_sparse A
cholmod_sparse* AC(0);
COMISO_GMM::gmm_to_cholmod(_A, AC, mp_cholmodCommon, 0, true);
std::cerr << "gmm_to_cholmod finished" << std::endl;
cholmod_print_sparse(AC, "AC", mp_cholmodCommon);
// 2. factorize A -> Q,R,P
UF_long econ = m;
cholmod_sparse *Q, *R;
// UF_long *P = new UF_long[n];
UF_long *P;
double rank = SuiteSparseQR<double>(ordering_, tolerance_, econ, AC, &Q, &R, &P, mp_cholmodCommon);
std::cerr << "factorization finished" << std::endl;
std::cerr << "rank: " << rank << std::endl;
cholmod_print_sparse(Q, "Q", mp_cholmodCommon);
// 3. convert Q,R,P -> _Q, _R, _P
COMISO_GMM::cholmod_to_gmm(*Q, _Q);
COMISO_GMM::cholmod_to_gmm(*R, _R);
std::cerr << "cholmod_to_gmm finished" << std::endl;
_P.clear(); _P.resize(n);
for( int i=0; i<n; ++i)
_P[i] = P[i];
std::cerr << "coy vector finished" << std::endl;
cholmod_l_free_sparse(&Q, mp_cholmodCommon);
cholmod_l_free_sparse(&R, mp_cholmodCommon);
std::cerr << "free1 finished" << std::endl;
// TODO: alloc or free P ???
cholmod_free(n, sizeof(UF_long), P, mp_cholmodCommon);
std::cerr << "free2 finished" << std::endl;
//// [Q,R,E] = qr(A), returning Q as a sparse matrix
//template <typename Entry> UF_long SuiteSparseQR // returns rank(A) estimate
//(
// int ordering, // all, except 3:given treated as 0:fixed
// double tol,
// UF_long econ,
// cholmod_sparse *A, // m-by-n sparse matrix
// // outputs
// cholmod_sparse **Q, // m-by-e sparse matrix where e=max(econ,rank(A))
// cholmod_sparse **R, // e-by-n sparse matrix
// UF_long **E, // permutation of 0:n-1, NULL if identity
// cholmod_common *cc // workspace and parameters
//) ;
std::cerr << " ############## QR Factorization Info #############\n";
std::cerr << " m: " << m << ", n: " << n << ", rank: " << rank << std::endl;
return rank;
}
}
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include "UMFPACKSolver.hh"
namespace COMISO {
UMFPACKSolver::UMFPACKSolver()
{
// initialize zero pointers
symbolic_ = 0;
numeric_ = 0;
tolerance_ = 1e-8;
// ordering_ = CHOLMOD_AMD;
show_timings_ = false;
}
//-----------------------------------------------------------------------------
UMFPACKSolver::~UMFPACKSolver()
{
if( symbolic_ )
{
umfpack_di_free_symbolic( &symbolic_);
symbolic_ = 0;
}
if( numeric_ )
{
umfpack_di_free_numeric( &numeric_);
numeric_ = 0;
}
}
//-----------------------------------------------------------------------------
bool UMFPACKSolver::calc_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values)
{
if(show_timings_) sw_.start();
colptr_ = _colptr;
rowind_ = _rowind;
values_ = _values;
int n = colptr_.size()-1;
// clean up
if( symbolic_ )
{
umfpack_di_free_symbolic( &symbolic_);
symbolic_ = 0;
}
if( numeric_ )
{
umfpack_di_free_numeric( &numeric_);
numeric_ = 0;
}
if(show_timings_)
{
std::cerr << " UMFPACK Timing cleanup: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
int status;
// symbolic factorization
status = umfpack_di_symbolic(n,n,&(colptr_[0]), &(rowind_[0]), &(values_[0]), &symbolic_, 0, 0);
if( status != UMFPACK_OK )
{
std::cout << "UMFPACK_symbolic failed" << std::endl;
print_error( status);
return false;
}
if(show_timings_)
{
std::cerr << " SuiteSparseQR_symbolic Timing: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
// numeric factorization
status = umfpack_di_numeric ( &(colptr_[0]), &(rowind_[0]), &(values_[0]), symbolic_, &numeric_, 0, 0);
if( status != UMFPACK_OK )
{
std::cout << "UMFPACK_numeric failed" << std::endl;
print_error( status);
return false;
}
if(show_timings_)
{
std::cerr << " SuiteSparseQR_numeric Timing: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
return true;
}
//-----------------------------------------------------------------------------
void UMFPACKSolver::print_error( int _status )
{
switch(_status)
{
case UMFPACK_OK : std::cerr << "UMFPACK Error: UMFPACK_OK\n";
case UMFPACK_ERROR_n_nonpositive : std::cerr << "UMFPACK Error: UMFPACK_ERROR_n_nonpositive\n";
case UMFPACK_ERROR_invalid_matrix : std::cerr << "UMFPACK Error: UMFPACK_ERROR_invalid_matrix\n";
case UMFPACK_ERROR_out_of_memory : std::cerr << "UMFPACK Error: UMFPACK_ERROR_out_of_memory\n";
case UMFPACK_ERROR_argument_missing: std::cerr << "UMFPACK Error: UMFPACK_ERROR_argument_missing\n";
case UMFPACK_ERROR_internal_error : std::cerr << "UMFPACK Error: UMFPACK_ERROR_internal_error\n";
default: std::cerr << "UMFPACK Error: UNSPECIFIED ERROR\n";
}
}
//-----------------------------------------------------------------------------
bool UMFPACKSolver::update_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values )
{
if( !symbolic_ )
return false;
if(show_timings_) sw_.start();
colptr_ = _colptr;
rowind_ = _rowind;
values_ = _values;
int n = colptr_.size()-1;
if( numeric_ )
{
umfpack_di_free_numeric( &numeric_);
numeric_ = 0;
}
if(show_timings_)
{
std::cerr << " UMFPACK Timing cleanup: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
int status;
// numeric factorization
status = umfpack_di_numeric ( &(colptr_[0]), &(rowind_[0]), &(values_[0]), symbolic_, &numeric_, 0, 0);
if( status != UMFPACK_OK )
{
std::cout << "UMFPACK_numeric failed" << std::endl;
print_error( status);
return false;
}
if(show_timings_)
{
std::cerr << " SuiteSparseQR_numeric Timing: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
return true;
}
//-----------------------------------------------------------------------------
bool UMFPACKSolver::solve( double * _x, double * _b)
{
if(show_timings_) sw_.start();
// const unsigned int n = colptr_.size() - 1;
int status = umfpack_di_solve( UMFPACK_A, &(colptr_[0]), &(rowind_[0]), &(values_[0]),
_x, _b, numeric_, 0, 0);
if( status != UMFPACK_OK )
{
std::cout << "UMFPACK_solve failed" << std::endl;
print_error( status);
return false;
}
if(show_timings_)
{
std::cerr << " UMFPACK_sove Timing: " << sw_.stop()/1000.0 << "s\n";
sw_.start();
}
return true;
}
//-----------------------------------------------------------------------------
}
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
//=============================================================================
//
// CLASS UMPFACKSolver
//
//=============================================================================
#ifndef COMISO_UMFPACK_SOLVER_HH
#define COMISO_UMFPACK_SOLVER_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "GMM_Tools.hh"
#include <CoMISo/Utils/StopWatch.hh>
#include <iostream>
#include <vector>
#include "umfpack.h"
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
class COMISODLLEXPORT UMFPACKSolver
{
public:
// _size is maximal size this instance can handle (smaller problems are possible!!!)
UMFPACKSolver();
~UMFPACKSolver();
bool calc_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool calc_system_gmm( const GMM_MatrixT& _mat);
bool update_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool update_system_gmm( const GMM_MatrixT& _mat);
bool solve ( double * _x0, double * _b);
bool solve ( std::vector<double>& _x0, std::vector<double>& _b)
{return solve( &(_x0[0]), &(_b[0]));}
double tolerance() { return tolerance_; }
int ordering() { return ordering_; }