TriangularSolverVector.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 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_TRIANGULAR_SOLVER_VECTOR_H
#define EIGEN_TRIANGULAR_SOLVER_VECTOR_H
 
// IWYU pragma: private
#include "../InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
 
template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> {
  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) {
    triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft,
                            ((Mode & Upper) == Upper ? Lower : Upper) | (Mode & UnitDiag), Conjugate,
                            StorageOrder == RowMajor ? ColMajor : RowMajor>::run(size, _lhs, lhsStride, rhs);
  }
};
 
// forward and backward substitution, row-major, rhs is a vector
template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> {
  enum { IsLower = ((Mode & Lower) == Lower) };
  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) {
    using LhsMap = Map<const Matrix<LhsScalar, Dynamic, Dynamic, RowMajor>, 0, OuterStride<> >;
    const LhsMap lhs(_lhs, size, size, OuterStride<>(lhsStride));
 
    using LhsMapper = const_blas_data_mapper<LhsScalar, Index, RowMajor>;
    using RhsMapper = const_blas_data_mapper<RhsScalar, Index, ColMajor>;
 
    std::conditional_t<Conjugate, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
                       const LhsMap&>
        cjLhs(lhs);
    static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
    for (Index pi = IsLower ? 0 : size; IsLower ? pi < size : pi > 0; IsLower ? pi += PanelWidth : pi -= PanelWidth) {
      Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
 
      Index r = IsLower ? pi : size - pi;  // remaining size
      if (r > 0) {
        // let's directly call the low level product function because:
        // 1 - it is faster to compile
        // 2 - it is slightly faster at runtime
        Index startRow = IsLower ? pi : pi - actualPanelWidth;
        Index startCol = IsLower ? 0 : pi;
 
        general_matrix_vector_product<Index, LhsScalar, LhsMapper, RowMajor, Conjugate, RhsScalar, RhsMapper,
                                      false>::run(actualPanelWidth, r,
                                                  LhsMapper(&lhs.coeffRef(startRow, startCol), lhsStride),
                                                  RhsMapper(rhs + startCol, 1), rhs + startRow, 1, RhsScalar(-1));
      }
 
      for (Index k = 0; k < actualPanelWidth; ++k) {
        Index i = IsLower ? pi + k : pi - k - 1;
        Index s = IsLower ? pi : i + 1;
        if (k > 0)
          rhs[i] -= (cjLhs.row(i).segment(s, k).transpose().cwiseProduct(
                         Map<const Matrix<RhsScalar, Dynamic, 1> >(rhs + s, k)))
                        .sum();
 
        EIGEN_IF_CONSTEXPR (!(Mode & UnitDiag)) {
          if (!is_identically_zero(rhs[i])) rhs[i] /= cjLhs(i, i);
        }
      }
    }
  }
};
 
// forward and backward substitution, column-major, rhs is a vector
template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> {
  enum { IsLower = ((Mode & Lower) == Lower) };
  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) {
    using LhsMap = Map<const Matrix<LhsScalar, Dynamic, Dynamic, ColMajor>, 0, OuterStride<> >;
    const LhsMap lhs(_lhs, size, size, OuterStride<>(lhsStride));
    using LhsMapper = const_blas_data_mapper<LhsScalar, Index, ColMajor>;
    using RhsMapper = const_blas_data_mapper<RhsScalar, Index, ColMajor>;
    std::conditional_t<Conjugate, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
                       const LhsMap&>
        cjLhs(lhs);
    static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
 
    for (Index pi = IsLower ? 0 : size; IsLower ? pi < size : pi > 0; IsLower ? pi += PanelWidth : pi -= PanelWidth) {
      Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
      Index startBlock = IsLower ? pi : pi - actualPanelWidth;
      Index endBlock = IsLower ? pi + actualPanelWidth : 0;
 
      for (Index k = 0; k < actualPanelWidth; ++k) {
        Index i = IsLower ? pi + k : pi - k - 1;
        if (!is_identically_zero(rhs[i])) {
          EIGEN_IF_CONSTEXPR (!(Mode & UnitDiag)) rhs[i] /= cjLhs.coeff(i, i);
 
          Index r = actualPanelWidth - k - 1;  // remaining size
          Index s = IsLower ? i + 1 : i - r;
          if (r > 0) Map<Matrix<RhsScalar, Dynamic, 1> >(rhs + s, r) -= rhs[i] * cjLhs.col(i).segment(s, r);
        }
      }
      Index r = IsLower ? size - endBlock : startBlock;  // remaining size
      if (r > 0) {
        // let's directly call the low level product function because:
        // 1 - it is faster to compile
        // 2 - it is slightly faster at runtime
        general_matrix_vector_product<Index, LhsScalar, LhsMapper, ColMajor, Conjugate, RhsScalar, RhsMapper,
                                      false>::run(r, actualPanelWidth,
                                                  LhsMapper(&lhs.coeffRef(endBlock, startBlock), lhsStride),
                                                  RhsMapper(rhs + startBlock, 1), rhs + endBlock, 1, RhsScalar(-1));
      }
    }
  }
};
 
}  // end namespace internal
 
}  // end namespace Eigen
 
#endif  // EIGEN_TRIANGULAR_SOLVER_VECTOR_H