[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [29235] branches/soc-2010-rohith291991/ intern/comiso/Examples/quadratic_solver: More updates to the solver.
Rohith B V
rohith291991 at gmail.com
Sat Jun 5 13:12:58 CEST 2010
Revision: 29235
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=29235
Author: rohith291991
Date: 2010-06-05 13:12:56 +0200 (Sat, 05 Jun 2010)
Log Message:
-----------
More updates to the solver. Not finished yet.
Modified Paths:
--------------
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.cc
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.h
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.cc
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.h
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/comiso.cc
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/comiso.h
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.h
Added Paths:
-----------
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.cc
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.cc
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.cc 2010-06-05 10:15:48 UTC (rev 29234)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.cc 2010-06-05 11:12:56 UTC (rev 29235)
@@ -26,47 +26,7 @@
#include "ConstrainedSolver.h"
-template<class IntegerT, class IntegerT2>
-void ConstrainedSolver::eliminate_vars_idx( const std::vector<IntegerT >& _evar,
- std::vector<IntegerT2>& _idx,
- IntegerT2 _dummy )
-{
- // sort input
- std::vector<IntegerT> evar( _evar );
- std::sort( evar.begin(), evar.end() );
- // precompute update
- std::vector<int> update_map( _idx.size() );
-
- typename std::vector<IntegerT>::iterator cur_var = evar.begin();
- unsigned int deleted_between = 0;
-
- for ( unsigned int i=0; i<update_map.size(); ++i )
- {
- if ( cur_var != evar.end() && (IntegerT)i == *cur_var )
- {
- // mark as deleted
- update_map[i] = _dummy;
-
- ++deleted_between;
- ++cur_var;
- }
- else
- {
- update_map[i] = i-deleted_between;
- }
- }
-
- for ( unsigned int i=0; i<_idx.size(); ++i )
- {
- if ( (IntegerT2)_idx[i] != _dummy )
- {
- _idx[i] = update_map[_idx[i]];
- }
- }
-}
-
-
void
ConstrainedSolver::
add_row_simultaneously( int _row_i,
@@ -196,7 +156,7 @@
int nc=constraints_c.cols();
// copy col
-int k,l;
+int l;
for(l=i+1 ;l<nr;l++)
@@ -220,132 +180,6 @@
}
-void ConstrainedSolver::eliminate_csc_vars2(
- const std::vector<int>& _evar,
- const std::vector<double>& _eval,
- MatrixXd& _A,
- std::vector<double>& _x,
- std::vector<double>& _rhs )
-{
- typedef unsigned int uint;
- // typedef typename gmm::csc_matrix<RealT> MatrixT;
-
- unsigned int nc = _A.cols(); // check
- unsigned int nr = _A.rows();//
-
- unsigned int n_new = nc - _evar.size();
-
- // modify rhs
- for ( unsigned int k=0; k<_evar.size(); ++k )
- {
- int i = _evar[k];
-
- // number of elements in this column
-
- uint r_idx = 0;
- double entry = 0.0;
-
- for( uint i_elem = 0; i_elem < nr; ++i_elem)
- if(_A(i_elem,i)!=0)
- {
-
- entry = _A(i_elem,i);
- _rhs[i_elem] -= _eval[k]*( entry );
- }
- }
-
- // sort vector
- std::vector<int> 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<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;*/
-}
-
-
-
-
-
void
ConstrainedSolver::
restore_eliminated_vars( MatrixXd& _constraints,
@@ -511,15 +345,16 @@
// BIG TODO
eliminate_csc_vars2( elim_varids, elim_varvals, _A, _x, _rhs);
- int dimn=_A.rows();
+ //_A=MatrixXd(_A.block(1,1,7,7));
+
+ int dimn=_A.rows();
int dimm=_A.cols();
-int i,j;
+ int i,j;
_Acsc.resize(dimn,dimm);
_Acsc.startFill(dimn*dimm);
-
for(i=0;i<dimn;i++)
{
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.h
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.h 2010-06-05 10:15:48 UTC (rev 29234)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.h 2010-06-05 11:12:56 UTC (rev 29235)
@@ -39,8 +39,8 @@
#include "StopWatch.h"
//#include "MISolver.hh"
-
-#include "H:\gsocsvn\blender\intern\comiso\Examples\quadratic_solver\common.h"
+#include "common.h"
+//#include "H:\gsocsvn\blender\intern\comiso\Examples\quadratic_solver\common.h"
#include <vector>
//== FORWARDDECLARATIONS ======================================================
@@ -59,12 +59,6 @@
{
public:
-
-template<class IntegerT, class IntegerT2>
-void eliminate_vars_idx( const std::vector<IntegerT >& _evar,
- std::vector<IntegerT2>& _idx,
- IntegerT2 _dummy );
-
void
ConstrainedSolver::
@@ -73,15 +67,6 @@
std::vector<int>& _idx_to_round,
std::vector<int>& _c_elim);
-
-void ConstrainedSolver::eliminate_csc_vars2(
- const std::vector<int>& _evar,
- const std::vector<double>& _eval,
- MatrixXd& _A,
- std::vector<double>& _x,
- std::vector<double>& _rhs );
-
-
void
ConstrainedSolver::
add_row_simultaneously( int _row_i,
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.cc
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.cc 2010-06-05 10:15:48 UTC (rev 29234)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.cc 2010-06-05 11:12:56 UTC (rev 29235)
@@ -23,15 +23,13 @@
\*===========================================================================*/
-
+#include "common.h"
#include "MISolver.h"
#include <float.h>
#define ROUND(x) ((x)<0?int((x)-0.5):int((x)+0.5))
-
-
// Constructor
MISolver::MISolver()
{
@@ -54,18 +52,37 @@
//-----------------------------------------------------------------------------
-bool MISolver::choleskySolve(SparseXdc A, VectorXd b, VectorXd& x)
+bool MISolver::choleskySolve(MatrixXd _A, std::vector<double>& _b, std::vector<double>& _x)
{
SparseLDLT<SparseXdc> cholSolver;
- x=b;
+ VectorXd x;
- cholSolver.compute(A);
+ std::cout<<_b.size()<<"\n\n";
+std::cout<<_A.rows()<<"\n\n";
- if(!cholSolver.solveInPlace(x))
+ x.resize(_b.size());
+ int i;
+
+ for(i=0;i<_b.size();i++)
{
+ x[i]=_b[i];
+ }
+
+
+ //cholSolver.compute(_A);
+
+ /*if(!cholSolver.solveInPlace(x))
+ {
std::cerr<<"Error occured during Cholesky Solve"<<"\n\n";
return false;
+ }*/
+
+ _x.resize(_b.size());
+
+ for(i=0;i<_b.size();i++)
+ {
+ _x[i]=x[i];
}
return true;
@@ -74,7 +91,7 @@
void
MISolver::solve(
- SparseXdc& _A,
+ MatrixXd& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round,
@@ -105,24 +122,31 @@
for(unsigned int i=0; i<old_idx.size(); ++i)
old_idx[i] = i;
- /*
+
// Setup Cholmod solver used for full solution
- ACG::CholmodSolver chol;
if( initial_full_solution_ || direct_rounding_)
{
if( noisy_ > 2) std::cerr << "initial full solution" << std::endl;
- chol.calc_system_gmm(_A);
- chol.solve(_x, _rhs);
+ choleskySolve(_A, _rhs, _x);
++n_full;
}
+
+
if( no_rounding_)
return;
+
+ // preconditioner for CG
- // preconditioner for CG
- gmm::identity_matrix PS, PR;
+ //Check dimensions
+ //MatrixXd PS, PR;
+ //PS.setIdentity();
+ //PR.setIdentity();
+
+
+
// neighbors for local optimization
Vecui neigh_i;
@@ -151,10 +175,10 @@
// eliminate vars from linear system
//TODO
- gmm::eliminate_csc_vars2( elim_i, elim_v, _A, xr, _rhs);
+ /* gmm::eliminate_csc_vars2( elim_i, elim_v, _A, xr, _rhs);*/
}
- else
+ else
// loop until solution computed
for(unsigned int i=0; i<to_round.size(); ++i)
@@ -162,8 +186,9 @@
if( noisy_ > 0)
{
std::cerr << "Integer DOF's left: " << to_round.size()-(i+1) << " ";
- if( noisy_ > 1)
- std::cerr << "residuum_norm: " << gmm::residuum_norm( _A, xr, _rhs) << std::endl;
+ //if( noisy_ > 1)
+ //TODO
+// std::cerr << "residuum_norm: " << gmm::residuum_norm( _A, xr, _rhs)<<std::endl;
}
// index to eliminate
@@ -192,6 +217,7 @@
}
}
+
// store rounded value
double rnd_x = ROUND(xr[i_best]);
_x[ old_idx[i_best] ] = rnd_x;
@@ -200,22 +226,22 @@
neigh_i.clear();
//TODO
- Col col = gmm::mat_const_col(_A, i_best);
- ColIter it = gmm::vect_const_begin( col);
- ColIter ite = gmm::vect_const_end ( col);
+
- for(; it!=ite; ++it)
- if(it.index() < i_best)
- neigh_i.push_back(it.index());
+ for(int l=0; l<_A.rows(); ++l)
+ if(_A(l,i_best)!=0)
+ {
+ if(l < i_best)
+ neigh_i.push_back(l);
else
- if(it.index() > i_best)
- neigh_i.push_back(it.index()-1);
+ if(l > i_best)
+ neigh_i.push_back(l-1);
+ }
-
// eliminate var
//----------TODO stuff--------
- gmm::eliminate_var(i_best, rnd_x, _A, xr, _rhs);
- gmm::eliminate_var_idx(i_best, to_round);
+ //gmm::eliminate_var(i_best, rnd_x, _A, xr, _rhs);
+ //gmm::eliminate_var_idx(i_best, to_round);
// update old_idx
old_idx.erase( old_idx.begin()+i_best );
@@ -233,7 +259,7 @@
//------------TODO-------------
- n_its = gmm::gauss_seidel_local(_A, xr, _rhs, neigh_i, max_local_iters_, max_local_error_);
+ // n_its = gmm::gauss_seidel_local(_A, xr, _rhs, neigh_i, max_local_iters_, max_local_error_);
++n_local;
}
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list