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_MATRIX_VECTOR_H 11*bf2c3715SXin Li #define EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H 12*bf2c3715SXin Li 13*bf2c3715SXin Li namespace internal { 14*bf2c3715SXin Li 15*bf2c3715SXin Li template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder> 16*bf2c3715SXin Li struct packed_triangular_matrix_vector_product; 17*bf2c3715SXin Li 18*bf2c3715SXin Li template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> 19*bf2c3715SXin Li struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor> 20*bf2c3715SXin Li { 21*bf2c3715SXin Li typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar; 22*bf2c3715SXin Li enum { 23*bf2c3715SXin Li IsLower = (Mode & Lower) ==Lower, 24*bf2c3715SXin Li HasUnitDiag = (Mode & UnitDiag)==UnitDiag, 25*bf2c3715SXin Li HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag 26*bf2c3715SXin Li }; 27*bf2c3715SXin Li static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) 28*bf2c3715SXin Li { 29*bf2c3715SXin Li internal::conj_if<ConjRhs> cj; 30*bf2c3715SXin Li typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 31*bf2c3715SXin Li typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; 32*bf2c3715SXin Li typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap; 33*bf2c3715SXin Li 34*bf2c3715SXin Li for (Index i=0; i<size; ++i) 35*bf2c3715SXin Li { 36*bf2c3715SXin Li Index s = IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; 37*bf2c3715SXin Li Index r = IsLower ? size-i: i+1; 38*bf2c3715SXin Li if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) 39*bf2c3715SXin Li ResMap(res+(IsLower ? s+i : 0),r) += alpha * cj(rhs[i]) * ConjLhsType(LhsMap(lhs+s,r)); 40*bf2c3715SXin Li if (HasUnitDiag) 41*bf2c3715SXin Li res[i] += alpha * cj(rhs[i]); 42*bf2c3715SXin Li lhs += IsLower ? size-i: i+1; 43*bf2c3715SXin Li } 44*bf2c3715SXin Li }; 45*bf2c3715SXin Li }; 46*bf2c3715SXin Li 47*bf2c3715SXin Li template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> 48*bf2c3715SXin Li struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor> 49*bf2c3715SXin Li { 50*bf2c3715SXin Li typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar; 51*bf2c3715SXin Li enum { 52*bf2c3715SXin Li IsLower = (Mode & Lower) ==Lower, 53*bf2c3715SXin Li HasUnitDiag = (Mode & UnitDiag)==UnitDiag, 54*bf2c3715SXin Li HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag 55*bf2c3715SXin Li }; 56*bf2c3715SXin Li static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) 57*bf2c3715SXin Li { 58*bf2c3715SXin Li internal::conj_if<ConjRhs> cj; 59*bf2c3715SXin Li typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 60*bf2c3715SXin Li typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; 61*bf2c3715SXin Li typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap; 62*bf2c3715SXin Li typedef typename conj_expr_if<ConjRhs,RhsMap>::type ConjRhsType; 63*bf2c3715SXin Li 64*bf2c3715SXin Li for (Index i=0; i<size; ++i) 65*bf2c3715SXin Li { 66*bf2c3715SXin Li Index s = !IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; 67*bf2c3715SXin Li Index r = IsLower ? i+1 : size-i; 68*bf2c3715SXin Li if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) 69*bf2c3715SXin Li res[i] += alpha * (ConjLhsType(LhsMap(lhs+s,r)).cwiseProduct(ConjRhsType(RhsMap(rhs+(IsLower ? 0 : s+i),r)))).sum(); 70*bf2c3715SXin Li if (HasUnitDiag) 71*bf2c3715SXin Li res[i] += alpha * cj(rhs[i]); 72*bf2c3715SXin Li lhs += IsLower ? i+1 : size-i; 73*bf2c3715SXin Li } 74*bf2c3715SXin Li }; 75*bf2c3715SXin Li }; 76*bf2c3715SXin Li 77*bf2c3715SXin Li } // end namespace internal 78*bf2c3715SXin Li 79*bf2c3715SXin Li #endif // EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H 80