TriangularSolverVector.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
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
13namespace Eigen {
14
15namespace internal {
16
17template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
18struct 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
30template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
31struct 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
84template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
85struct 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