00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_TRIANGULAR_SOLVER_VECTOR_H
00026 #define EIGEN_TRIANGULAR_SOLVER_VECTOR_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
00033 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
00034 {
00035 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00036 {
00037 triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
00038 ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
00039 Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
00040 >::run(size, _lhs, lhsStride, rhs);
00041 }
00042 };
00043
00044
00045 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
00046 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
00047 {
00048 enum {
00049 IsLower = ((Mode&Lower)==Lower)
00050 };
00051 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00052 {
00053 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
00054 const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
00055 typename internal::conditional<
00056 Conjugate,
00057 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
00058 const LhsMap&>
00059 ::type cjLhs(lhs);
00060 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00061 for(Index pi=IsLower ? 0 : size;
00062 IsLower ? pi<size : pi>0;
00063 IsLower ? pi+=PanelWidth : pi-=PanelWidth)
00064 {
00065 Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
00066
00067 Index r = IsLower ? pi : size - pi;
00068 if (r > 0)
00069 {
00070
00071
00072
00073 Index startRow = IsLower ? pi : pi-actualPanelWidth;
00074 Index startCol = IsLower ? 0 : pi;
00075
00076 general_matrix_vector_product<Index,LhsScalar,RowMajor,Conjugate,RhsScalar,false>::run(
00077 actualPanelWidth, r,
00078 &lhs.coeffRef(startRow,startCol), lhsStride,
00079 rhs + startCol, 1,
00080 rhs + startRow, 1,
00081 RhsScalar(-1));
00082 }
00083
00084 for(Index k=0; k<actualPanelWidth; ++k)
00085 {
00086 Index i = IsLower ? pi+k : pi-k-1;
00087 Index s = IsLower ? pi : i+1;
00088 if (k>0)
00089 rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum();
00090
00091 if(!(Mode & UnitDiag))
00092 rhs[i] /= cjLhs(i,i);
00093 }
00094 }
00095 }
00096 };
00097
00098
00099 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
00100 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
00101 {
00102 enum {
00103 IsLower = ((Mode&Lower)==Lower)
00104 };
00105 static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00106 {
00107 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
00108 const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
00109 typename internal::conditional<Conjugate,
00110 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
00111 const LhsMap&
00112 >::type cjLhs(lhs);
00113 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00114
00115 for(Index pi=IsLower ? 0 : size;
00116 IsLower ? pi<size : pi>0;
00117 IsLower ? pi+=PanelWidth : pi-=PanelWidth)
00118 {
00119 Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
00120 Index startBlock = IsLower ? pi : pi-actualPanelWidth;
00121 Index endBlock = IsLower ? pi + actualPanelWidth : 0;
00122
00123 for(Index k=0; k<actualPanelWidth; ++k)
00124 {
00125 Index i = IsLower ? pi+k : pi-k-1;
00126 if(!(Mode & UnitDiag))
00127 rhs[i] /= cjLhs.coeff(i,i);
00128
00129 Index r = actualPanelWidth - k - 1;
00130 Index s = IsLower ? i+1 : i-r;
00131 if (r>0)
00132 Map<Matrix<RhsScalar,Dynamic,1> >(rhs+s,r) -= rhs[i] * cjLhs.col(i).segment(s,r);
00133 }
00134 Index r = IsLower ? size - endBlock : startBlock;
00135 if (r > 0)
00136 {
00137
00138
00139
00140 general_matrix_vector_product<Index,LhsScalar,ColMajor,Conjugate,RhsScalar,false>::run(
00141 r, actualPanelWidth,
00142 &lhs.coeffRef(endBlock,startBlock), lhsStride,
00143 rhs+startBlock, 1,
00144 rhs+endBlock, 1, RhsScalar(-1));
00145 }
00146 }
00147 }
00148 };
00149
00150 }
00151
00152 }
00153
00154 #endif // EIGEN_TRIANGULAR_SOLVER_VECTOR_H