LU decomposition
(src/algebra/linear/lu.hpp)

• View this file on GitHub
• Last update: 2021-04-15 16:09:50+09:00
• Include: #include "src/algebra/linear/lu.hpp"

Depends on

• Matrix (src/algebra/linear/matrix.hpp)

Verified with

• test/library-checker/matrix_det.test.cpp
• test/library-checker/system_of_linear_equations.test.cpp

Code

#pragma once

/**
 * @file lu.hpp
 * @brief LU decomposition
 */

#include <numeric>

#include "matrix.hpp"

namespace workspace {

template <class _Matrix> class lu_decomposition : public _Matrix {
 public:
  using value_type = typename _Matrix::value_type;
  using size_type = typename _Matrix::size_type;

  lu_decomposition() = default;
  lu_decomposition(const _Matrix &__x) : _Matrix(__x) { run(); }
  lu_decomposition(_Matrix &&__x) : _Matrix(std::move(__x)) { run(); }

 protected:
  std::vector<size_type> __perm, __pivots;
  bool sgn;

  void run() {
    __perm.resize(_Matrix::rows());
    std::iota(__perm.begin(), __perm.end(), 0);
    sgn = false;
    __pivots.clear();

    for (size_type __c = 0;
         __c != _Matrix::cols() && __pivots.size() != _Matrix::rows(); ++__c) {
      auto __max = (*this)[__pivots.size()][__c];
      auto __pos = __pivots.size();

      if constexpr (std::is_floating_point<
                        value_type>::value) {  // Find the biggest absolute
                                               // value in the column.
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows(); ++__r)
          if (std::abs(__max) < std::abs((*this)[__r][__c]))
            __max = (*this)[__pos = __r][__c];
      }

      else if (__max ==
               static_cast<value_type>(
                   0))  // Find the first non-zero element in the column.
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows(); ++__r)
          if ((__max = (*this)[__r][__c]) != static_cast<value_type>(0)) {
            __pos = __r;
            break;
          }

      if (__pos != __pivots.size()) {  // Swap 2 rows.
        sgn = !sgn;
        std::swap(__perm[__pos], __perm[__pivots.size()]);
        std::swap((*this)[__pos], (*this)[__pivots.size()]);
      }

      if (__max != static_cast<value_type>(0)) {  // Forward elimination
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows();
             ++__r) {
          auto __m = (*this)[__r][__c] / __max;
          (*this)[__r][__c] = 0;
          (*this)[__r][__pivots.size()] = __m;

          for (size_type __i = __c + 1; __i != _Matrix::cols(); ++__i)
            (*this)[__r][__i] -= (*this)[__pivots.size()][__i] * __m;
        }

        __pivots.emplace_back(__c);
      }
    }
  }

 public:
  size_type rank() const { return __pivots.size(); }

  value_type det() const {
    assert(_Matrix::rows() == _Matrix::cols());

    value_type __d = sgn ? -1 : 1;
    for (size_type __i = 0; __i != _Matrix::rows(); ++__i)
      __d *= (*this)[__i][__i];
    return __d;
  }

  _Matrix lower() const;

  _Matrix upper() const;

  _Matrix inverse() const {
    assert(_Matrix::rows() == _Matrix::cols());

    _Matrix __inv;
    // add solve(e_i) to __inv for i=0, ..., rows()
    return __inv;
  }

  // O(dim(ker) * size)
  _Matrix kernel() const {
    _Matrix __ker(_Matrix::cols() - rank(), _Matrix::cols());

    for (size_type __c = 0, __i = 0; __c != _Matrix::cols(); ++__c) {
      if (__i != __pivots.size() && __pivots[__i] == __c) {
        ++__i;
        continue;
      }

      auto &__v = __ker[__c - __i];
      __v[__c] = 1;
      for (size_type __r = 0; __r != __i; ++__r) __v[__r] = -(*this)[__r][__c];

      for (size_type __j = __i; __j--;) {
        auto __x = __v[__j] / (*this)[__j][__pivots[__j]];
        __v[__j] = 0;
        __v[__pivots[__j]] = __x;

        for (size_type __r = 0; __r != __j; ++__r)
          __v[__r] -= (*this)[__r][__pivots[__j]] * __x;
      }
    }

    return __ker;
  }

  template <class _Vec> std::pair<bool, _Vec> solve(const _Vec &__b) const {
    assert(!(__b.size() < _Matrix::rows()));

    // Solution
    _Vec __y(_Matrix::rows()), __x(_Matrix::cols());

    // Backward substitution with L
    for (size_type __c = 0; __c != _Matrix::rows(); ++__c) {
      __y[__c] += __b[__perm[__c]];

      for (size_type __r = __c + 1; __r != _Matrix::rows(); ++__r)
        __y[__r] -= __y[__c] * (*this)[__r][__c];
    }

    // Backward substitution with U
    for (size_type __i = rank(); __i != _Matrix::rows(); ++__i)
      if (__y[__i] != static_cast<value_type>(0))
        return std::make_pair(false, __x);

    for (size_type __i = rank(); __i--;) {
      auto __c = __pivots[__i];

      __x[__c] = __y[__i] / (*this)[__i][__c];

      for (size_type __r = 0; __r != __i; ++__r)
        __y[__r] -= __x[__c] * (*this)[__r][__c];
    }

    return std::make_pair(true, __x);
  }
};

}  // namespace workspace

#line 2 "src/algebra/linear/lu.hpp"

/**
 * @file lu.hpp
 * @brief LU decomposition
 */

#include <numeric>

#line 2 "src/algebra/linear/matrix.hpp"

/**
 * @file matrix.hpp
 * @brief Matrix
 * @date 2021-02-15
 *
 *
 */

#include <cassert>
#include <valarray>

namespace workspace {

/**
 * @brief Fixed size matrix.
 *
 * @tparam _Scalar
 * @tparam _Rows Number of rows
 * @tparam _Cols Number of columns
 */
template <class _Scalar, std::size_t _Rows = 0, std::size_t _Cols = _Rows>
class matrix {
 public:
  _Scalar __data[_Rows][_Cols] = {};

  using value_type = _Scalar;
  using size_type = std::size_t;

  constexpr static matrix eye() {
    static_assert(_Rows == _Cols);

    matrix __e;
    for (size_type __d = 0; __d != _Rows; ++__d) __e.__data[__d][__d] = 1;
    return __e;
  }

  constexpr operator decltype((__data))() { return __data; }
  constexpr operator decltype(std::declval<const matrix>().__data)
      const&() const {
    return __data;
  }

  constexpr auto begin() { return __data; }
  constexpr auto begin() const { return __data; }

  constexpr auto end() { return __data + _Rows; }
  constexpr auto end() const { return __data + _Rows; }

  constexpr size_type rows() const { return _Rows; }

  constexpr size_type cols() const { return _Cols; }

  constexpr auto transpose() const {
    matrix<_Scalar, _Cols, _Rows> __t;

    for (size_type __r = 0; __r != _Rows; ++__r)
      for (size_type __c = 0; __c != _Cols; ++__c)
        __t.__data[__c][__r] = __data[__r][__c];

    return __t;
  }

  constexpr matrix operator+() const { return *this; }

  constexpr matrix operator-() const {
    matrix __cp = *this;

    for (auto& __v : __cp.__data)
      for (auto& __e : __v) __e = -__e;

    return __cp;
  }

  template <class _Matrix> constexpr matrix& operator+=(const _Matrix& __x) {
    auto __m = std::min(_Rows, __x.rows());
    auto __n = std::min(_Cols, __x.cols());

    for (size_type __r = 0; __r != __m; ++__r)
      for (size_type __c = 0; __c != __n; ++__c)
        __data[__r][__c] += __x[__r][__c];

    return *this;
  }

  template <class _Matrix>
  constexpr matrix operator+(const _Matrix& __x) const {
    return matrix(*this) += __x;
  }

  template <class _Matrix> constexpr matrix& operator-=(const _Matrix& __x) {
    auto __m = std::min(_Rows, __x.rows());
    auto __n = std::min(_Cols, __x.cols());

    for (size_type __r = 0; __r != __m; ++__r)
      for (size_type __c = 0; __c != __n; ++__c)
        __data[__r][__c] -= __x[__r][__c];

    return *this;
  }

  template <class _Matrix>
  constexpr matrix operator-(const _Matrix& __x) const {
    return matrix(*this) -= __x;
  }

  template <class _Scalar2>
  constexpr matrix& operator*=(const matrix<_Scalar2, _Cols, _Cols>& __x) {
    if (this == &__x) return operator=(operator*(__x));

    for (auto& __r : __data) {
      _Scalar __tmp[_Cols] = {};

      auto __v = *__x.__data;
      for (auto& __w : __tmp) {
        auto __i = __v++;
        for (const auto& __e : __r) __w += __e * *__i, __i += _Cols;
      }

      auto __w = __tmp;
      for (auto& __e : __r) __e = std::move(*__w++);
    }

    return *this;
  }

  template <class _Scalar2, size_type _Rows2, size_type _Cols2>
  constexpr auto operator*(const matrix<_Scalar2, _Rows2, _Cols2>& __x) const {
    matrix<typename std::common_type<_Scalar, _Scalar2>::type, _Rows, _Cols2>
        __m;

    auto __w = *__m.__data;
    for (const auto& __r : __data)
      for (auto __v = *__x.__data, __v_end = __v + _Cols2; __v != __v_end;
           ++__w) {
        auto __i = __v++;
        for (auto __e = __r; __e != __r + std::min(_Cols, _Rows2); ++__e)
          *__w += *__e * *__i, __i += _Cols2;
      }

    return __m;
  }

  template <class _Matrix>
  constexpr
      typename std::enable_if<!std::is_convertible<_Matrix, value_type>::value,
                              matrix<_Scalar>>::type
      operator*(const _Matrix& __x) const {
    matrix<_Scalar> __m(_Rows, __x.cols());

    for (size_type __r = 0; __r != _Rows; ++__r)
      for (size_type __i = 0; __i != __x.cols(); ++__i)
        for (size_type __c = 0; __c != std::min(_Cols, __x.rows()); ++__c)
          __m[__r][__i] += __data[__r][__c] * __x[__c][__i];

    return __m;
  }

  constexpr matrix& operator*=(const value_type& __x) {
    for (auto& __v : __data)
      for (auto& __e : __v) __e *= __x;

    return *this;
  }

  constexpr matrix operator*(const value_type& __x) const {
    return matrix(*this) *= __x;
  }

  constexpr matrix& operator/=(const value_type& __x) {
    assert(__x != value_type(0));

    for (auto& __v : __data)
      for (auto& __e : __v) __e /= __x;

    return *this;
  }

  constexpr matrix operator/(const value_type& __x) const {
    return matrix(*this) /= __x;
  }

  template <class _Int> constexpr matrix pow(_Int __e) const {
    assert(0 <= __e);

    matrix __m = eye();
    for (matrix __cp = *this; __e; __cp *= __cp, __e >>= 1)
      if (__e & 1) __m *= __cp;

    return __m;
  }

  template <class _Os>
  constexpr friend _Os& operator<<(_Os& __os, const matrix& __x) {
    for (auto __i = __x.begin(); __i != __x.end(); ++__i, __os << '\n')
      for (size_type __c = 0; __c != _Cols; ++__c)
        __c ? void(__os << ' ') : (void)0, __os << *(*__i + __c);

    return __os;
  }
};  // namespace workspace

/**
 * @brief Dynamic matrix.
 *
 * @tparam _Scalar
 * @tparam _Rows Number of rows
 * @tparam _Cols Number of columns
 */
template <class _Scalar>
class matrix<_Scalar, 0, 0> : public std::valarray<std::valarray<_Scalar>> {
  using base = std::valarray<std::valarray<_Scalar>>;
  using row_type = typename base::value_type;

 public:
  using value_type = _Scalar;
  using size_type = std::size_t;

  using base::operator[];

  static matrix eye(size_type __n) {
    matrix __e(__n, __n);
    for (size_type __d = 0; __d != __n; ++__d) __e[__d][__d] = 1;
    return __e;
  }

  matrix() = default;

  matrix(size_type __n) : matrix(__n, __n) {}

  matrix(size_type __m, size_type __n) : base(row_type(__n), __m) {}

  template <class _Tp, typename = typename std::enable_if<
                           std::is_constructible<base, _Tp>::value &&
                           !std::is_constructible<size_type, _Tp>::value>::type>
  matrix(_Tp&& __x) : base(__x) {}

  matrix(std::initializer_list<row_type> __x) : base(__x) {}

  size_type rows() const { return base::size(); }

  size_type cols() const { return rows() ? operator[](0).size() : 0; }

  matrix transpose() const {
    matrix __t(cols(), rows());

    for (size_type __r = 0; __r != rows(); ++__r)
      for (size_type __c = 0; __c != cols(); ++__c)
        __t[__c][__r] = operator[](__r)[__c];

    return __t;
  }

  void resize(size_type __m, size_type __n) {
    matrix __t(__m, __n);

    if (rows() < __m) __m = rows();
    if (cols() < __n) __n = cols();

    for (size_type __r = 0; __r != __m; ++__r)
      for (size_type __c = 0; __c != __n; ++__c)
        __t[__r][__c] = std::move(operator[](__r)[__c]);

    base::swap(__t);
  }

  // binary operators {{

  template <class _Matrix, typename = void>
  struct is_valarray_based : std::false_type {};

  template <class _Matrix>
  struct is_valarray_based<
      _Matrix,
      typename std::enable_if<std::is_same<
          row_type, typename std::decay<decltype(
                        std::declval<_Matrix>()[0])>::type>::value>::type>
      : std::true_type {};

  template <class _Matrix>
  typename std::enable_if<!std::is_convertible<_Matrix, value_type>::value,
                          matrix&>::type
  operator*=(_Matrix&& __x) {
    return *this = operator*(std::forward<_Matrix>(__x));
  }

  template <class _Matrix>
  typename std::enable_if<!std::is_convertible<_Matrix, value_type>::value,
                          matrix>::type
  operator*(const _Matrix& __x) const {
    matrix __m(rows(), __x.cols());

    if constexpr (is_valarray_based<_Matrix>::value)
      for (size_type __r = 0; __r != rows(); ++__r)
        for (size_type __c = 0; __c != std::min(cols(), __x.rows()); ++__c)
          __m[__r] += operator[](__r)[__c] * __x[__c];

    else
      for (size_type __r = 0; __r != rows(); ++__r)
        for (size_type __i = 0; __i != __x.cols(); ++__i)
          for (size_type __c = 0; __c != std::min(cols(), __x.rows()); ++__c)
            __m[__r][__i] += operator[](__r)[__c] * __x[__c][__i];

    return __m;
  }

  matrix& operator*=(const value_type& __x) {
    for (size_type __r = 0; __r != rows(); ++__r)
      operator[](__r).operator*=(__x);

    return *this;
  }

  matrix operator*(const value_type& __x) const { return matrix(*this) *= __x; }

  friend matrix operator*(const value_type& __x, matrix __i) {
    for (size_type __r = 0; __r != __i.rows(); ++__r)
      __i.operator[](__r) = __x * __i.operator[](__r);

    return __i;
  }

  matrix& operator/=(const value_type& __x) {
    assert(__x != value_type(0));

    for (size_type __r = 0; __r != rows(); ++__r)
      operator[](__r).operator/=(__x);

    return *this;
  }

  matrix operator/(const value_type& __x) const { return matrix(*this) /= __x; }

  // }} binary operators

  template <class _Int> matrix pow(_Int __e) const {
    assert(0 <= __e);

    matrix __m = eye(rows());
    for (matrix __cp = *this; __e; __cp *= __cp, __e >>= 1)
      if (__e & 1) __m *= __cp;

    return __m;
  }

  // template <class _Is> friend _Is& operator>>(_Is& __is, matrix& __x) {
  //   for (size_type __r = 0; __r != __x.rows(); ++__r)
  //     for (size_type __c = 0; __c != __x.cols(); ++__c)
  //       __is >> __x.operator[](__r).operator[](__c);

  //   return __is;
  // }

  template <class _Os> friend _Os& operator<<(_Os& __os, const matrix& __x) {
    for (size_type __r = 0; __r != __x.rows(); ++__r, __os << '\n')
      for (size_type __c = 0; __c != __x.cols(); ++__c)
        __c ? void(__os << ' ') : (void)0,
            __os << __x.operator[](__r).operator[](__c);

    return __os;
  }
};

}  // namespace workspace
#line 11 "src/algebra/linear/lu.hpp"

namespace workspace {

template <class _Matrix> class lu_decomposition : public _Matrix {
 public:
  using value_type = typename _Matrix::value_type;
  using size_type = typename _Matrix::size_type;

  lu_decomposition() = default;
  lu_decomposition(const _Matrix &__x) : _Matrix(__x) { run(); }
  lu_decomposition(_Matrix &&__x) : _Matrix(std::move(__x)) { run(); }

 protected:
  std::vector<size_type> __perm, __pivots;
  bool sgn;

  void run() {
    __perm.resize(_Matrix::rows());
    std::iota(__perm.begin(), __perm.end(), 0);
    sgn = false;
    __pivots.clear();

    for (size_type __c = 0;
         __c != _Matrix::cols() && __pivots.size() != _Matrix::rows(); ++__c) {
      auto __max = (*this)[__pivots.size()][__c];
      auto __pos = __pivots.size();

      if constexpr (std::is_floating_point<
                        value_type>::value) {  // Find the biggest absolute
                                               // value in the column.
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows(); ++__r)
          if (std::abs(__max) < std::abs((*this)[__r][__c]))
            __max = (*this)[__pos = __r][__c];
      }

      else if (__max ==
               static_cast<value_type>(
                   0))  // Find the first non-zero element in the column.
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows(); ++__r)
          if ((__max = (*this)[__r][__c]) != static_cast<value_type>(0)) {
            __pos = __r;
            break;
          }

      if (__pos != __pivots.size()) {  // Swap 2 rows.
        sgn = !sgn;
        std::swap(__perm[__pos], __perm[__pivots.size()]);
        std::swap((*this)[__pos], (*this)[__pivots.size()]);
      }

      if (__max != static_cast<value_type>(0)) {  // Forward elimination
        for (size_type __r = __pivots.size() + 1; __r != _Matrix::rows();
             ++__r) {
          auto __m = (*this)[__r][__c] / __max;
          (*this)[__r][__c] = 0;
          (*this)[__r][__pivots.size()] = __m;

          for (size_type __i = __c + 1; __i != _Matrix::cols(); ++__i)
            (*this)[__r][__i] -= (*this)[__pivots.size()][__i] * __m;
        }

        __pivots.emplace_back(__c);
      }
    }
  }

 public:
  size_type rank() const { return __pivots.size(); }

  value_type det() const {
    assert(_Matrix::rows() == _Matrix::cols());

    value_type __d = sgn ? -1 : 1;
    for (size_type __i = 0; __i != _Matrix::rows(); ++__i)
      __d *= (*this)[__i][__i];
    return __d;
  }

  _Matrix lower() const;

  _Matrix upper() const;

  _Matrix inverse() const {
    assert(_Matrix::rows() == _Matrix::cols());

    _Matrix __inv;
    // add solve(e_i) to __inv for i=0, ..., rows()
    return __inv;
  }

  // O(dim(ker) * size)
  _Matrix kernel() const {
    _Matrix __ker(_Matrix::cols() - rank(), _Matrix::cols());

    for (size_type __c = 0, __i = 0; __c != _Matrix::cols(); ++__c) {
      if (__i != __pivots.size() && __pivots[__i] == __c) {
        ++__i;
        continue;
      }

      auto &__v = __ker[__c - __i];
      __v[__c] = 1;
      for (size_type __r = 0; __r != __i; ++__r) __v[__r] = -(*this)[__r][__c];

      for (size_type __j = __i; __j--;) {
        auto __x = __v[__j] / (*this)[__j][__pivots[__j]];
        __v[__j] = 0;
        __v[__pivots[__j]] = __x;

        for (size_type __r = 0; __r != __j; ++__r)
          __v[__r] -= (*this)[__r][__pivots[__j]] * __x;
      }
    }

    return __ker;
  }

  template <class _Vec> std::pair<bool, _Vec> solve(const _Vec &__b) const {
    assert(!(__b.size() < _Matrix::rows()));

    // Solution
    _Vec __y(_Matrix::rows()), __x(_Matrix::cols());

    // Backward substitution with L
    for (size_type __c = 0; __c != _Matrix::rows(); ++__c) {
      __y[__c] += __b[__perm[__c]];

      for (size_type __r = __c + 1; __r != _Matrix::rows(); ++__r)
        __y[__r] -= __y[__c] * (*this)[__r][__c];
    }

    // Backward substitution with U
    for (size_type __i = rank(); __i != _Matrix::rows(); ++__i)
      if (__y[__i] != static_cast<value_type>(0))
        return std::make_pair(false, __x);

    for (size_type __i = rank(); __i--;) {
      auto __c = __pivots[__i];

      __x[__c] = __y[__i] / (*this)[__i][__c];

      for (size_type __r = 0; __r != __i; ++__r)
        __y[__r] -= __x[__c] * (*this)[__r][__c];
    }

    return std::make_pair(true, __x);
  }
};

}  // namespace workspace

Back to top page