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 89