IncompleteLU.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_INCOMPLETE_LU_H
#define EIGEN_INCOMPLETE_LU_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename Scalar_>
class IncompleteLU : public SparseSolverBase<IncompleteLU<Scalar_> > {
 protected:
  using Base = SparseSolverBase<IncompleteLU<Scalar_>>;
  using Base::m_isInitialized;
 
  using Scalar = Scalar_;
  using Vector = Matrix<Scalar, Dynamic, 1>;
  using Index = typename Vector::Index;
  using FactorType = SparseMatrix<Scalar, RowMajor>;
 
 public:
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic>;
 
  IncompleteLU() {}
 
  template <typename MatrixType>
  IncompleteLU(const MatrixType& mat) {
    compute(mat);
  }
 
  Index rows() const { return m_lu.rows(); }
  Index cols() const { return m_lu.cols(); }
 
  template <typename MatrixType>
  IncompleteLU& compute(const MatrixType& mat) {
    m_lu = mat;
    Index size = mat.cols();
    Vector diag(size);
    for (Index i = 0; i < size; ++i) {
      typename FactorType::InnerIterator k_it(m_lu, i);
      for (; k_it && k_it.index() < i; ++k_it) {
        Index k = k_it.index();
        k_it.valueRef() /= diag(k);
 
        typename FactorType::InnerIterator j_it(k_it);
        typename FactorType::InnerIterator kj_it(m_lu, k);
        while (kj_it && kj_it.index() <= k) ++kj_it;
        for (++j_it; j_it;) {
          if (kj_it.index() == j_it.index()) {
            j_it.valueRef() -= k_it.value() * kj_it.value();
            ++j_it;
            ++kj_it;
          } else if (kj_it.index() < j_it.index())
            ++kj_it;
          else
            ++j_it;
        }
      }
      if (k_it && k_it.index() == i)
        diag(i) = k_it.value();
      else
        diag(i) = 1;
    }
    m_isInitialized = true;
    return *this;
  }
 
  template <typename Rhs, typename Dest>
  void _solve_impl(const Rhs& b, Dest& x) const {
    x = m_lu.template triangularView<UnitLower>().solve(b);
    x = m_lu.template triangularView<Upper>().solve(x);
  }
 
 protected:
  FactorType m_lu;
};
 
}  // end namespace Eigen
 
#endif  // EIGEN_INCOMPLETE_LU_H