| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> |
| // |
| // 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_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
| #define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
| |
| namespace internal { |
| |
| template <typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> |
| struct packed_triangular_solve_vector; |
| |
| // forward and backward substitution, row-major, rhs is a vector |
| template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
| struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> { |
| enum { IsLower = (Mode & Lower) == Lower }; |
| static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) { |
| internal::conj_if<Conjugate> cj; |
| typedef Map<const Matrix<LhsScalar, Dynamic, 1> > LhsMap; |
| typedef typename conj_expr_if<Conjugate, LhsMap>::type ConjLhsType; |
| |
| lhs += IsLower ? 0 : (size * (size + 1) >> 1) - 1; |
| for (Index pi = 0; pi < size; ++pi) { |
| Index i = IsLower ? pi : size - pi - 1; |
| Index s = IsLower ? 0 : 1; |
| if (pi > 0) |
| rhs[i] -= (ConjLhsType(LhsMap(lhs + s, pi)) |
| .cwiseProduct(Map<const Matrix<RhsScalar, Dynamic, 1> >(rhs + (IsLower ? 0 : i + 1), pi))) |
| .sum(); |
| if (!(Mode & UnitDiag)) rhs[i] /= cj(lhs[IsLower ? i : 0]); |
| IsLower ? lhs += pi + 1 : lhs -= pi + 2; |
| } |
| } |
| }; |
| |
| // forward and backward substitution, column-major, rhs is a vector |
| template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
| struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> { |
| enum { IsLower = (Mode & Lower) == Lower }; |
| static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) { |
| internal::conj_if<Conjugate> cj; |
| typedef Map<const Matrix<LhsScalar, Dynamic, 1> > LhsMap; |
| typedef typename conj_expr_if<Conjugate, LhsMap>::type ConjLhsType; |
| |
| lhs += IsLower ? 0 : size * (size - 1) >> 1; |
| for (Index pi = 0; pi < size; ++pi) { |
| Index i = IsLower ? pi : size - pi - 1; |
| Index r = size - pi - 1; |
| if (!(Mode & UnitDiag)) rhs[i] /= cj(lhs[IsLower ? 0 : i]); |
| if (r > 0) |
| Map<Matrix<RhsScalar, Dynamic, 1> >(rhs + (IsLower ? i + 1 : 0), r) -= |
| rhs[i] * ConjLhsType(LhsMap(lhs + (IsLower ? 1 : 0), r)); |
| IsLower ? lhs += size - pi : lhs -= r; |
| } |
| } |
| }; |
| |
| template <typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> |
| struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> { |
| static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) { |
| packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, |
| ((Mode & Upper) == Upper ? Lower : Upper) | (Mode & UnitDiag), Conjugate, |
| StorageOrder == RowMajor ? ColMajor : RowMajor>::run(size, lhs, rhs); |
| } |
| }; |
| |
| } // end namespace internal |
| |
| #endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |