blob: 014af24e4af1a9a6f570a5fc89b5be1a8026b6a9 [file] [log] [blame]
// 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/.
#ifndef EIGEN_BAND_TRIANGULARSOLVER_H
#define EIGEN_BAND_TRIANGULARSOLVER_H
namespace internal {
/* \internal
* Solve Ax=b with A a band triangular matrix
* TODO: extend it to matrices for x abd b */
template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
struct band_solve_triangular_selector;
template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
struct band_solve_triangular_selector<Index, Mode, LhsScalar, ConjLhs, RhsScalar, RowMajor> {
typedef Map<const Matrix<LhsScalar, Dynamic, Dynamic, RowMajor>, 0, OuterStride<> > LhsMap;
typedef Map<Matrix<RhsScalar, Dynamic, 1> > RhsMap;
enum { IsLower = (Mode & Lower) ? 1 : 0 };
static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) {
const LhsMap lhs(_lhs, size, k + 1, OuterStride<>(lhsStride));
RhsMap other(_other, size, 1);
std::conditional_t<ConjLhs, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
const LhsMap&>
cjLhs(lhs);
for (int col = 0; col < other.cols(); ++col) {
for (int ii = 0; ii < size; ++ii) {
int i = IsLower ? ii : size - ii - 1;
int actual_k = (std::min)(k, ii);
int actual_start = IsLower ? k - actual_k : 1;
if (actual_k > 0)
other.coeffRef(i, col) -= cjLhs.row(i)
.segment(actual_start, actual_k)
.transpose()
.cwiseProduct(other.col(col).segment(IsLower ? i - actual_k : i + 1, actual_k))
.sum();
if ((Mode & UnitDiag) == 0) other.coeffRef(i, col) /= cjLhs(i, IsLower ? k : 0);
}
}
}
};
template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
struct band_solve_triangular_selector<Index, Mode, LhsScalar, ConjLhs, RhsScalar, ColMajor> {
typedef Map<const Matrix<LhsScalar, Dynamic, Dynamic, ColMajor>, 0, OuterStride<> > LhsMap;
typedef Map<Matrix<RhsScalar, Dynamic, 1> > RhsMap;
enum { IsLower = (Mode & Lower) ? 1 : 0 };
static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) {
const LhsMap lhs(_lhs, k + 1, size, OuterStride<>(lhsStride));
RhsMap other(_other, size, 1);
std::conditional_t<ConjLhs, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
const LhsMap&>
cjLhs(lhs);
for (int col = 0; col < other.cols(); ++col) {
for (int ii = 0; ii < size; ++ii) {
int i = IsLower ? ii : size - ii - 1;
int actual_k = (std::min)(k, size - ii - 1);
int actual_start = IsLower ? 1 : k - actual_k;
if ((Mode & UnitDiag) == 0) other.coeffRef(i, col) /= cjLhs(IsLower ? 0 : k, i);
if (actual_k > 0)
other.col(col).segment(IsLower ? i + 1 : i - actual_k, actual_k) -=
other.coeff(i, col) * cjLhs.col(i).segment(actual_start, actual_k);
}
}
}
};
} // end namespace internal
#endif // EIGEN_BAND_TRIANGULARSOLVER_H