TriangularSolverVector.h
Go to the documentation of this file.
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@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_TRIANGULAR_SOLVER_VECTOR_H
26 #define EIGEN_TRIANGULAR_SOLVER_VECTOR_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
32 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
33 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
34 {
35  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
36  {
37  triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
38  ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
39  Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
40  >::run(size, _lhs, lhsStride, rhs);
41  }
42 };
43 
44 // forward and backward substitution, row-major, rhs is a vector
45 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
46 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
47 {
48  enum {
49  IsLower = ((Mode&Lower)==Lower)
50  };
51  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
52  {
53  typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
54  const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
55  typename internal::conditional<
56  Conjugate,
57  const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
58  const LhsMap&>
59  ::type cjLhs(lhs);
60  static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
61  for(Index pi=IsLower ? 0 : size;
62  IsLower ? pi<size : pi>0;
63  IsLower ? pi+=PanelWidth : pi-=PanelWidth)
64  {
65  Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
66 
67  Index r = IsLower ? pi : size - pi; // remaining size
68  if (r > 0)
69  {
70  // let's directly call the low level product function because:
71  // 1 - it is faster to compile
72  // 2 - it is slighlty faster at runtime
73  Index startRow = IsLower ? pi : pi-actualPanelWidth;
74  Index startCol = IsLower ? 0 : pi;
75 
77  actualPanelWidth, r,
78  &lhs.coeffRef(startRow,startCol), lhsStride,
79  rhs + startCol, 1,
80  rhs + startRow, 1,
81  RhsScalar(-1));
82  }
83 
84  for(Index k=0; k<actualPanelWidth; ++k)
85  {
86  Index i = IsLower ? pi+k : pi-k-1;
87  Index s = IsLower ? pi : i+1;
88  if (k>0)
89  rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum();
90 
91  if(!(Mode & UnitDiag))
92  rhs[i] /= cjLhs(i,i);
93  }
94  }
95  }
96 };
97 
98 // forward and backward substitution, column-major, rhs is a vector
99 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
100 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
101 {
102  enum {
103  IsLower = ((Mode&Lower)==Lower)
104  };
105  static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
106  {
107  typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
108  const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
109  typename internal::conditional<Conjugate,
110  const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
111  const LhsMap&
112  >::type cjLhs(lhs);
113  static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
114 
115  for(Index pi=IsLower ? 0 : size;
116  IsLower ? pi<size : pi>0;
117  IsLower ? pi+=PanelWidth : pi-=PanelWidth)
118  {
119  Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
120  Index startBlock = IsLower ? pi : pi-actualPanelWidth;
121  Index endBlock = IsLower ? pi + actualPanelWidth : 0;
122 
123  for(Index k=0; k<actualPanelWidth; ++k)
124  {
125  Index i = IsLower ? pi+k : pi-k-1;
126  if(!(Mode & UnitDiag))
127  rhs[i] /= cjLhs.coeff(i,i);
128 
129  Index r = actualPanelWidth - k - 1; // remaining size
130  Index s = IsLower ? i+1 : i-r;
131  if (r>0)
132  Map<Matrix<RhsScalar,Dynamic,1> >(rhs+s,r) -= rhs[i] * cjLhs.col(i).segment(s,r);
133  }
134  Index r = IsLower ? size - endBlock : startBlock; // remaining size
135  if (r > 0)
136  {
137  // let's directly call the low level product function because:
138  // 1 - it is faster to compile
139  // 2 - it is slighlty faster at runtime
141  r, actualPanelWidth,
142  &lhs.coeffRef(endBlock,startBlock), lhsStride,
143  rhs+startBlock, 1,
144  rhs+endBlock, 1, RhsScalar(-1));
145  }
146  }
147  }
148 };
149 
150 } // end namespace internal
151 
152 } // end namespace Eigen
153 
154 #endif // EIGEN_TRIANGULAR_SOLVER_VECTOR_H