| Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame^] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> |
| 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_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
| 11 | #define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
| 12 | |
| 13 | namespace internal { |
| 14 | |
| 15 | template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> |
| 16 | struct packed_triangular_solve_vector; |
| 17 | |
| 18 | // forward and backward substitution, row-major, rhs is a vector |
| 19 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
| 20 | struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> |
| 21 | { |
| 22 | enum { |
| 23 | IsLower = (Mode&Lower)==Lower |
| 24 | }; |
| 25 | static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
| 26 | { |
| 27 | internal::conj_if<Conjugate> cj; |
| 28 | typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; |
| 29 | typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; |
| 30 | |
| 31 | lhs += IsLower ? 0 : (size*(size+1)>>1)-1; |
| 32 | for(Index pi=0; pi<size; ++pi) |
| 33 | { |
| 34 | Index i = IsLower ? pi : size-pi-1; |
| 35 | Index s = IsLower ? 0 : 1; |
| 36 | if (pi>0) |
| 37 | rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi)) |
| 38 | .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum(); |
| 39 | if (!(Mode & UnitDiag)) |
| 40 | rhs[i] /= cj(lhs[IsLower ? i : 0]); |
| 41 | IsLower ? lhs += pi+1 : lhs -= pi+2; |
| 42 | } |
| 43 | } |
| 44 | }; |
| 45 | |
| 46 | // forward and backward substitution, column-major, rhs is a vector |
| 47 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
| 48 | struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> |
| 49 | { |
| 50 | enum { |
| 51 | IsLower = (Mode&Lower)==Lower |
| 52 | }; |
| 53 | static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
| 54 | { |
| 55 | internal::conj_if<Conjugate> cj; |
| 56 | typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; |
| 57 | typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; |
| 58 | |
| 59 | lhs += IsLower ? 0 : size*(size-1)>>1; |
| 60 | for(Index pi=0; pi<size; ++pi) |
| 61 | { |
| 62 | Index i = IsLower ? pi : size-pi-1; |
| 63 | Index r = size - pi - 1; |
| 64 | if (!(Mode & UnitDiag)) |
| 65 | rhs[i] /= cj(lhs[IsLower ? 0 : i]); |
| 66 | if (r>0) |
| 67 | Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -= |
| 68 | rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r)); |
| 69 | IsLower ? lhs += size-pi : lhs -= r; |
| 70 | } |
| 71 | } |
| 72 | }; |
| 73 | |
| 74 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> |
| 75 | struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> |
| 76 | { |
| 77 | static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
| 78 | { |
| 79 | packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft, |
| 80 | ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), |
| 81 | Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor |
| 82 | >::run(size, lhs, rhs); |
| 83 | } |
| 84 | }; |
| 85 | |
| 86 | } // end namespace internal |
| 87 | |
| 88 | #endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |