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> 19 struct band_solve_triangular_selector; 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 98