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 | {
|
---|
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
|
---|