1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_BAND_TRIANGULARSOLVER_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_BAND_TRIANGULARSOLVER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* \internal 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Solve Ax=b with A a band triangular matrix 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * TODO: extend it to matrices for x abd b */ 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder> 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct band_solve_triangular_selector; 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar> 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,RowMajor> 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap; 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap; 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { IsLower = (Mode&Lower) ? 1 : 0 }; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride)); 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RhsMap other(_other,size,1); 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::conditional< 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ConjLhs, 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>, 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const LhsMap&> 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::type cjLhs(lhs); 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int ii=0; ii<size; ++ii) 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i = IsLower ? ii : size-ii-1; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int actual_k = (std::min)(k,ii); 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int actual_start = IsLower ? k-actual_k : 1; 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(actual_k>0) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose() 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum(); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((Mode&UnitDiag)==0) 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar> 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ColMajor> 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { IsLower = (Mode&Lower) ? 1 : 0 }; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride)); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RhsMap other(_other,size,1); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::conditional< 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ConjLhs, 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>, 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const LhsMap&> 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::type cjLhs(lhs); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int ii=0; ii<size; ++ii) 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i = IsLower ? ii : size-ii-1; 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int actual_k = (std::min)(k,size-ii-1); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int actual_start = IsLower ? 1 : k-actual_k; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((Mode&UnitDiag)==0) 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) /= cjLhs(IsLower ? 0 : k, i); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(actual_k>0) 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k) 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_BAND_TRIANGULARSOLVER_H 98