Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 1 | // 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 Larsen | 0403750 | 2024-04-09 13:58:59 -0700 | [diff] [blame] | 13 | namespace Eigen { |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 14 | namespace internal { |
| 15 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 16 | /* \internal |
| 17 | * Solve Ax=b with A a band triangular matrix |
| 18 | * TODO: extend it to matrices for x abd b */ |
| 19 | template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder> |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 20 | struct band_solve_triangular_selector; |
| 21 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 22 | template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar> |
| 23 | struct 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); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 33 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 34 | 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; |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 39 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 40 | 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(); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 46 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 47 | if ((Mode & UnitDiag) == 0) other.coeffRef(i, col) /= cjLhs(i, IsLower ? k : 0); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 48 | } |
| 49 | } |
| 50 | } |
| 51 | }; |
| 52 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 53 | template <typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar> |
| 54 | struct 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); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 64 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 65 | 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; |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 70 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 71 | 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 Larsen | 0403750 | 2024-04-09 13:58:59 -0700 | [diff] [blame] | 81 | } // namespace internal |
| 82 | } // namespace Eigen |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 83 | |
| 84 | #endif // EIGEN_BAND_TRIANGULARSOLVER_H |