[eigen] Iterative solver module

[Will S <billy.baloop@xxxxxxxxx>]

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
Description: Binary data