Commit 1f6c4c85 authored by Nabil Chouika's avatar Nabil Chouika
Browse files

Linear algebra documentation partons/core/partons#4.

parent 8615f98b
...@@ -18,12 +18,30 @@ namespace NumA { ...@@ -18,12 +18,30 @@ namespace NumA {
/** /**
* @class EigenUtils * @class EigenUtils
* @brief Tools for the Eigen wrapper. * @brief Tools for the Eigen wrapper.
*
* Workaround to not have PARTONS depend on Eigen.
* The methods defined here are only used in the cpp files of NumA++. Do not include this file in PARTONS.
*/ */
class EigenUtils { class EigenUtils {
public: public:
/**
* Convert a matrix from NumA type to Eigen type.
* @param A NumA matrix.
* @return Eigen matrix.
*/
static Eigen::MatrixXd convertToEigen(const MatrixD & A); static Eigen::MatrixXd convertToEigen(const MatrixD & A);
/**
* Convert a vector from NumA type to Eigen type.
* @param V NumA vector.
* @return Eigen vector.
*/
static Eigen::VectorXd convertToEigen(const VectorD & V); static Eigen::VectorXd convertToEigen(const VectorD & V);
/**
* Convert a vector from Eigen type to NumA type.
* @param V Eigen vector.
* @return NumA vector.
*/
static VectorD convertToVectorD(const Eigen::VectorXd & V); static VectorD convertToVectorD(const Eigen::VectorXd & V);
}; };
......
...@@ -26,61 +26,84 @@ namespace NumA { ...@@ -26,61 +26,84 @@ namespace NumA {
* @brief LSMR solves Ax = b or min ||Ax - b|| with or without damping, * @brief LSMR solves Ax = b or min ||Ax - b|| with or without damping,
* using the iterative algorithm of David Fong and Michael Saunders: * using the iterative algorithm of David Fong and Michael Saunders:
* http://www.stanford.edu/group/SOL/software/lsmr.html * http://www.stanford.edu/group/SOL/software/lsmr.html
*
* Adapted for C++ from the scipy implementation:
* https://github.com/scipy/scipy/blob/v0.18.1/scipy/sparse/linalg/isolve/lsmr.py
*/ */
class LSMRSolver { class LSMRSolver {
public: public:
/**
* Constructor.
* @param damp Damping factor for regularized least-squares.
* @param atol Stopping tolerance.
* @param btol Stopping tolerance.
* @param conlim LSMR terminates if an estimate of `cond(A)` exceeds conlim.
* @param maxiter LSMR terminates if the number of iterations reaches maxiter.
* @param output Pointer to the stream where to write the output.
*/
LSMRSolver(double damp = 0., double atol = 1.e-6, double btol = 1.e-6, LSMRSolver(double damp = 0., double atol = 1.e-6, double btol = 1.e-6,
double conlim = 1.e8, size_t maxiter = 0, double conlim = 1.e8, size_t maxiter = 0,
ElemUtils::Formatter* output = 0); ElemUtils::Formatter* output = 0);
/**
* Default destructor.
*/
virtual ~LSMRSolver(); virtual ~LSMRSolver();
/** /**
* @brief Solves the linear system \f$ A X = B \f$ and returns \f$ X \f$ in the least-squares sense. * @brief Solves the linear system \f$ A X = B \f$ and returns \f$ X \f$ in the least-squares sense.
* Iterative method for large sparse matrices. * Iterative method for large sparse matrices.
* @param A MatrixD * @param A MatrixD.
* @param B VectorD * @param B VectorD.
* @return VectorD Solution X. * @return VectorD Solution X.
*/ */
VectorD solve(const MatrixD & A, const VectorD & B); VectorD solve(const MatrixD & A, const VectorD & B);
// ##### GETTERS & SETTERS ##### // ##### GETTERS & SETTERS #####
double getToleranceA() const; double getToleranceA() const; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
void setToleranceA(double atol); void setToleranceA(double atol); ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
double getToleranceB() const; double getToleranceB() const; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
void setToleranceB(double btol); void setToleranceB(double btol); ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
double getConditionNumberLimit() const; double getConditionNumberLimit() const; ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
void setConditionNumberLimit(double conlim); void setConditionNumberLimit(double conlim); ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
double getDamping() const; double getDamping() const; ///< Damping factor for regularized least-squares.
void setDamping(double damp); void setDamping(double damp); ///< Damping factor for regularized least-squares.
size_t getMaximumIerations() const; size_t getMaximumIerations() const; ///< LSMR terminates if the number of iterations reaches maxiter.
void setMaximumIerations(size_t maxiter); void setMaximumIerations(size_t maxiter); ///< LSMR terminates if the number of iterations reaches maxiter.
/**
*
* @return Reference to the stream where the output is written.
*/
ElemUtils::Formatter& getOutput() const; ElemUtils::Formatter& getOutput() const;
/**
*
* @param output Reference to the stream where the output is written.
*/
void setOutput(ElemUtils::Formatter& output); void setOutput(ElemUtils::Formatter& output);
double getConditionNumber() const; double getConditionNumber() const;
unsigned int getStoppingCase() const; unsigned int getStoppingCase() const; ///< Reason for stopping the iterations.
size_t getNumberIterations() const; size_t getNumberIterations() const; ///< Number of iterations.
double getNormA() const; double getNormA() const;
double getNormAr() const; double getNormAr() const;
double getNormr() const; double getNormr() const;
double getNormX() const; double getNormX() const;
const std::string& getStoppingMessage() const; const std::string& getStoppingMessage() const; ///< Reason for stopping the iterations.
private: private:
// ##### Parameters of the solver ##### // ##### Parameters of the solver #####
double m_damp; double m_damp; ///< Damping factor for regularized least-squares.
double m_atol; double m_atol; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
double m_btol; double m_btol; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
double m_conlim; double m_conlim; ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
size_t m_maxiter; size_t m_maxiter; ///< LSMR terminates if the number of iterations reaches maxiter.
ElemUtils::Formatter* m_output; ElemUtils::Formatter* m_output; ///< Pointer to the stream where to write the output.
// ##### Results of the solver ##### // ##### Results of the solver #####
unsigned int m_istop; unsigned int m_istop; ///< Reason for stopping the iterations.
size_t m_itn; size_t m_itn; ///< Number of iterations.
double m_normr, m_normar, m_normA, m_condA, m_normx; double m_normr, m_normar, m_normA, m_condA, m_normx;
std::string m_stopMessage; std::string m_stopMessage; ///< Reason for stopping the iterations.
// Miscellaneous // Miscellaneous
std::vector<std::string> m_msg; std::vector<std::string> m_msg;
......
...@@ -22,6 +22,28 @@ class VectorD; ...@@ -22,6 +22,28 @@ class VectorD;
* @class MatrixD * @class MatrixD
* *
* @brief Represents a two-dimensional array of double. * @brief Represents a two-dimensional array of double.
*
* The indices start at 0 up to `size`-1.
*
* Examples:
* ```cpp
* NumA::MatrixD A(5,5); // 5x5 matrix with coefficients equal to 0.
* A.at(1,3) = 10.; // Defines one coefficient (second row, fourth column).
* NumA::MatrixD B = NumA::MatrixD::Id(5); // Identity matrix.
* B.at(4,4) = 5.; // Last row, last column.
* // Matrix operations:
* NumA::MatrixD C;
* C = A*B + B/2.;
* C.toString();
* ```
* This returns:
* ```
* 0.5 0 0 0 0
* 0 0.5 0 10 0
* 0 0 0.5 0 0
* 0 0 0 0.5 0
* 0 0 0 0 2.5
* ```
*/ */
class MatrixD { class MatrixD {
...@@ -33,7 +55,7 @@ public: ...@@ -33,7 +55,7 @@ public:
MatrixD(); MatrixD();
/** /**
* Create matrix with specific dimension and fill it with given values separate by comma. * Creates a matrix with specific dimensions and fills it with given values separated by comma.
* *
* TODO: Write something better? * TODO: Write something better?
* *
...@@ -48,7 +70,7 @@ public: ...@@ -48,7 +70,7 @@ public:
double first_value, ...); double first_value, ...);
/** /**
* Create matrix with specific dimension and fill it with 0. * Creates matrix with specific dimensions and fills it with 0.
* *
* @param _rowsNumber : number of rows * @param _rowsNumber : number of rows
* @param _columnsNumber : number of columns * @param _columnsNumber : number of columns
...@@ -58,14 +80,14 @@ public: ...@@ -58,14 +80,14 @@ public:
/** /**
* Copy constructor. * Copy constructor.
* *
* @param rhs * @param rhs Matrix to be copied.
*/ */
MatrixD(const MatrixD& rhs); MatrixD(const MatrixD& rhs);
/** /**
* VectorD constructor (creates a one-column matrix). * VectorD constructor (creates a one-column matrix).
* *
* @param rhs VectorD * @param rhs Vector to be converted to a matrix.
*/ */
MatrixD(const VectorD& rhs); MatrixD(const VectorD& rhs);
...@@ -75,62 +97,129 @@ public: ...@@ -75,62 +97,129 @@ public:
virtual ~MatrixD(); virtual ~MatrixD();
/** /**
* Assign a new matrix with all coefficients set to value. * Assigns a new matrix with all coefficients set to `value`.
* @param _rowsNumber * @param _rowsNumber Number of rows.
* @param _columnsNumber * @param _columnsNumber Number of columns.
* @param value * @param value Default value for all the coefficients.
*/ */
void assign(const size_t _rowsNumber, const size_t _columnsNumber, void assign(const size_t _rowsNumber, const size_t _columnsNumber,
double value = 0.); double value = 0.);
/**
* Sets an existing line of the matrix to new values given by a vector.
* @param i Index of the line to be modified.
* @param line Vector to copy in the matrix.
*/
void setLine(size_t i, const NumA::VectorD& line); void setLine(size_t i, const NumA::VectorD& line);
/**
* Adds a new line in the matrix (its dimension is therefore incremented) with the values given by a vector.
* @param i Index the row where the line will be inserted.
* @param line Vector to copy in the matrix.
*/
void addLine(size_t i, const NumA::VectorD& line); void addLine(size_t i, const NumA::VectorD& line);
/**
* Adds a new line at the end of the matrix (its dimension is therefore incremented) with the values given by a vector.
* @param line Vector to copy in the matrix.
*/
void appendLine(const NumA::VectorD& line); void appendLine(const NumA::VectorD& line);
/**
* Extracts a single line from the matrix.
* @param lineIndex Index of the row to be extracted.
* @return Vector corresponding to the row.
*/
NumA::VectorD getLine(const size_t lineIndex) const; NumA::VectorD getLine(const size_t lineIndex) const;
// Matrix/vector operations // Matrix/vector operations
/** /**
* Matrix vector multiplication. * Matrix vector multiplication.
* @param rhs * @param rhs Vector.
* @return Matrix . rhs * @return Vector \f$ M \cdot V \f$.
*/ */
NumA::VectorD operator*(const NumA::VectorD& rhs) const; NumA::VectorD operator*(const NumA::VectorD& rhs) const;
//std::vector<double> operator*=(const std::vector<double>& rhs); //std::vector<double> operator*=(const std::vector<double>& rhs);
// Matrix/matrix operations // Matrix/matrix operations
/**
* Matrix-matrix multiplication.
* @param rhs Matrix.
* @return Product of matrices.
*/
NumA::MatrixD operator*(const NumA::MatrixD& rhs) const; NumA::MatrixD operator*(const NumA::MatrixD& rhs) const;
/**
* Matrix addition.
* @param rhs Matrix.
* @return Sum of matrices.
*/
NumA::MatrixD operator+(const NumA::MatrixD& rhs) const; NumA::MatrixD operator+(const NumA::MatrixD& rhs) const;
/**
* Matrix subtraction.
* @param rhs Matrix.
* @return Difference of matrices.
*/
NumA::MatrixD operator-(const NumA::MatrixD& rhs) const; NumA::MatrixD operator-(const NumA::MatrixD& rhs) const;
// Matrix/double operations // Matrix/double operations
/**
* Multiplication of all the coefficients by a scalar.
* @param rhs Scalar.
* @return Matrix.
*/
NumA::MatrixD operator*(double rhs) const; NumA::MatrixD operator*(double rhs) const;
/**
* Addition coefficient by coefficient with a scalar.
* @param rhs Scalar.
* @return Matrix.
*/
NumA::MatrixD operator+(double rhs) const; NumA::MatrixD operator+(double rhs) const;
/**
* Subtraction coefficient by coefficient with a scalar.
* @param rhs Scalar.
* @return Matrix.
*/
NumA::MatrixD operator-(double rhs) const; NumA::MatrixD operator-(double rhs) const;
/**
* Division of all the coefficients by a scalar.
* @param rhs Scalar.
* @return Matrix.
*/
NumA::MatrixD operator/(double rhs) const; NumA::MatrixD operator/(double rhs) const;
/** /**
* Matrix transpose. * Matrix transpose.
* @return MatrixD transposed * @return Matrix transposed.
*/ */
NumA::MatrixD transpose() const; NumA::MatrixD transpose() const;
/** /**
* Matrix identity. * Matrix identity.
* @param n Dimension. * @param n Dimension.
* @return MatrixD id * @return Matrix identity.
*/ */
static NumA::MatrixD Id(unsigned int n); static NumA::MatrixD Id(unsigned int n);
/** /**
* Update the value at the given coordinates. * Update the value at the given indices.
* *
* @param i : row coordinate * @param i Row index.
* @param j : column coordinate * @param j Column index.
* @param value * @param value New value for the coefficient.
*/ */
void update(const size_t i, const size_t j, const double value); void update(const size_t i, const size_t j, const double value);
/**
* Element access.
* Can be used to set the value of a coefficient.
* @param i Row index.
* @param j Column index.
* @return Reference to a certain coefficient of the matrix.
*/
double & at(const size_t i, const size_t j); double & at(const size_t i, const size_t j);
/**
* Element access.
* @param i Row index.
* @param j Column index.
* @return Constant reference to a certain coefficient of the matrix.
*/
const double & at(const size_t i, const size_t j) const; const double & at(const size_t i, const size_t j) const;
/** /**
...@@ -143,19 +232,19 @@ public: ...@@ -143,19 +232,19 @@ public:
// ##### GETTERS & SETTERS ##### // ##### GETTERS & SETTERS #####
/** /**
* *
* @return Number of columns * @return Number of columns.
*/ */
size_t cols() const; size_t cols() const;
/** /**
* *
* @return Number of rows * @return Number of rows.
*/ */
size_t rows() const; size_t rows() const;
private: private:
std::vector<double> m_matrix; ///< a flat std::vector to represent the matrix std::vector<double> m_matrix; ///< A flat std::vector to represent the matrix.
size_t m_rowsNumber; ///< number of rows size_t m_rowsNumber; ///< Number of rows.
size_t m_columnsNumber; ///< number of columns size_t m_columnsNumber; ///< Number of columns.
}; };
} /* namespace NumA */ } /* namespace NumA */
......
...@@ -14,29 +14,76 @@ namespace NumA { ...@@ -14,29 +14,76 @@ namespace NumA {
/** /**
* @class Vector2D * @class Vector2D
*
* Object representing a two-dimensional vector.
*/ */
class Vector2D { class Vector2D {
public: public:
/**
* Default constructor.
*/
Vector2D(); Vector2D();
/**
* Constructor.
* @param x x-coordinate.
* @param y y-coordinate.
*/
Vector2D(double x, double y); Vector2D(double x, double y);
/**
* Default destructor.
*/
~Vector2D(); ~Vector2D();
// ##### GETTERS & SETTERS ##### // ##### GETTERS & SETTERS #####
/**
*
* @return x-coordinate.
*/
double getX() const; double getX() const;
/**
*
* @param x x-coordinate.
*/
void setX(double x); void setX(double x);
/**
*
* @return y-coordinate.
*/
double getY() const; double getY() const;
/**
*
* @param y y-coordinate.
*/
void setY(double y); void setY(double y);
/**
* Addition of two vectors.
* @param rhs Vector2D.
* @return Sum of two Vector2D.
*/
Vector2D operator+(Vector2D const &rhs); Vector2D operator+(Vector2D const &rhs);
/**
* Addition with a scalar.
* @param rhs Scalar.
*/
void operator+=(Vector2D const &rhs); void operator+=(Vector2D const &rhs);
/**
* Subtraction of two vectors
* @param rhs Vector2D.
* @return Subtraction of two Vector2D.
*/
Vector2D operator-(Vector2D const &rhs); Vector2D operator-(Vector2D const &rhs);
/**
* Subtraction with a scalar.
* @param rhs Scalar.
*/
void operator-=(Vector2D const &rhs); void operator-=(Vector2D const &rhs);
private: private:
double m_x; double m_x; ///< x-coordinate.
double m_y; double m_y; ///< y-coordinate.
}; };
} /* namespace NumA */ } /* namespace NumA */
......
...@@ -18,25 +18,57 @@ namespace NumA { ...@@ -18,25 +18,57 @@ namespace NumA {
/** /**
* @class Vector3D * @class Vector3D
*
* Object representing a three-dimensional vector.
*/ */
class Vector3D: public Vector2D { class Vector3D: public Vector2D {
public: public:
/**
* Default constructor.
*/
Vector3D(); Vector3D();
/**
* Constructor.
* @param x x-coordinate.
* @param y y-coordinate.
* @param z z-coordinate.
*/
Vector3D(double x, double y, double z); Vector3D(double x, double y, double z);
/**
* Default destructor.
*/
~Vector3D(); ~Vector3D();
// Scalar operation // Scalar operation
/**
* Scalar product.
* @param rhs Vector3D.