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
13namespace internal {
14
15 /* \internal
16  * Solve Ax=b with A a band triangular matrix
17  * TODO: extend it to matrices for x abd b */
18template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
19struct band_solve_triangular_selector;
20
21
22template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
23struct 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
58template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
59struct 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
98