• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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
140