BandTriangularSolver.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) 2011 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_BAND_TRIANGULARSOLVER_H
11 #define EIGEN_BAND_TRIANGULARSOLVER_H
12 
13 namespace internal {
14 
15  /* \internal
16  * Solve Ax=b with A a band triangular matrix
17  * TODO: extend it to matrices for x abd b */
18 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
20 
21 
22 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
23 struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,RowMajor>
24 {
25  typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
26  typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
27  enum { IsLower = (Mode&Lower) ? 1 : 0 };
28  static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
29  {
30  const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride));
31  RhsMap other(_other,size,1);
32  typename internal::conditional<
33  ConjLhs,
34  const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
35  const LhsMap&>
36  ::type cjLhs(lhs);
37 
38  for(int col=0 ; col<other.cols() ; ++col)
39  {
40  for(int ii=0; ii<size; ++ii)
41  {
42  int i = IsLower ? ii : size-ii-1;
43  int actual_k = (std::min)(k,ii);
44  int actual_start = IsLower ? k-actual_k : 1;
45 
46  if(actual_k>0)
47  other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose()
48  .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum();
49 
50  if((Mode&UnitDiag)==0)
51  other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0);
52  }
53  }
54  }
55 
56 };
57 
58 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
59 struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ColMajor>
60 {
61  typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
62  typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
63  enum { IsLower = (Mode&Lower) ? 1 : 0 };
64  static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
65  {
66  const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride));
67  RhsMap other(_other,size,1);
68  typename internal::conditional<
69  ConjLhs,
70  const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
71  const LhsMap&>
72  ::type cjLhs(lhs);
73 
74  for(int col=0 ; col<other.cols() ; ++col)
75  {
76  for(int ii=0; ii<size; ++ii)
77  {
78  int i = IsLower ? ii : size-ii-1;
79  int actual_k = (std::min)(k,size-ii-1);
80  int actual_start = IsLower ? 1 : k-actual_k;
81 
82  if((Mode&UnitDiag)==0)
83  other.coeffRef(i,col) /= cjLhs(IsLower ? 0 : k, i);
84 
85  if(actual_k>0)
86  other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k)
87  -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k);
88 
89  }
90  }
91  }
92 };
93 
94 
95 } // end namespace internal
96 
97 #endif // EIGEN_BAND_TRIANGULARSOLVER_H
col
m col(1)
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, RowMajor >::LhsMap
Map< const Matrix< LhsScalar, Dynamic, Dynamic, RowMajor >, 0, OuterStride<> > LhsMap
Definition: BandTriangularSolver.h:25
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, ColMajor >::LhsMap
Map< const Matrix< LhsScalar, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > LhsMap
Definition: BandTriangularSolver.h:61
internal::band_solve_triangular_selector
Definition: BandTriangularSolver.h:19
type
Definition: pytypes.h:1525
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, RowMajor >::run
static void run(Index size, Index k, const LhsScalar *_lhs, Index lhsStride, RhsScalar *_other)
Definition: BandTriangularSolver.h:28
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, ColMajor >::RhsMap
Map< Matrix< RhsScalar, Dynamic, 1 > > RhsMap
Definition: BandTriangularSolver.h:62
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, ColMajor >::run
static void run(Index size, Index k, const LhsScalar *_lhs, Index lhsStride, RhsScalar *_other)
Definition: BandTriangularSolver.h:64
Eigen::Lower
@ Lower
Definition: Constants.h:209
internal::band_solve_triangular_selector< Index, Mode, LhsScalar, ConjLhs, RhsScalar, RowMajor >::RhsMap
Map< Matrix< RhsScalar, Dynamic, 1 > > RhsMap
Definition: BandTriangularSolver.h:26
min
#define min(a, b)
Definition: datatypes.h:19
internal
Definition: BandTriangularSolver.h:13
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:319
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
pybind_wrapper_test_script.other
other
Definition: pybind_wrapper_test_script.py:42
Eigen::UnitDiag
@ UnitDiag
Definition: Constants.h:213
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:31:59