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_MATRIX_VECTOR_H
11#define EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H
12
13namespace internal {
14
15template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder>
16struct packed_triangular_matrix_vector_product;
17
18template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
19struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor>
20{
21  typedef typename ScalarBinaryOpTraits<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
47template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
48struct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor>
49{
50  typedef typename ScalarBinaryOpTraits<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