linear_algebra/least_squares/LSMRSolver.h

0001 /**
0002  * @file LSMRSolver.h
0003  * @author Nabil Chouika (Irfu/SPhN - CEA Saclay)
0004  * @date 2 déc. 2016
0005  * @version 1.0
0006  */
0007
0008 #ifndef LSMRSOLVER_H_
0009 #define LSMRSOLVER_H_
0010
0011 #include <string>
0012 #include <vector>
0013 #include <cstddef>
0014
0015 #include "../matrix/MatrixD.h"
0016 #include "../vector/VectorD.h"
0017
0018 namespace ElemUtils {
0019 class Formatter;
0020 } /* namespace ElemUtils */
0021
0022 namespace NumA {
0023
0024 /**
0025  * @class LSMRSolver
0026  * @brief LSMR solves Ax = b or min ||Ax - b|| with or without damping,
0027  * using the iterative algorithm of David Fong and Michael Saunders:
0028  *     http://www.stanford.edu/group/SOL/software/lsmr.html
0029  *
0030  * Adapted for C++ from the scipy implementation:
0031  * https://github.com/scipy/scipy/blob/v0.18.1/scipy/sparse/linalg/isolve/lsmr.py
0032  */
0033
0034 class LSMRSolver {
0035 public:
0036     /**
0037      * Constructor.
0038      * @param damp Damping factor for regularized least-squares.
0039      * @param atol Stopping tolerance.
0040      * @param btol Stopping tolerance.
0041      * @param conlim LSMR terminates if an estimate of `cond(A)` exceeds conlim.
0042      * @param maxiter LSMR terminates if the number of iterations reaches maxiter.
0043      * @param output Pointer to the stream where to write the output.
0044      */
0045     LSMRSolver(double damp = 0., double atol = 1.e-6, double btol = 1.e-6,
0046             double conlim = 1.e8, size_t maxiter = 0,
0047             ElemUtils::Formatter* output = 0);
0048     /**
0049      * Default destructor.
0050      */
0051     virtual ~LSMRSolver();
0052
0053     /**
0054      * @brief Solves the linear system \f$ A X = B \f$ and returns \f$ X \f$ in the least-squares sense.
0055      * Iterative method for large sparse matrices.
0056      * @param A MatrixD.
0057      * @param B VectorD.
0058      * @return VectorD Solution X.
0059      */
0060     VectorD solve(const MatrixD & A, const VectorD & B);
0061
0062     // ##### GETTERS & SETTERS #####
0063
0064     double getToleranceA() const; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
0065     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.
0066     double getToleranceB() const; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
0067     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.
0068     double getConditionNumberLimit() const; ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
0069     void setConditionNumberLimit(double conlim); ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
0070     double getDamping() const; ///< Damping factor for regularized least-squares.
0071     void setDamping(double damp); ///< Damping factor for regularized least-squares.
0072     size_t getMaximumIerations() const; ///< LSMR terminates if the number of iterations reaches maxiter.
0073     void setMaximumIerations(size_t maxiter); ///< LSMR terminates if the number of iterations reaches maxiter.
0074      /**
0075       *
0076       * @return Reference to the stream where the output is written.
0077       */
0078     ElemUtils::Formatter& getOutput() const;
0079     /**
0080      *
0081      * @param output Reference to the stream where the output is written.
0082      */
0083     void setOutput(ElemUtils::Formatter& output);
0084     double getConditionNumber() const;
0085     unsigned int getStoppingCase() const; ///< Reason for stopping the iterations.
0086     size_t getNumberIterations() const;  ///< Number of iterations.
0087     double getNormA() const;
0088     double getNormAr() const;
0089     double getNormr() const;
0090     double getNormX() const;
0091     const std::string& getStoppingMessage() const; ///< Reason for stopping the iterations.
0092
0093 private:
0094     // ##### Parameters of the solver #####
0095     double m_damp; ///< Damping factor for regularized least-squares.
0096     double m_atol; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
0097     double m_btol; ///< Stopping tolerance. LSMR continues iterations until a certain backward error estimate is smaller than some quantity depending on atol and btol.
0098     double m_conlim; ///< LSMR terminates if an estimate of `cond(A)` exceeds conlim.
0099     size_t m_maxiter; ///< LSMR terminates if the number of iterations reaches maxiter.
0100     ElemUtils::Formatter* m_output; ///< Pointer to the stream where to write the output.
0101
0102     // ##### Results of the solver #####
0103     unsigned int m_istop; ///< Reason for stopping the iterations.
0104     size_t m_itn; ///< Number of iterations.
0105     double m_normr, m_normar, m_normA, m_condA, m_normx;
0106     std::string m_stopMessage; ///< Reason for stopping the iterations.
0107
0108     // Miscellaneous
0109     std::vector<std::string> m_msg;
0110     std::string m_hdg1;
0111     std::string m_hdg2;
0112     unsigned int m_pfreq;   // print frequency (for repeating the heading)
0113
0114     /**
0115      * Routine of rotation from scipy.
0116      * @param a
0117      * @param b
0118      * @param c
0119      * @param s
0120      * @param r
0121      */
0122     void _sym_ortho(double a, double b, double& c, double& s, double& r) const;
0123 };
0124
0125 } /* namespace NumA */
0126
0127 #endif /* LSMRSOLVER_H_ */