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