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