[136] | 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
|
---|