1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) inria.fr> 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_TRIANGULAR_SOLVER_VECTOR_H 11 #define EIGEN_TRIANGULAR_SOLVER_VECTOR_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> 18 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> 19 { 20 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) 21 { 22 triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft, 23 ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), 24 Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor 25 >::run(size, _lhs, lhsStride, rhs); 26 } 27 }; 28 29 // forward and backward substitution, row-major, rhs is a vector 30 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> 31 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> 32 { 33 enum { 34 IsLower = ((Mode&Lower)==Lower) 35 }; 36 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) 37 { 38 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap; 39 const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride)); 40 typename internal::conditional< 41 Conjugate, 42 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>, 43 const LhsMap&> 44 ::type cjLhs(lhs); 45 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; 46 for(Index pi=IsLower ? 0 : size; 47 IsLower ? pi<size : pi>0; 48 IsLower ? pi+=PanelWidth : pi-=PanelWidth) 49 { 50 Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth); 51 52 Index r = IsLower ? pi : size - pi; // remaining size 53 if (r > 0) 54 { 55 // let's directly call the low level product function because: 56 // 1 - it is faster to compile 57 // 2 - it is slighlty faster at runtime 58 Index startRow = IsLower ? pi : pi-actualPanelWidth; 59 Index startCol = IsLower ? 0 : pi; 60 61 general_matrix_vector_product<Index,LhsScalar,RowMajor,Conjugate,RhsScalar,false>::run( 62 actualPanelWidth, r, 63 &lhs.coeffRef(startRow,startCol), lhsStride, 64 rhs + startCol, 1, 65 rhs + startRow, 1, 66 RhsScalar(-1)); 67 } 68 69 for(Index k=0; k<actualPanelWidth; ++k) 70 { 71 Index i = IsLower ? pi+k : pi-k-1; 72 Index s = IsLower ? pi : i+1; 73 if (k>0) 74 rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum(); 75 76 if(!(Mode & UnitDiag)) 77 rhs[i] /= cjLhs(i,i); 78 } 79 } 80 } 81 }; 82 83 // forward and backward substitution, column-major, rhs is a vector 84 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> 85 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> 86 { 87 enum { 88 IsLower = ((Mode&Lower)==Lower) 89 }; 90 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs) 91 { 92 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap; 93 const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride)); 94 typename internal::conditional<Conjugate, 95 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>, 96 const LhsMap& 97 >::type cjLhs(lhs); 98 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; 99 100 for(Index pi=IsLower ? 0 : size; 101 IsLower ? pi<size : pi>0; 102 IsLower ? pi+=PanelWidth : pi-=PanelWidth) 103 { 104 Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth); 105 Index startBlock = IsLower ? pi : pi-actualPanelWidth; 106 Index endBlock = IsLower ? pi + actualPanelWidth : 0; 107 108 for(Index k=0; k<actualPanelWidth; ++k) 109 { 110 Index i = IsLower ? pi+k : pi-k-1; 111 if(!(Mode & UnitDiag)) 112 rhs[i] /= cjLhs.coeff(i,i); 113 114 Index r = actualPanelWidth - k - 1; // remaining size 115 Index s = IsLower ? i+1 : i-r; 116 if (r>0) 117 Map<Matrix<RhsScalar,Dynamic,1> >(rhs+s,r) -= rhs[i] * cjLhs.col(i).segment(s,r); 118 } 119 Index r = IsLower ? size - endBlock : startBlock; // remaining size 120 if (r > 0) 121 { 122 // let's directly call the low level product function because: 123 // 1 - it is faster to compile 124 // 2 - it is slighlty faster at runtime 125 general_matrix_vector_product<Index,LhsScalar,ColMajor,Conjugate,RhsScalar,false>::run( 126 r, actualPanelWidth, 127 &lhs.coeffRef(endBlock,startBlock), lhsStride, 128 rhs+startBlock, 1, 129 rhs+endBlock, 1, RhsScalar(-1)); 130 } 131 } 132 } 133 }; 134 135 } // end namespace internal 136 137 } // end namespace Eigen 138 139 #endif // EIGEN_TRIANGULAR_SOLVER_VECTOR_H 140
