blob: ce2d74daa95df3c760870425e20306b6f869273a [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2011 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_BAND_TRIANGULARSOLVER_H #define EIGEN_BAND_TRIANGULARSOLVER_H namespace internal { /* \internal * Solve Ax=b with A a band triangular matrix * TODO: extend it to matrices for x abd b */ template struct band_solve_triangular_selector; template struct band_solve_triangular_selector { typedef Map, 0, OuterStride<> > LhsMap; typedef Map > RhsMap; enum { IsLower = (Mode&Lower) ? 1 : 0 }; static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) { const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride)); RhsMap other(_other,size,1); typename internal::conditional< ConjLhs, const CwiseUnaryOp,LhsMap>, const LhsMap&> ::type cjLhs(lhs); for(int col=0 ; col0) other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose() .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum(); if((Mode&UnitDiag)==0) other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0); } } } }; template struct band_solve_triangular_selector { typedef Map, 0, OuterStride<> > LhsMap; typedef Map > RhsMap; enum { IsLower = (Mode&Lower) ? 1 : 0 }; static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) { const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride)); RhsMap other(_other,size,1); typename internal::conditional< ConjLhs, const CwiseUnaryOp,LhsMap>, const LhsMap&> ::type cjLhs(lhs); for(int col=0 ; col0) other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k) -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k); } } } }; } // end namespace internal #endif // EIGEN_BAND_TRIANGULARSOLVER_H