blob: f9bfdc1d46d6cfa377c16b72974d8c1a13b711e3 [file] [log] [blame]
Googler45874d82019-08-21 12:06:47 -07001// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_BAND_TRIANGULARSOLVER_H
11#define EIGEN_BAND_TRIANGULARSOLVER_H
12
Rasmus Munk Larsen04037502024-04-09 13:58:59 -070013namespace Eigen {
Googler45874d82019-08-21 12:06:47 -070014namespace internal {
15
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080016/* \internal
17 * Solve Ax=b with A a band triangular matrix
18 * TODO: extend it to matrices for x abd b */
19template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
Googler45874d82019-08-21 12:06:47 -070020struct band_solve_triangular_selector;
21
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080022template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
23struct band_solve_triangular_selector<Index, Mode, LhsScalar, ConjLhs, RhsScalar, RowMajor> {
24 typedef Map<const Matrix<LhsScalar, Dynamic, Dynamic, RowMajor>, 0, OuterStride<> > LhsMap;
25 typedef Map<Matrix<RhsScalar, Dynamic, 1> > RhsMap;
26 enum { IsLower = (Mode & Lower) ? 1 : 0 };
27 static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) {
28 const LhsMap lhs(_lhs, size, k + 1, OuterStride<>(lhsStride));
29 RhsMap other(_other, size, 1);
30 std::conditional_t<ConjLhs, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
31 const LhsMap&>
32 cjLhs(lhs);
Googler45874d82019-08-21 12:06:47 -070033
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080034 for (int col = 0; col < other.cols(); ++col) {
35 for (int ii = 0; ii < size; ++ii) {
36 int i = IsLower ? ii : size - ii - 1;
37 int actual_k = (std::min)(k, ii);
38 int actual_start = IsLower ? k - actual_k : 1;
Googler45874d82019-08-21 12:06:47 -070039
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080040 if (actual_k > 0)
41 other.coeffRef(i, col) -= cjLhs.row(i)
42 .segment(actual_start, actual_k)
43 .transpose()
44 .cwiseProduct(other.col(col).segment(IsLower ? i - actual_k : i + 1, actual_k))
45 .sum();
Googler45874d82019-08-21 12:06:47 -070046
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080047 if ((Mode & UnitDiag) == 0) other.coeffRef(i, col) /= cjLhs(i, IsLower ? k : 0);
Googler45874d82019-08-21 12:06:47 -070048 }
49 }
50 }
51};
52
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080053template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
54struct band_solve_triangular_selector<Index, Mode, LhsScalar, ConjLhs, RhsScalar, ColMajor> {
55 typedef Map<const Matrix<LhsScalar, Dynamic, Dynamic, ColMajor>, 0, OuterStride<> > LhsMap;
56 typedef Map<Matrix<RhsScalar, Dynamic, 1> > RhsMap;
57 enum { IsLower = (Mode & Lower) ? 1 : 0 };
58 static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) {
59 const LhsMap lhs(_lhs, k + 1, size, OuterStride<>(lhsStride));
60 RhsMap other(_other, size, 1);
61 std::conditional_t<ConjLhs, const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>, LhsMap>,
62 const LhsMap&>
63 cjLhs(lhs);
Googler45874d82019-08-21 12:06:47 -070064
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080065 for (int col = 0; col < other.cols(); ++col) {
66 for (int ii = 0; ii < size; ++ii) {
67 int i = IsLower ? ii : size - ii - 1;
68 int actual_k = (std::min)(k, size - ii - 1);
69 int actual_start = IsLower ? 1 : k - actual_k;
Googler45874d82019-08-21 12:06:47 -070070
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080071 if ((Mode & UnitDiag) == 0) other.coeffRef(i, col) /= cjLhs(IsLower ? 0 : k, i);
72
73 if (actual_k > 0)
74 other.col(col).segment(IsLower ? i + 1 : i - actual_k, actual_k) -=
75 other.coeff(i, col) * cjLhs.col(i).segment(actual_start, actual_k);
76 }
77 }
78 }
79};
80
Rasmus Munk Larsen04037502024-04-09 13:58:59 -070081} // namespace internal
82} // namespace Eigen
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080083
84#endif // EIGEN_BAND_TRIANGULARSOLVER_H