[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2008-2010 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_TRIANGULAR_SOLVER_VECTOR_H
|
---|
| 11 | #define EIGEN_TRIANGULAR_SOLVER_VECTOR_H
|
---|
| 12 |
|
---|
| 13 | namespace Eigen {
|
---|
| 14 |
|
---|
| 15 | namespace internal {
|
---|
| 16 |
|
---|
| 17 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
|
---|
| 18 | struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
|
---|
| 19 | {
|
---|
| 20 | static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
|
---|
| 21 | {
|
---|
| 22 | triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
|
---|
| 23 | ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
|
---|
| 24 | Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
|
---|
| 25 | >::run(size, _lhs, lhsStride, rhs);
|
---|
| 26 | }
|
---|
| 27 | };
|
---|
| 28 |
|
---|
| 29 | // forward and backward substitution, row-major, rhs is a vector
|
---|
| 30 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
|
---|
| 31 | struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
|
---|
| 32 | {
|
---|
| 33 | enum {
|
---|
| 34 | IsLower = ((Mode&Lower)==Lower)
|
---|
| 35 | };
|
---|
| 36 | static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
|
---|
| 37 | {
|
---|
| 38 | typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
|
---|
| 39 | const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
|
---|
| 40 | typename internal::conditional<
|
---|
| 41 | Conjugate,
|
---|
| 42 | const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
|
---|
| 43 | const LhsMap&>
|
---|
| 44 | ::type cjLhs(lhs);
|
---|
| 45 | static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
|
---|
| 46 | for(Index pi=IsLower ? 0 : size;
|
---|
| 47 | IsLower ? pi<size : pi>0;
|
---|
| 48 | IsLower ? pi+=PanelWidth : pi-=PanelWidth)
|
---|
| 49 | {
|
---|
| 50 | Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
|
---|
| 51 |
|
---|
| 52 | Index r = IsLower ? pi : size - pi; // remaining size
|
---|
| 53 | if (r > 0)
|
---|
| 54 | {
|
---|
| 55 | // let's directly call the low level product function because:
|
---|
| 56 | // 1 - it is faster to compile
|
---|
| 57 | // 2 - it is slighlty faster at runtime
|
---|
| 58 | Index startRow = IsLower ? pi : pi-actualPanelWidth;
|
---|
| 59 | Index startCol = IsLower ? 0 : pi;
|
---|
| 60 |
|
---|
| 61 | general_matrix_vector_product<Index,LhsScalar,RowMajor,Conjugate,RhsScalar,false>::run(
|
---|
| 62 | actualPanelWidth, r,
|
---|
| 63 | &lhs.coeffRef(startRow,startCol), lhsStride,
|
---|
| 64 | rhs + startCol, 1,
|
---|
| 65 | rhs + startRow, 1,
|
---|
| 66 | RhsScalar(-1));
|
---|
| 67 | }
|
---|
| 68 |
|
---|
| 69 | for(Index k=0; k<actualPanelWidth; ++k)
|
---|
| 70 | {
|
---|
| 71 | Index i = IsLower ? pi+k : pi-k-1;
|
---|
| 72 | Index s = IsLower ? pi : i+1;
|
---|
| 73 | if (k>0)
|
---|
| 74 | rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum();
|
---|
| 75 |
|
---|
| 76 | if(!(Mode & UnitDiag))
|
---|
| 77 | rhs[i] /= cjLhs(i,i);
|
---|
| 78 | }
|
---|
| 79 | }
|
---|
| 80 | }
|
---|
| 81 | };
|
---|
| 82 |
|
---|
| 83 | // forward and backward substitution, column-major, rhs is a vector
|
---|
| 84 | template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
|
---|
| 85 | struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
|
---|
| 86 | {
|
---|
| 87 | enum {
|
---|
| 88 | IsLower = ((Mode&Lower)==Lower)
|
---|
| 89 | };
|
---|
| 90 | static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
|
---|
| 91 | {
|
---|
| 92 | typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
|
---|
| 93 | const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
|
---|
| 94 | typename internal::conditional<Conjugate,
|
---|
| 95 | const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
|
---|
| 96 | const LhsMap&
|
---|
| 97 | >::type cjLhs(lhs);
|
---|
| 98 | static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
|
---|
| 99 |
|
---|
| 100 | for(Index pi=IsLower ? 0 : size;
|
---|
| 101 | IsLower ? pi<size : pi>0;
|
---|
| 102 | IsLower ? pi+=PanelWidth : pi-=PanelWidth)
|
---|
| 103 | {
|
---|
| 104 | Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
|
---|
| 105 | Index startBlock = IsLower ? pi : pi-actualPanelWidth;
|
---|
| 106 | Index endBlock = IsLower ? pi + actualPanelWidth : 0;
|
---|
| 107 |
|
---|
| 108 | for(Index k=0; k<actualPanelWidth; ++k)
|
---|
| 109 | {
|
---|
| 110 | Index i = IsLower ? pi+k : pi-k-1;
|
---|
| 111 | if(!(Mode & UnitDiag))
|
---|
| 112 | rhs[i] /= cjLhs.coeff(i,i);
|
---|
| 113 |
|
---|
| 114 | Index r = actualPanelWidth - k - 1; // remaining size
|
---|
| 115 | Index s = IsLower ? i+1 : i-r;
|
---|
| 116 | if (r>0)
|
---|
| 117 | Map<Matrix<RhsScalar,Dynamic,1> >(rhs+s,r) -= rhs[i] * cjLhs.col(i).segment(s,r);
|
---|
| 118 | }
|
---|
| 119 | Index r = IsLower ? size - endBlock : startBlock; // remaining size
|
---|
| 120 | if (r > 0)
|
---|
| 121 | {
|
---|
| 122 | // let's directly call the low level product function because:
|
---|
| 123 | // 1 - it is faster to compile
|
---|
| 124 | // 2 - it is slighlty faster at runtime
|
---|
| 125 | general_matrix_vector_product<Index,LhsScalar,ColMajor,Conjugate,RhsScalar,false>::run(
|
---|
| 126 | r, actualPanelWidth,
|
---|
| 127 | &lhs.coeffRef(endBlock,startBlock), lhsStride,
|
---|
| 128 | rhs+startBlock, 1,
|
---|
| 129 | rhs+endBlock, 1, RhsScalar(-1));
|
---|
| 130 | }
|
---|
| 131 | }
|
---|
| 132 | }
|
---|
| 133 | };
|
---|
| 134 |
|
---|
| 135 | } // end namespace internal
|
---|
| 136 |
|
---|
| 137 | } // end namespace Eigen
|
---|
| 138 |
|
---|
| 139 | #endif // EIGEN_TRIANGULAR_SOLVER_VECTOR_H
|
---|