17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder> 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct packed_triangular_matrix_vector_product; 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor> 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IsLower = (Mode & Lower) ==Lower, 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HasUnitDiag = (Mode & UnitDiag)==UnitDiag, 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::conj_if<ConjRhs> cj; 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap; 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i=0; i<size; ++i) 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index s = IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index r = IsLower ? size-i: i+1; 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ResMap(res+(IsLower ? s+i : 0),r) += alpha * cj(rhs[i]) * ConjLhsType(LhsMap(lhs+s,r)); 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (HasUnitDiag) 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res[i] += alpha * cj(rhs[i]); 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lhs += IsLower ? size-i: i+1; 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs> 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct packed_triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor> 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IsLower = (Mode & Lower) ==Lower, 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HasUnitDiag = (Mode & UnitDiag)==UnitDiag, 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HasZeroDiag = (Mode & ZeroDiag)==ZeroDiag 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Index size, const LhsScalar* lhs, const RhsScalar* rhs, ResScalar* res, ResScalar alpha) 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::conj_if<ConjRhs> cj; 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename conj_expr_if<ConjLhs,LhsMap>::type ConjLhsType; 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap; 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename conj_expr_if<ConjRhs,RhsMap>::type ConjRhsType; 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i=0; i<size; ++i) 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index s = !IsLower&&(HasUnitDiag||HasZeroDiag) ? 1 : 0; 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index r = IsLower ? i+1 : size-i; 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (EIGEN_IMPLIES(HasUnitDiag||HasZeroDiag, (--r)>0)) 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res[i] += alpha * (ConjLhsType(LhsMap(lhs+s,r)).cwiseProduct(ConjRhsType(RhsMap(rhs+(IsLower ? 0 : s+i),r)))).sum(); 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (HasUnitDiag) 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res[i] += alpha * cj(rhs[i]); 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lhs += IsLower ? i+1 : size-i; 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_PACKED_TRIANGULAR_MATRIX_VECTOR_H 80