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 // 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  {
39  const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
40  typename internal::conditional<
41  Conjugate,
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 
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  {
93  const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
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
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
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
static void run(Index size, const LhsScalar *_lhs, Index lhsStride, RhsScalar *rhs)
#define EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH
Definition: Settings.h:38
void rhs(const real_t *x, real_t *f)
static void run(Index size, const LhsScalar *_lhs, Index lhsStride, RhsScalar *rhs)
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:59
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
Convenience specialization of Stride to specify only an outer stride See class Map for some examples...
Definition: Stride.h:97


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:14