Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Matrix.hpp"
       4             : 
       5             : /** 
       6             :  * @brief   LU factorization without pivoting.
       7             :  * 
       8             :  * Factorizes a square matrix into a lower triangular and an upper-triangular
       9             :  * factor.
      10             :  * 
      11             :  * This version does not use row pivoting, and is not rank-revealing.
      12             :  * 
      13             :  * @warning **Never** use this factorization, it is not numerically stable and
      14             :  *          will fail completely if a zero pivot is encountered. This algorithm
      15             :  *          is included for educational purposes only. Use a pivoted LU 
      16             :  *          factorization or a QR factorization instead.
      17             :  * 
      18             :  * @ingroup Factorizations
      19             :  */
      20             : class NoPivotLU {
      21             :   public:
      22             :     /// @name Constructors
      23             :     /// @{
      24             : 
      25             :     /// Default constructor.
      26             :     NoPivotLU() = default;
      27             :     /// Factorize the given matrix.
      28          11 :     NoPivotLU(const SquareMatrix &matrix) { compute(matrix); }
      29             :     /// Factorize the given matrix.
      30           1 :     NoPivotLU(SquareMatrix &&matrix) { compute(std::move(matrix)); }
      31             : 
      32             :     /// @}
      33             : 
      34             :   public:
      35             :     /// @name Factorization
      36             :     /// @{
      37             : 
      38             :     /// Perform the LU factorization of the given matrix.
      39             :     void compute(SquareMatrix &&matrix);
      40             :     /// Perform the LU factorization of the given matrix.
      41             :     void compute(const SquareMatrix &matrix);
      42             : 
      43             :     /// @}
      44             : 
      45             :   public:
      46             :     /// @name   Retrieving the L factor
      47             :     /// @{
      48             : 
      49             :     /// Get the lower-triangular matrix L, reusing the internal storage.
      50             :     /// @warning    After calling this function, the NoPivotLU object is no
      51             :     ///             longer valid, because this function steals its storage.
      52             :     SquareMatrix &&steal_L();
      53             : 
      54             :     /// Copy the lower-triangular matrix L to the given matrix.
      55             :     void get_L_inplace(Matrix &L) const;
      56             :     /// Get a copy of the lower-triangular matrix L.
      57             :     SquareMatrix get_L() const &;
      58             :     /// Get the lower-triangular matrix L.
      59           2 :     SquareMatrix &&get_L() && { return steal_L(); }
      60             : 
      61             :     /// @}
      62             : 
      63             :   public:
      64             :     /// @name   Retrieving the U factor
      65             :     /// @{
      66             : 
      67             :     /// Get the upper-triangular matrix U, reusing the internal storage.
      68             :     /// @warning    After calling this function, the NoPivotLU object is no
      69             :     ///             longer valid, because this function steals its storage.
      70             :     SquareMatrix &&steal_U();
      71             : 
      72             :     /// Copy the upper-triangular matrix U to the given matrix.
      73             :     void get_U_inplace(Matrix &U) const;
      74             :     /// Get a copy of the upper-triangular matrix U.
      75             :     SquareMatrix get_U() const &;
      76             :     /// Get the upper-triangular matrix U.
      77           2 :     SquareMatrix &&get_U() && { return steal_U(); }
      78             : 
      79             :     /// @}
      80             : 
      81             :   public:
      82             :     /// @name   Solving systems of equations problems
      83             :     /// @{
      84             : 
      85             :     /// Solve the system AX = B or LUX = B.
      86             :     /// Matrix B is overwritten with the result X.
      87             :     void solve_inplace(Matrix &B) const;
      88             :     /// Solve the system AX = B or LUX = B.
      89             :     Matrix solve(const Matrix &B) const;
      90             :     /// Solve the system AX = B or LUX = B.
      91             :     Matrix &&solve(Matrix &&B) const;
      92             :     /// Solve the system Ax = b or LUx = b.
      93             :     Vector solve(const Vector &B) const;
      94             :     /// Solve the system Ax = b or LUx = b.
      95             :     Vector &&solve(Vector &&B) const;
      96             : 
      97             :     /// @}
      98             : 
      99             :   public:
     100             :     /// @name   Access to internal representation
     101             :     /// @{
     102             : 
     103             :     /// Check if this object contains a factorization.
     104          13 :     bool is_factored() const { return state == Factored; }
     105             : 
     106             :     /// Check if this object contains valid L and U factors.
     107           8 :     bool has_LU() const { return is_factored(); }
     108             : 
     109             :     /// Get the internal storage of the upper-triangular matrix U and the strict
     110             :     /// lower-triangular part of matrix L.
     111             :     /// @warning    After calling this function, the LU object is no longer
     112             :     ///             valid, because this function steals its storage.
     113             :     SquareMatrix &&steal_LU();
     114             : 
     115             :     /// Get a copy of the internal storage of the upper-triangular matrix U and
     116             :     /// the strict lower-triangular part of matrix L.
     117           1 :     const SquareMatrix &get_LU() const & { return LU; }
     118             :     /// Get the internal storage of the upper-triangular matrix U and the strict
     119             :     /// lower-triangular part of matrix L.
     120           1 :     SquareMatrix &&get_LU() && { return steal_LU(); }
     121             : 
     122             :     /// @}
     123             : 
     124             :   private:
     125             :     /// The actual LU factorization algorithm.
     126             :     void compute_factorization();
     127             :     /// Back substitution algorithm for solving upper-triangular systems UX = B.
     128             :     void back_subs(const Matrix &B, Matrix &X) const;
     129             :     /// Forward substitution algorithm for solving lower-triangular systems
     130             :     /// LX = B.
     131             :     void forward_subs(const Matrix &B, Matrix &X) const;
     132             : 
     133             :   private:
     134             :     /// Result of a LU factorization: stores the upper-triangular
     135             :     /// matrix U and the strict lower-triangular part of matrix L. The diagonal
     136             :     /// elements of L are implicitly 1.
     137             :     SquareMatrix LU;
     138             : 
     139             :     enum State {
     140             :         NotFactored = 0,
     141             :         Factored = 1,
     142             :     } state = NotFactored;
     143             : };
     144             : 
     145             : /// Print the L and U matrices of an NoPivotLU object.
     146             : /// @related    NoPivotLU
     147             : std::ostream &operator<<(std::ostream &os, const NoPivotLU &lu);