Commit e18b4ac5 authored by Martin Marinov's avatar Martin Marinov
Browse files

Fixed some CoMISo issues with the latest merge from VCI.

parent e773b81a
......@@ -11,8 +11,7 @@
#define COMISO_HSL_AVAILABLE 0
#define COMISO_PETSC_AVAILABLE 0
#define COMISO_TAO_AVAILABLE 0
#define COMISO_IPOPT_AVAILABLE 0
#define COMISO_IPOPTLEAN_AVAILABLE 1
#define COMISO_IPOPT_AVAILABLE 1
#define COMISO_MUMPS_AVAILABLE 1
#define COMISO_METIS_AVAILABLE 1
#define COMISO_CGAL_AVAILABLE 1
......
......@@ -26,6 +26,7 @@ SET(my_headers
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterface.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterfaceADOLC.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterfaceDCO.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemIPOPT.hh
${CMAKE_CURRENT_SOURCE_DIR}/NPTiming.hh
${CMAKE_CURRENT_SOURCE_DIR}/QuadraticProblem.hh
${CMAKE_CURRENT_SOURCE_DIR}/SuperSparseMatrixT.hh
......@@ -67,6 +68,7 @@ SET(my_sources
${CMAKE_CURRENT_SOURCE_DIR}/NPDerivativeChecker.cc
${CMAKE_CURRENT_SOURCE_DIR}/NPLinearConstraints.cc
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterface.cc
${CMAKE_CURRENT_SOURCE_DIR}/NProblemIPOPTc.cc
${CMAKE_CURRENT_SOURCE_DIR}/NPTiming.cc
${CMAKE_CURRENT_SOURCE_DIR}/TAOSolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/TapeIDSingleton.cc
......
......@@ -17,6 +17,7 @@ SET(my_t_impls
${CMAKE_CURRENT_SOURCE_DIR}/ConstrainedSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/EigenLDLTSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/Eigen_ToolsT.cc
${CMAKE_CURRENT_SOURCE_DIR}/GMM_ToolsT.cc
${CMAKE_CURRENT_SOURCE_DIR}/IterativeSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/SparseQRSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/MISolverT.cc
......
......@@ -30,1257 +30,9 @@
//
//=============================================================================
#define COMISO_GMM_TOOLS_C
//== INCLUDES =================================================================
#include "GMM_Tools.hh"
#define GMM_USES_LAPACK
#include <gmm/gmm_lapack_interface.h>
#include <queue>
#include <CoMISo/Utils/VSToolsT.hh>
//== NAMESPACES ===============================================================
namespace COMISO_GMM
{
//== IMPLEMENTATION ==========================================================
//-----------------------------------------------------------------------------
// too many if-statements, seems slower in most cases
// use eliminate_csc_vars2
template<class ScalarT, class VectorT, class RealT, class IntegerT>
void eliminate_csc_vars(
const std::vector<IntegerT>& _evar,
const std::vector<ScalarT>& _eval,
typename gmm::csc_matrix<RealT>& _A,
VectorT& _x,
VectorT& _rhs )
{
typedef unsigned int uint;
typedef typename gmm::csc_matrix<RealT> MatrixT;
unsigned int nc = _A.nc;
unsigned int nr = _A.nr;
unsigned int n_new = nc - _evar.size();
// modify rhs
for ( unsigned int k=0; k<_evar.size(); ++k )
{
IntegerT i = _evar[k];
// number of elements in this column
uint n_elem = _A.jc[i+1]-_A.jc[i];
uint r_idx = 0;
RealT entry = 0.0;
for( uint i_elem = _A.jc[i]; i_elem < _A.jc[i+1]; ++i_elem)
{
r_idx = _A.ir[i_elem];
entry = _A.pr[i_elem];
_rhs[r_idx] -= _eval[k]*( entry );
}
}
// sort vector
std::vector<IntegerT> evar( _evar );
std::sort( evar.begin(), evar.end() );
evar.push_back( std::numeric_limits<IntegerT>::max() );
// build subindex set and update rhs
std::vector<size_t> si( n_new );
int cur_evar_idx=0;
for ( unsigned int i=0; i<nc; ++i )
{
unsigned int next_i = evar[cur_evar_idx];
if ( i != next_i )
{
_rhs[i-cur_evar_idx] = _rhs[i];
_x [i-cur_evar_idx] = _x [i];
si [i-cur_evar_idx] = i;
}
else
++cur_evar_idx;
}
// delete last elements
_rhs.resize( n_new );
_x.resize( n_new );
// csc erasing rows and columns
uint offset(0);
uint next_evar_row(0);
uint next_evar_col(0);
uint last_jc = _A.jc[0];
uint col_offset(0);
uint offset_update(0);
uint last_evar_idx=evar.size()-2;
for( uint c = 0; c < nc; ++c)
{
uint el_span( _A.jc[c+1]-last_jc);
offset_update=0;
if( c != (uint) evar[next_evar_col] )
{
next_evar_row=0;
for( uint e = last_jc; e < _A.jc[c+1]; ++e)
{
while( _A.ir[e] > (uint) evar[next_evar_row])
{
++next_evar_row;
++offset_update;
}
_A.pr[e-offset] = _A.pr[e];
_A.ir[e-offset] = _A.ir[e]-offset_update;
if( _A.ir[e] == (uint)evar[next_evar_row])
{
++offset;
}
}
last_jc = _A.jc[c+1];
_A.jc[c+1-col_offset] = _A.jc[c+1]-offset;
}
else
{
++col_offset;
offset+=el_span;
++next_evar_col;
last_jc = _A.jc[c+1];
}
}
_A.nc = nc - evar.size()+1;
_A.nr = nr - evar.size()+1;
}
//-----------------------------------------------------------------------------
template<class ScalarT, class VectorT, class RealT, class IntegerT>
void eliminate_csc_vars2(
const std::vector<IntegerT>& _evar,
const std::vector<ScalarT>& _eval,
typename gmm::csc_matrix<RealT>& _A,
VectorT& _x,
VectorT& _rhs )
{
typedef unsigned int uint;
typedef typename gmm::csc_matrix<RealT> MatrixT;
unsigned int nc = _A.nc;
unsigned int nr = _A.nr;
unsigned int n_new = nc - _evar.size();
// modify rhs
for ( unsigned int k=0; k<_evar.size(); ++k )
{
IntegerT i = _evar[k];
// number of elements in this column
// uint n_elem = _A.jc[i+1]-_A.jc[i];
uint r_idx = 0;
RealT entry = 0.0;
for( uint i_elem = _A.jc[i]; i_elem < _A.jc[i+1]; ++i_elem)
{
r_idx = _A.ir[i_elem];
entry = _A.pr[i_elem];
_rhs[r_idx] -= _eval[k]*( entry );
}
}
// sort vector
std::vector<IntegerT> evar( _evar );
std::sort( evar.begin(), evar.end() );
evar.push_back( std::numeric_limits<IntegerT>::max() );
// build subindex set and update rhs
std::vector<size_t> si( n_new );
int cur_evar_idx=0;
for ( unsigned int i=0; i<nc; ++i )
{
unsigned int next_i = evar[cur_evar_idx];
if ( i != next_i )
{
_rhs[i-cur_evar_idx] = _rhs[i];
_x [i-cur_evar_idx] = _x [i];
si [i-cur_evar_idx] = i;
}
else
{
++cur_evar_idx;
}
}
// delete last elements
_rhs.resize( n_new );
_x.resize( n_new );
std::vector< int > new_row_idx_map( nr, -1);
// build re-indexing map
// -1 means deleted
{
int next_evar_idx(0);
int offset(0);
for( int i = 0; i < (int)nr; ++i)
{
if( i == (int)evar[next_evar_idx])
{
++next_evar_idx;
++offset;
}
else
{
new_row_idx_map[i] = i-offset;
}
}
}
// csc erasing rows and columns
int read(0), write(0), evar_col(0), last_jc(0);
for( unsigned int c = 0; c < nc; ++c)
{
if( c == (unsigned int)evar[evar_col] )
{
++evar_col;
read += _A.jc[c+1]-last_jc;
}
else
{
while( read < (int)_A.jc[c+1] )
{
int new_idx = new_row_idx_map[_A.ir[read]];
if( new_idx != -1)
{
_A.pr[write] = _A.pr[read];
_A.ir[write] = new_idx;
++write;
}
++read;
}
}
last_jc = _A.jc[c+1];
_A.jc[c+1-evar_col] = write;
}
_A.nc = nc - evar.size()+1;
_A.nr = nr - evar.size()+1;
}
//-----------------------------------------------------------------------------
template<class MatrixT, class REALT, class INTT>
void get_ccs_symmetric_data( const MatrixT& _mat,
const char _uplo,
std::vector<REALT>& _values,
std::vector<INTT>& _rowind,
std::vector<INTT>& _colptr )
{
unsigned int m = gmm::mat_nrows( _mat );
unsigned int n = gmm::mat_ncols( _mat );
gmm::csc_matrix<REALT> csc_mat( m,n );
gmm::copy( _mat, csc_mat );
_values.resize( 0 );
_rowind.resize( 0 );
_colptr.resize( 0 );
INTT iv( 0 );
switch ( _uplo )
{
case 'l':
case 'L':
// for all columns
for ( unsigned int i=0; i<n; ++i )
{
_colptr.push_back( iv );
for ( unsigned int j=( csc_mat.jc )[i]; j<( csc_mat.jc )[i+1]; ++j )
if (( csc_mat.ir )[j] >= i )
{
++iv;
_values.push_back(( csc_mat.pr )[j] );
_rowind.push_back(( csc_mat.ir )[j] );
}
}
_colptr.push_back( iv );
break;
case 'u':
case 'U':
// for all columns
for ( unsigned int i=0; i<n; ++i )
{
_colptr.push_back( iv );
for ( unsigned int j=( csc_mat.jc )[i]; j<( csc_mat.jc )[i+1]; ++j )
if (( csc_mat.ir )[j] <= i )
{
++iv;
_values.push_back(( csc_mat.pr )[j] );
_rowind.push_back(( csc_mat.ir )[j] );
}
}
_colptr.push_back( iv );
break;
case 'c':
case 'C':
// for all columns
for ( unsigned int i=0; i<n; ++i )
{
_colptr.push_back( iv );
for ( unsigned int j=( csc_mat.jc )[i]; j<( csc_mat.jc )[i+1]; ++j )
// if (( csc_mat.ir )[j] <= i )
{
++iv;
_values.push_back(( csc_mat.pr )[j] );
_rowind.push_back(( csc_mat.ir )[j] );
}
}
_colptr.push_back( iv );
break;
default:
std::cerr << "ERROR: parameter uplo must bei either 'U' or 'L' or 'C'!!!\n";
break;
}
}
//-----------------------------------------------------------------------------
template<class ScalarT, class VectorT, class SVT>
void eliminate_var( const unsigned int _j,
const ScalarT _value_j,
gmm::col_matrix<SVT>& _A,
VectorT& _x,
VectorT& _rhs )
{
//ACG::StopWatch sw1;
//sw1.start();
typedef typename gmm::linalg_traits< gmm::col_matrix< SVT > >::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
unsigned int m = gmm::mat_nrows( _A );
unsigned int n = gmm::mat_ncols( _A );
// update rhs
ColT col = mat_const_col( _A, _j );
CIter it = gmm::vect_const_begin( col );
CIter ite = gmm::vect_const_end( col );
for ( ; it!=ite; ++it )
{
_rhs[it.index()] -= _value_j*( *it );
}
_rhs.erase( _rhs.begin() + _j );
_x.erase( _x.begin() + _j );
//std::cerr << " FIRST BLOCK " << sw1.stop()/1000.0 << std::endl;
/*
for(unsigned int i=_j; i<m-1; ++i)
{
_rhs[i] = _rhs[i+1];
_x [i] = _x [i+1];
}
// delete last element
_rhs.resize( _rhs.size()-1);
_x.resize( _x.size()-1);
*/
/*
MatrixT A_temp( m, n);
gmm::copy(_A, A_temp);
gmm::clear(_A);
for(unsigned int i=0; i<m; ++i)
{
// skip _j'th column
if( i==_j) continue;
ColT col = mat_const_col(A_temp, i);
CIter it = gmm::vect_const_begin(col);
CIter ite = gmm::vect_const_end(col);
// compute new index
unsigned int i_new = i;
if( i>_j) --i_new;
for(; it!=ite; ++it)
{
unsigned int j_new = it.index();
if( j_new == _j) continue;
if( j_new > _j) --j_new;
_A(j_new,i_new) = *it;
// if( it.index() != j_new || i_new != i)
// _A(it.index(), i) = 0.0;
}
}
gmm::resize( _A, m-1, n-1);
*/
//sw1.start();
typedef typename gmm::linalg_traits<SVT>::const_iterator SIter;
for ( unsigned int i=0; i<m; ++i )
{
// skip _j'th column
if ( i==_j ) continue;
// compute new index
unsigned int i_new = i;
if ( i>_j ) --i_new;
SVT* cur_col = &( _A.col( i ) );
if ( i == i_new )
{
cur_col = new SVT( _A.col( i ) );
//gmm::copy( _A.col(i), *cur_col);
}
SIter it = gmm::vect_const_begin( *cur_col );
SIter ite = gmm::vect_const_end( *cur_col );
//mat_col(_A, i_new).clear();
SVT& new_col = _A.col( i_new );
new_col.clear();
for ( ; it!=ite; ++it )
{
unsigned int j_new = it.index();
if ( j_new == _j ) continue;
if ( j_new > _j ) --j_new;
//_A(j_new,i_new) = *it;
new_col[j_new] = *it;
// if( it.index() != j_new || i_new != i)
// _A(it.index(), i) = 0.0;
}
if ( i == i_new )
{
delete cur_col;
}
}
//std::cerr << " BLOCK TWO " << sw1.stop()/1000.0 << std::endl;
//sw1.start();
gmm::resize( _A, m-1, n-1 );
//std::cerr << "RESIZE " << sw1.stop()/1000.0 << std::endl;
}
//-----------------------------------------------------------------------------
template<class ScalarT, class VectorT, class RealT>
void eliminate_var( const unsigned int _i,
const ScalarT _xi,
typename gmm::csc_matrix<RealT>& _A,
VectorT& _x,
VectorT& _rhs )
{
unsigned int n = _A.nc;
unsigned int iv_old( 0 );
unsigned int iv_new( 0 );
// for all columns
for ( unsigned int j=0; j<n; ++j )
{
// update x and rhs
if ( j > _i )
{
_rhs[j-1] = _rhs[j];
_x [j-1] = _x [j];
}
if ( j == _i )
{
// update rhs
for ( unsigned int i=_A.jc[j]; i<_A.jc[j+1]; ++i )
{
_rhs[_A.ir[iv_old]] -= _xi*_A.pr[iv_old];
++iv_old;
}
}
else
{
// store index to update colptr
unsigned int iv_new_save = iv_new;
for ( unsigned int i=_A.jc[j]; i<_A.jc[j+1]; ++i )
{
if ( _A.ir[iv_old] < _i )
{
_A.ir[iv_new] = _A.ir[iv_old];
_A.pr[iv_new] = _A.pr[iv_old];
++iv_new;
}
else if ( _A.ir[iv_old] > _i )
{
_A.ir[iv_new] = _A.ir[iv_old]-1;
_A.pr[iv_new] = _A.pr[iv_old];
++iv_new;
}
++iv_old;
}
if ( j<_i )
_A.jc[j] = iv_new_save;
else
if ( j>_i )
_A.jc[j-1] = iv_new_save;
}
}
// store index to end
_A.jc[n-1] = iv_new;
// resize matrix and vectors
_A.nc = n-1;
_A.nr = n-1;
_x.resize( _x.size()-1 );
_rhs.resize( _rhs.size()-1 );
}
//-----------------------------------------------------------------------------
template<class IntegerT, class ScalarT, class VectorT, class MatrixT>
void eliminate_vars( const std::vector<IntegerT>& _evar,
const std::vector<ScalarT>& _eval,
MatrixT& _A,
VectorT& _x,
VectorT& _rhs )
{
std::cerr << __FUNCTION__ << std::endl;
typedef typename gmm::linalg_traits<MatrixT>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
// unsigned int m = gmm::mat_nrows( _A);
unsigned int n = gmm::mat_ncols( _A );
unsigned int n_new = n - _evar.size();
// modify rhs
for ( unsigned int k=0; k<_evar.size(); ++k )
{
IntegerT i = _evar[k];
ColT col = mat_const_col( _A, i );
CIter it = gmm::vect_const_begin( col );
CIter ite = gmm::vect_const_end( col );
for ( ; it!=ite; ++it )
{
_rhs[it.index()] -= _eval[k]*( *it );
}
}
// sort vector
std::vector<IntegerT> evar( _evar );
std::sort( evar.begin(), evar.end() );
evar.push_back( std::numeric_limits<int>::max() );
// build subindex set and update rhs
std::vector<size_t> si( n_new );
int cur_evar_idx=0;
for ( unsigned int i=0; i<n; ++i )
{