1*bf2c3715SXin Li // This file is part of Eigen, a lightweight C++ template library 2*bf2c3715SXin Li // for linear algebra. 3*bf2c3715SXin Li // 4*bf2c3715SXin Li // Copyright (C) 2012 Chen-Pang He <[email protected]> 5*bf2c3715SXin Li // 6*bf2c3715SXin Li // This Source Code Form is subject to the terms of the Mozilla 7*bf2c3715SXin Li // Public License v. 2.0. If a copy of the MPL was not distributed 8*bf2c3715SXin Li // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9*bf2c3715SXin Li 10*bf2c3715SXin Li #ifndef EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H 11*bf2c3715SXin Li #define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H 12*bf2c3715SXin Li 13*bf2c3715SXin Li namespace internal { 14*bf2c3715SXin Li 15*bf2c3715SXin Li template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> 16*bf2c3715SXin Li struct packed_triangular_solve_vector; 17*bf2c3715SXin Li 18*bf2c3715SXin Li // forward and backward substitution, row-major, rhs is a vector 19*bf2c3715SXin Li template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> 20*bf2c3715SXin Li struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> 21*bf2c3715SXin Li { 22*bf2c3715SXin Li enum { 23*bf2c3715SXin Li IsLower = (Mode&Lower)==Lower 24*bf2c3715SXin Li }; 25*bf2c3715SXin Li static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) 26*bf2c3715SXin Li { 27*bf2c3715SXin Li internal::conj_if<Conjugate> cj; 28*bf2c3715SXin Li typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 29*bf2c3715SXin Li typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; 30*bf2c3715SXin Li 31*bf2c3715SXin Li lhs += IsLower ? 0 : (size*(size+1)>>1)-1; 32*bf2c3715SXin Li for(Index pi=0; pi<size; ++pi) 33*bf2c3715SXin Li { 34*bf2c3715SXin Li Index i = IsLower ? pi : size-pi-1; 35*bf2c3715SXin Li Index s = IsLower ? 0 : 1; 36*bf2c3715SXin Li if (pi>0) 37*bf2c3715SXin Li rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi)) 38*bf2c3715SXin Li .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum(); 39*bf2c3715SXin Li if (!(Mode & UnitDiag)) 40*bf2c3715SXin Li rhs[i] /= cj(lhs[IsLower ? i : 0]); 41*bf2c3715SXin Li IsLower ? lhs += pi+1 : lhs -= pi+2; 42*bf2c3715SXin Li } 43*bf2c3715SXin Li } 44*bf2c3715SXin Li }; 45*bf2c3715SXin Li 46*bf2c3715SXin Li // forward and backward substitution, column-major, rhs is a vector 47*bf2c3715SXin Li template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> 48*bf2c3715SXin Li struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> 49*bf2c3715SXin Li { 50*bf2c3715SXin Li enum { 51*bf2c3715SXin Li IsLower = (Mode&Lower)==Lower 52*bf2c3715SXin Li }; 53*bf2c3715SXin Li static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) 54*bf2c3715SXin Li { 55*bf2c3715SXin Li internal::conj_if<Conjugate> cj; 56*bf2c3715SXin Li typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 57*bf2c3715SXin Li typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; 58*bf2c3715SXin Li 59*bf2c3715SXin Li lhs += IsLower ? 0 : size*(size-1)>>1; 60*bf2c3715SXin Li for(Index pi=0; pi<size; ++pi) 61*bf2c3715SXin Li { 62*bf2c3715SXin Li Index i = IsLower ? pi : size-pi-1; 63*bf2c3715SXin Li Index r = size - pi - 1; 64*bf2c3715SXin Li if (!(Mode & UnitDiag)) 65*bf2c3715SXin Li rhs[i] /= cj(lhs[IsLower ? 0 : i]); 66*bf2c3715SXin Li if (r>0) 67*bf2c3715SXin Li Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -= 68*bf2c3715SXin Li rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r)); 69*bf2c3715SXin Li IsLower ? lhs += size-pi : lhs -= r; 70*bf2c3715SXin Li } 71*bf2c3715SXin Li } 72*bf2c3715SXin Li }; 73*bf2c3715SXin Li 74*bf2c3715SXin Li template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> 75*bf2c3715SXin Li struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> 76*bf2c3715SXin Li { 77*bf2c3715SXin Li static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) 78*bf2c3715SXin Li { 79*bf2c3715SXin Li packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft, 80*bf2c3715SXin Li ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), 81*bf2c3715SXin Li Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor 82*bf2c3715SXin Li >::run(size, lhs, rhs); 83*bf2c3715SXin Li } 84*bf2c3715SXin Li }; 85*bf2c3715SXin Li 86*bf2c3715SXin Li } // end namespace internal 87*bf2c3715SXin Li 88*bf2c3715SXin Li #endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H 89