PackedTriangularSolverVector.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) 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_SOLVER_VECTOR_H
11 #define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
12 
13 namespace internal {
14 
15 template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
17 
18 // forward and backward substitution, row-major, rhs is a vector
19 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
20 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
21 {
22  enum {
23  IsLower = (Mode&Lower)==Lower
24  };
25  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
26  {
27  internal::conj_if<Conjugate> cj;
28  typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
29  typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
30 
31  lhs += IsLower ? 0 : (size*(size+1)>>1)-1;
32  for(Index pi=0; pi<size; ++pi)
33  {
34  Index i = IsLower ? pi : size-pi-1;
35  Index s = IsLower ? 0 : 1;
36  if (pi>0)
37  rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi))
38  .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum();
39  if (!(Mode & UnitDiag))
40  rhs[i] /= cj(lhs[IsLower ? i : 0]);
41  IsLower ? lhs += pi+1 : lhs -= pi+2;
42  }
43  }
44 };
45 
46 // forward and backward substitution, column-major, rhs is a vector
47 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
48 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
49 {
50  enum {
51  IsLower = (Mode&Lower)==Lower
52  };
53  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
54  {
55  internal::conj_if<Conjugate> cj;
56  typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
57  typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
58 
59  lhs += IsLower ? 0 : size*(size-1)>>1;
60  for(Index pi=0; pi<size; ++pi)
61  {
62  Index i = IsLower ? pi : size-pi-1;
63  Index r = size - pi - 1;
64  if (!(Mode & UnitDiag))
65  rhs[i] /= cj(lhs[IsLower ? 0 : i]);
66  if (r>0)
67  Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -=
68  rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r));
69  IsLower ? lhs += size-pi : lhs -= r;
70  }
71  }
72 };
73 
74 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
75 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
76 {
77  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
78  {
80  ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
81  Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
82  >::run(size, lhs, rhs);
83  }
84 };
85 
86 } // end namespace internal
87 
88 #endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
RealScalar s


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:59