[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [29271] branches/soc-2010-rohith291991/ intern/comiso/Examples/quadratic_solver: everything compiles.
Rohith B V
rohith291991 at gmail.com
Sun Jun 6 13:43:47 CEST 2010
Revision: 29271
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=29271
Author: rohith291991
Date: 2010-06-06 13:43:47 +0200 (Sun, 06 Jun 2010)
Log Message:
-----------
everything compiles. CG remains in the solver. Output is believable though correctness has to be checked.
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/common.cc
branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.h
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-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.cc 2010-06-06 11:43:47 UTC (rev 29271)
@@ -183,7 +183,7 @@
void
ConstrainedSolver::
restore_eliminated_vars( MatrixXd& _constraints,
- VectorXd& _x,
+ std::vector<double>& _x,
std::vector<int>& _c_elim,
std::vector<int>& _new_idx)
{
@@ -217,8 +217,15 @@
double cur_val = _constraints(i, cur_var);
_constraints(i, cur_var) = 0.0;
+double dot=0;
- _x[cur_var] = _x.dot(_constraints.row(i))/cur_val;
+for(int l=0;l<_x.size();l++)
+ {
+
+ dot+=_x[l]*_constraints.row(i)[l];
+ }
+
+ _x[cur_var] = dot/cur_val;
}
}
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-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/ConstrainedSolver.h 2010-06-06 11:43:47 UTC (rev 29271)
@@ -78,7 +78,7 @@
void
ConstrainedSolver::
restore_eliminated_vars( MatrixXd& _constraints,
- VectorXd& _x,
+ std::vector<double>& _x,
std::vector<int>& _c_elim,
std::vector<int>& _new_idx);
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-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.cc 2010-06-06 11:43:47 UTC (rev 29271)
@@ -26,6 +26,7 @@
#include "common.h"
#include "MISolver.h"
#include <float.h>
+#include <queue>
#define ROUND(x) ((x)<0?int((x)-0.5):int((x)+0.5))
@@ -58,25 +59,49 @@
SparseLDLT<SparseXdc> cholSolver;
VectorXd x;
- std::cout<<_b.size()<<"\n\n";
-std::cout<<_A.rows()<<"\n\n";
+
+
x.resize(_b.size());
- int i;
+ int i,j;
for(i=0;i<_b.size();i++)
{
x[i]=_b[i];
}
+int dimn=_A.rows();
+ int dimm=_A.cols();
- //cholSolver.compute(_A);
- /*if(!cholSolver.solveInPlace(x))
+ SparseXdc Acsc;
+Acsc.resize(dimn,dimm);
+
+Acsc.startFill(dimn*dimm);
+
+ for(i=0;i<dimn;i++)
+ {
+
+ for(j=0;j<dimm;j++)
+
+ {
+
+ Acsc.fill(j,i)=_A(j,i);
+
+ }
+
+ }
+
+Acsc.endFill();
+
+
+ cholSolver.compute(Acsc);
+
+ if(!cholSolver.solveInPlace(x))
{
std::cerr<<"Error occured during Cholesky Solve"<<"\n\n";
return false;
- }*/
+ }
_x.resize(_b.size());
@@ -90,6 +115,60 @@
+int MISolver::gauss_seidel_local( MatrixXd& _A,
+ std::vector<double>& _x,
+ std::vector<double>& _rhs,
+ std::vector<unsigned int> _idxs,
+ int _max_iter,
+ double _tolerance )
+{
+
+
+ double t2 = _tolerance*_tolerance;
+
+ std::vector<unsigned int> i_temp;
+ std::queue<unsigned int> q;
+
+ for ( unsigned int i=0; i<_idxs.size(); ++i )
+ q.push( _idxs[i] );
+
+ int it_count = 0;
+
+ while ( !q.empty() && it_count < _max_iter )
+ {
+ ++it_count;
+ unsigned int cur_i = q.front();
+ q.pop();
+ i_temp.clear();
+
+
+ double res_i = -_rhs[cur_i];
+ double x_i_new = _rhs[cur_i];
+
+
+ for (int l=0 ; l<_A.rows(); ++l )
+ if(_A(l,cur_i)!=0)
+ {
+ res_i += ( _A(l,cur_i) ) * _x[l];
+ x_i_new -= ( _A(l,cur_i) ) * _x[l];
+ i_temp.push_back( l );
+ }
+
+ if ( res_i*res_i > t2 )
+ {
+ _x[cur_i] += x_i_new/_A( cur_i, cur_i );
+
+ for ( unsigned int j=0; j<i_temp.size(); ++j )
+ q.push( i_temp[j] );
+ }
+ }
+
+
+ return it_count;
+}
+
+
+
void
MISolver::solve(
MatrixXd& _A,
@@ -142,7 +221,8 @@
// preconditioner for CG
//Check dimensions
- //MatrixXd PS, PR;
+ MatrixXd PS(_A.rows(),_A.cols());
+ MatrixXd PR(_A.rows(),_A.cols());
//PS.setIdentity();
//PR.setIdentity();
@@ -176,7 +256,7 @@
// eliminate vars from linear system
//TODO
- /* gmm::eliminate_csc_vars2( elim_i, elim_v, _A, xr, _rhs);*/
+ eliminate_csc_vars2( elim_i, elim_v, _A, xr, _rhs);
}
else
@@ -238,9 +318,9 @@
}
// eliminate var
- //----------TODO stuff--------
- eliminate_var(i_best, rnd_x, _A, xr, _rhs);
- eliminate_var_idx(i_best, to_round,-1);
+
+ eliminate_var(i_best, rnd_x, _A, xr, _rhs);
+ eliminate_var_idx(i_best, to_round, -1);
// update old_idx
old_idx.erase( old_idx.begin()+i_best );
@@ -257,8 +337,8 @@
if( noisy_ > 2)std::cerr << "use local iteration ";
- //------------TODO-------------
- // n_its = gmm::gauss_seidel_local(_A, xr, _rhs, neigh_i, max_local_iters_, max_local_error_);
+ n_its = gauss_seidel_local(_A, xr, _rhs, neigh_i, max_local_iters_, max_local_error_);
+
++n_local;
}
@@ -299,8 +379,7 @@
{
//TODO replace
- // chol.calc_system_gmm(_A);
- //chol.solve(xr,_rhs);
+ choleskySolve(_A, _rhs, xr);
++n_full;
}
@@ -309,7 +388,6 @@
if( noisy_ > 2)std::cerr << std::endl;
}
- std::cerr<<"TEST\n";
// final full solution?
if( final_full_solution_ || direct_rounding_)
@@ -320,8 +398,7 @@
{
//-----------TODO replace------------
- // chol.calc_system_gmm(_A);
- // chol.solve( xr, _rhs);
+ choleskySolve(_A, _rhs, xr);
++n_full;
}
}
@@ -330,6 +407,7 @@
for(unsigned int i=0; i<old_idx.size(); ++i)
{
_x[ old_idx[i] ] = xr[i];
+
}
// output statistics
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.h
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.h 2010-06-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/MISolver.h 2010-06-06 11:43:47 UTC (rev 29271)
@@ -89,6 +89,15 @@
Veci& _to_round,
bool _fixed_order = false );
+
+
+int gauss_seidel_local( MatrixXd& _A,
+ std::vector<double>& _x,
+ std::vector<double>& _rhs,
+ std::vector<unsigned int> _idxs,
+ int _max_iter,
+ double _tolerance );
+
bool choleskySolve(MatrixXd A, std::vector<double>& b, std::vector<double>& x);
/** @name Get/Set functions for algorithm parameters
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/comiso.cc
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/comiso.cc 2010-06-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/comiso.cc 2010-06-06 11:43:47 UTC (rev 29271)
@@ -144,7 +144,7 @@
int i;
int j;
int m = 9;
- int n = 10+1;
+ int n = 5+1;
std::cout << "---------- 1) setup an (m x n) sparse row matrix B (i.e. the B in the system ((Bx)^T)Bx)" << std::endl;
@@ -160,8 +160,14 @@
std::cout << A << " " << std::endl;
std::cout << std::endl;
+
+ for(int l=0;l<b.size();l++)
+ {
+ std::cout<<b[l]<<"\n";
+
+ }
std::cout << "---------- 3) define a set of linear constraints as an (c x n) row matrix C" << std::endl;
- int c = 3;
+ int c = 2;
MatrixXd C;
simple_constraint_row_matrix( C, c, n);// <- Replace with actual matrix
/*C(0,0)=2;
@@ -192,7 +198,8 @@
// vector of indices to round (this is the mixed-integer part)
std::vector< int > ids_to_round;
// lets say we want to round the third variable
- ids_to_round.push_back(2);
+ //ids_to_round.push_back(4);
+//ids_to_round.push_back(2);
// vector of independent variables to be eliminated (computed by the make_constraints_independent function)
std::vector< int > ids_to_elim;
@@ -262,18 +269,24 @@
std::cout << "---------- ---------- 4.4) now the solution must be re-indexed to the expected/original form/size...." << std::endl;
+
+
- VectorXd fx(x.size());
- for(i=0;i<x.size();i++)
+
+ cs.restore_eliminated_vars( Ccpy,x, ids_to_elim, new_idx);
+
+ std::cout << " fullsize solution vector x is\n" ;
+
+ for(int l=0;l<x.size();l++)
{
- fx[i]=x[i];
-
+ std::cout<<x[l]<<"\n";
+
}
- cs.restore_eliminated_vars( Ccpy,fx, ids_to_elim, new_idx);
+
+
- std::cout << " fullsize solution vector x is\n" << fx << std::endl << std::endl;
-
+ std::cout<<A<<"\n";
getchar();
return -1;
}
Modified: branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.cc
===================================================================
--- branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.cc 2010-06-06 09:40:25 UTC (rev 29270)
+++ branches/soc-2010-rohith291991/intern/comiso/Examples/quadratic_solver/common.cc 2010-06-06 11:43:47 UTC (rev 29271)
@@ -1,150 +1,144 @@
#include "common.h"
-
-double residuum_norm( MatrixXd& _A, std::vector<double>& _x, std::vector<double>& _rhs )
+void eliminate_var( const unsigned int _j,
+ const double _value_j,
+ MatrixXd& _A,
+ std::vector<double>& _x,
+ std::vector<double>& _rhs )
{
- if ( _A.cols() != _x.size() )
- std::cerr << "DIM ERROR (residuum_norm): " << _A.cols() << " " << _x.size() << std::endl;
- if ( _rhs.size() != _x.size() )
- std::cerr << "DIM ERROR 2 (residuum_norm): " << _rhs.size() << " " << _x.size() << std::endl;
- // temp vectors
- VectorXd Ax( _x.size() );
- VectorXd res( _x.size() );
- VectorXd x(_x.size());
- VectorXd rhs(_rhs.size());
+ unsigned int m = _A.rows();
+ unsigned int n = _A.cols();
- for(int i=0;i<_x.size();i++)
- {
- x[i]=_x[i];
- rhs[i]=_rhs[i];
- }
- Ax= _A*x;
- res= Ax- rhs;
+ // update rhs
- return res.norm();
-}
+ for (int i=0 ; i<m; ++i )
+ if(_A(i,_j)!=0)
+ {
+ _rhs[i] -= _value_j*( _A(i,_j) );
+ }
+ _rhs.erase( _rhs.begin() + _j );
+ _x.erase( _x.begin() + _j );
-void eliminate_csc_vars2(
- const std::vector<int>& _evar,
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list