Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Matrix.hpp"
       4             : #include "PermutationMatrix.hpp"
       5             : 
       6             : /** 
       7             :  * @brief   LU factorization with row pivoting.
       8             :  * 
       9             :  * Factorizes a square matrix into a lower triangular and an upper-triangular
      10             :  * factor.
      11             :  * 
      12             :  * This version uses row pivoting, but it is not rank-revealing.
      13             :  * 
      14             :  * @ingroup Factorizations
      15             :  */
      16             : class RowPivotLU {
      17             :   public:
      18             :     /// @name Constructors
      19             :     /// @{
      20             : 
      21             :     /// Default constructor.
      22             :     RowPivotLU() = default;
      23             :     /// Factorize the given matrix.
      24          13 :     RowPivotLU(const SquareMatrix &matrix) { compute(matrix); }
      25             :     /// Factorize the given matrix.
      26           1 :     RowPivotLU(SquareMatrix &&matrix) { compute(std::move(matrix)); }
      27             : 
      28             :     /// @}
      29             : 
      30             :   public:
      31             :     /// @name Factorization
      32             :     /// @{
      33             : 
      34             :     /// Perform the LU factorization of the given matrix.
      35             :     void compute(SquareMatrix &&matrix);
      36             :     /// Perform the LU factorization of the given matrix.
      37             :     void compute(const SquareMatrix &matrix);
      38             : 
      39             :     /// @}
      40             : 
      41             :   public:
      42             :     /// @name   Retrieving the L factor
      43             :     /// @{
      44             : 
      45             :     /// Get the lower-triangular matrix L, reusing the internal storage.
      46             :     /// @warning    After calling this function, the LU object is no
      47             :     ///             longer valid, because this function steals its storage.
      48             :     ///             Stealing both L and P is allowed (if you do not steal U,
      49             :     ///             because it shares storage with L).
      50             :     SquareMatrix &&steal_L();
      51             : 
      52             :     /// Copy the lower-triangular matrix L to the given matrix.
      53             :     void get_L_inplace(Matrix &L) const;
      54             :     /// Get a copy of the lower-triangular matrix L.
      55             :     SquareMatrix get_L() const &;
      56             :     /// Get the lower-triangular matrix L.
      57           2 :     SquareMatrix &&get_L() && { return steal_L(); }
      58             : 
      59             :     /// @}
      60             : 
      61             :   public:
      62             :     /// @name   Retrieving the U factor
      63             :     /// @{
      64             : 
      65             :     /// Get the upper-triangular matrix U, reusing the internal storage.
      66             :     /// @warning    After calling this function, the LU object is no
      67             :     ///             longer valid, because this function steals its storage.
      68             :     ///             Stealing both U and P is allowed (if you do not steal L,
      69             :     ///             because it shares storage with U).
      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   Retrieving the P factor
      83             :     /// @{
      84             : 
      85             :     /// Get the permutation matrix P, reusing the internal storage.
      86             :     /// @warning    After calling this function, the LU object is no
      87             :     ///             longer valid, because this function steals its storage.
      88             :     ///             Stealing P and either L or U (not both) is allowed.
      89             :     PermutationMatrix &&steal_P();
      90             : 
      91             :     /// Get a copy of the permutation matrix P.
      92           2 :     PermutationMatrix get_P() const & { return P; }
      93             :     /// Get the permutation matrix P.
      94           1 :     PermutationMatrix &&get_P() && { return steal_P(); }
      95             : 
      96             :     /// @}
      97             : 
      98             :   public:
      99             :     /// @name   Solving systems of equations problems
     100             :     /// @{
     101             : 
     102             :     /// Solve the system AX = B or LUX = B.
     103             :     /// Matrix B is overwritten with the result X.
     104             :     void solve_inplace(Matrix &B) const;
     105             :     /// Solve the system AX = B or LUX = B.
     106             :     Matrix solve(const Matrix &B) const;
     107             :     /// Solve the system AX = B or LUX = B.
     108             :     Matrix &&solve(Matrix &&B) const;
     109             :     /// Solve the system Ax = b or LUx = b.
     110             :     Vector solve(const Vector &B) const;
     111             :     /// Solve the system Ax = b or LUx = b.
     112             :     Vector &&solve(Vector &&B) const;
     113             : 
     114             :     /// @}
     115             : 
     116             :   public:
     117             :     /// @name   Access to internal representation
     118             :     /// @{
     119             : 
     120             :     /// Check if this object contains a factorization.
     121           5 :     bool is_factored() const { return state == Factored; }
     122             : 
     123             :     /// Check if this object contains valid L and U factors.
     124          10 :     bool has_LU() const { return valid_LU; }
     125             : 
     126             :     /// Check if this object contains a valid permutation matrix P.
     127             :     bool has_P() const { return valid_P; }
     128             : 
     129             :     /// Get the internal storage of the upper-triangular matrix U and the strict
     130             :     /// lower-triangular part of matrix L.
     131             :     /// @warning    After calling this function, the LU object is no longer
     132             :     ///             valid, because this function steals its storage.
     133             :     ///             Stealing both LU and P is allowed (but not L or U
     134             :     ///             individually).
     135             :     SquareMatrix &&steal_LU();
     136             : 
     137             :     /// Get a copy of the internal storage of the upper-triangular matrix U and
     138             :     /// the strict lower-triangular part of matrix L.
     139           1 :     const SquareMatrix &get_LU() const & { return LU; }
     140             :     /// Get the internal storage of the upper-triangular matrix U and the strict
     141             :     /// lower-triangular part of matrix L.
     142           1 :     SquareMatrix &&get_LU() && { return steal_LU(); }
     143             : 
     144             :     /// @}
     145             : 
     146             :   private:
     147             :     /// The actual LU factorization algorithm.
     148             :     void compute_factorization();
     149             :     /// Back substitution algorithm for solving upper-triangular systems UX = B.
     150             :     void back_subs(const Matrix &B, Matrix &X) const;
     151             :     /// Forward substitution algorithm for solving lower-triangular systems
     152             :     /// LX = B.
     153             :     void forward_subs(const Matrix &B, Matrix &X) const;
     154             : 
     155             :   private:
     156             :     /// Result of a LU factorization: stores the upper-triangular
     157             :     /// matrix U and the strict lower-triangular part of matrix L. The diagonal
     158             :     /// elements of L are implicitly 1.
     159             :     SquareMatrix LU;
     160             :     /// The permutation of A that maximizes pivot size.
     161             :     PermutationMatrix P = PermutationMatrix::RowPermutation;
     162             : 
     163             :     enum State {
     164             :         NotFactored = 0,
     165             :         Factored = 1,
     166             :     } state = NotFactored;
     167             : 
     168             :     bool valid_LU = false;
     169             :     bool valid_P = false;
     170             : };
     171             : 
     172             : /// Print the L and U matrices of an LU object.
     173             : /// @related    RowPivotLU
     174             : std::ostream &operator<<(std::ostream &os, const RowPivotLU &lu);