[eigen] Iterative solver module |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: [eigen] Iterative solver module*From*: Will S <billy.baloop@xxxxxxxxx>*Date*: Wed, 1 Sep 2010 10:57:45 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:date:message-id :subject:from:to:content-type; bh=6Uk3CWwa/Vyqvlzx3D1votGhERZF2j/iy6igwNFsCRI=; b=jxaRxKi9RC/NPBr1xlA2baIHhffbUHAka0VnHa2eFOco0/UiylDAGSBcPAsZ1RInym qYgbzO0qZv+vFLM4SRRslFyWU9pFgGrIAmdA1x1uZWrkP7TWkfY4cv2H3R8KdAlrE3H0 ytE1+iwDlwmb8uaQhRUlQMgj9TeSe5paSCjLc=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; b=IeK5gaY5YeZ2i4ox1ML9opWHN3e2hntIDJmvY9Eurg9bOVF+AE2J7bP5zRnc6Hvd/E ONRbB/RnF+ot/09V9VDGhIREZiXOUjuYdje6K+UxLNGyTedZi+eW1LfrbW+EkY7oBzXV Iw24y0gj8o5lfKOCbbDQWwTCSxBTv4RMeqfYY=

I needed an iterative solver for my project, so I implemented Jacobi, Gauss-Seidel, and Southwell iterative solvers. They work well, but they probably can be made much faster. Consider only the API. Mentioned in http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/06/msg00269.html My notes: - I tried to follow the Eigen model of providing meaningful names, and separate concrete types for the different solver objects. - I think passing matrix A as a pointer breaks some of the Eigen3 conventions, but is necessary to avoid copying large matrices (unless we accept const Eigen::MatrixXd & A, and retain a pointer to it that way). - I think the use of inheritance is good for forcing consistency between the solver classes but potentially bad for efficiency Comments? #ifndef ITERATIVE_SOLVER_H #define ITERATIVE_SOLVER_H #include <Eigen/Eigen> #define real double class IterativeSolver { public: enum SolutionState { /// Solver hasn't even /// started yet Indeterminate, /// The solution is in progress /// and has neither converged nor /// diverged yet, but it looks /// like it is converging as the /// error is decreasing IsConverging, /// The solution has completed /// and we have a good "x" solution Converged, /// The solution is going backwards /// and error is increasing with each /// iteration IsDiverging, /// Diverged. Got an infinity somewhere. /// The iterative solver /// selected cannot actually solve /// the system specified Diverged } ; protected: SolutionState solutionState ; /// The main system to solve. /// A copy of the system A is /// NOT made, because that's a /// huge waste of memory. /// The IterativeSolver /// series of classes will never /// modify the system A const Eigen::MatrixXd * A ; /// The solution vector that we /// want A*x = b. const Eigen::VectorXd * b ; /// A has N rows, the solution vector /// has N columns. int N ; /// 'curRow' is an important state variable /// its the row we are currently relaxing long long curRow ; /// 'iteration' is different from I. /// It is the MAJOR outer loop iteration. /// For Gauss-Seidel and Jacobi, one iteration /// constitutes a FULL PASS through the /// all the rows of the A matrix, /// but for the Southwell algorithm, its /// relaxation of a single row of the A /// matrix. At each "iteration", the /// Southwell needs to find the maximal error /// row, so, its concept of an "iteration" is /// different. long long iteration ; // long long.. support a huge number of iterations! /// how far each element may be from 0.0 /// in order for xGuess to be close enough /// to the true value of "x" for us to call /// it quits on the iterative process /// Solving to an iterative solver /// necessarily is accompanied by /// "to what precision"? real tolerance ; /// Every iterative technique has /// a residual that needs to be /// updated after each _major_ iteration Eigen::VectorXd residual ; /// xGuess -- our progressively /// (iteratively) refined estimate /// of what "x" should be Eigen::VectorXd xGuess ; /// Maintain copies of the length /// squared of the xGuess vector /// for this iteration AND the /// previous iteration. Difference /// between them can be used to /// identify whether solution is /// in fact converging or possibly diverging real norm2, lastNorm2 ; public: IterativeSolver( const Eigen::MatrixXd * A_System, const Eigen::VectorXd * b_solution, real iTolerance ) ; /// returns the real getLastNorm() ; /// Gives you the estimated /// solution vector x as computed /// by the iterative solver Eigen::VectorXd& getSolutionX() ; /// The "residual" is a vector that /// measures "how far off" our xG /// currently is from the true solution /// that should produce A*xG - b = 0 void updateResidual() ; /// Updates the solutionState member /// and returns to you the current /// SolutionState SolutionState checkConvergence() ; /// Returns true if a Eigen::MatrixXd is /// diagonally dominant. DD is /// a sufficient condition for /// convergence with the Gauss-Seidel method. bool isDiagonallyDominant( const Eigen::MatrixXd * mat ) ; /// Tells you if the matrix is DD, /// and returns the first row that /// was not DD bool isDiagonallyDominant( const Eigen::MatrixXd * mat, int *rowNotDD ) ; /// Completely solves the system to /// the tolerance precision you have /// specified in the ctor call virtual void solve() = 0 ; /// Things to be done prior to an /// iterate() call virtual void preIterate() = 0 ; /// Performs an actual iteration virtual void iterate() = 0 ; /// Things to be done following an /// iterate call() virtual void postIterate() = 0 ; } ; #endif

**Attachment:
IterativeSolver_interface_draft_1.h**

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Eigen/src/plugins/BlockMethods.h absent from 'make install' target** - Next by Date:
**Re: [eigen] eigen3 SelfAdjointEigenSolver<Matrix3f> still "buggy"** - Previous by thread:
**Re: [eigen] SSE floor/ceil** - Next by thread:
**[eigen] eigen3 compile error on jaunty (gcc 4.3.3)**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |