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_TRIANGULARMATRIXVECTOR_H
00026 #define EIGEN_TRIANGULARMATRIXVECTOR_H
00027
00028 namespace internal {
00029
00030 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder>
00031 struct product_triangular_matrix_vector;
00032
00033 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
00034 struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor>
00035 {
00036 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00037 enum {
00038 IsLower = ((Mode&Lower)==Lower),
00039 HasUnitDiag = (Mode & UnitDiag)==UnitDiag
00040 };
00041 static EIGEN_DONT_INLINE void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
00042 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
00043 {
00044 EIGEN_UNUSED_VARIABLE(resIncr);
00045 eigen_assert(resIncr==1);
00046
00047 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00048
00049 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
00050 const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
00051 typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
00052
00053 typedef Map<const Matrix<RhsScalar,Dynamic,1>, 0, InnerStride<> > RhsMap;
00054 const RhsMap rhs(_rhs,cols,InnerStride<>(rhsIncr));
00055 typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
00056
00057 typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap;
00058 ResMap res(_res,rows);
00059
00060 for (Index pi=0; pi<cols; pi+=PanelWidth)
00061 {
00062 Index actualPanelWidth = std::min(PanelWidth, cols-pi);
00063 for (Index k=0; k<actualPanelWidth; ++k)
00064 {
00065 Index i = pi + k;
00066 Index s = IsLower ? (HasUnitDiag ? i+1 : i ) : pi;
00067 Index r = IsLower ? actualPanelWidth-k : k+1;
00068 if ((!HasUnitDiag) || (--r)>0)
00069 res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
00070 if (HasUnitDiag)
00071 res.coeffRef(i) += alpha * cjRhs.coeff(i);
00072 }
00073 Index r = IsLower ? cols - pi - actualPanelWidth : pi;
00074 if (r>0)
00075 {
00076 Index s = IsLower ? pi+actualPanelWidth : 0;
00077 general_matrix_vector_product<Index,LhsScalar,ColMajor,ConjLhs,RhsScalar,ConjRhs>::run(
00078 r, actualPanelWidth,
00079 &lhs.coeffRef(s,pi), lhsStride,
00080 &rhs.coeffRef(pi), rhsIncr,
00081 &res.coeffRef(s), resIncr, alpha);
00082 }
00083 }
00084 }
00085 };
00086
00087 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
00088 struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor>
00089 {
00090 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00091 enum {
00092 IsLower = ((Mode&Lower)==Lower),
00093 HasUnitDiag = (Mode & UnitDiag)==UnitDiag
00094 };
00095 static void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
00096 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
00097 {
00098 eigen_assert(rhsIncr==1);
00099 EIGEN_UNUSED_VARIABLE(rhsIncr);
00100
00101 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00102
00103 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
00104 const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
00105 typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
00106
00107 typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap;
00108 const RhsMap rhs(_rhs,cols);
00109 typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
00110
00111 typedef Map<Matrix<ResScalar,Dynamic,1>, 0, InnerStride<> > ResMap;
00112 ResMap res(_res,rows,InnerStride<>(resIncr));
00113
00114 for (Index pi=0; pi<cols; pi+=PanelWidth)
00115 {
00116 Index actualPanelWidth = std::min(PanelWidth, cols-pi);
00117 for (Index k=0; k<actualPanelWidth; ++k)
00118 {
00119 Index i = pi + k;
00120 Index s = IsLower ? pi : (HasUnitDiag ? i+1 : i);
00121 Index r = IsLower ? k+1 : actualPanelWidth-k;
00122 if ((!HasUnitDiag) || (--r)>0)
00123 res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum();
00124 if (HasUnitDiag)
00125 res.coeffRef(i) += alpha * cjRhs.coeff(i);
00126 }
00127 Index r = IsLower ? pi : cols - pi - actualPanelWidth;
00128 if (r>0)
00129 {
00130 Index s = IsLower ? 0 : pi + actualPanelWidth;
00131 general_matrix_vector_product<Index,LhsScalar,RowMajor,ConjLhs,RhsScalar,ConjRhs>::run(
00132 actualPanelWidth, r,
00133 &lhs.coeffRef(pi,s), lhsStride,
00134 &rhs.coeffRef(s), rhsIncr,
00135 &res.coeffRef(pi), resIncr, alpha);
00136 }
00137 }
00138 }
00139 };
00140
00141
00142
00143
00144
00145 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
00146 struct traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true> >
00147 : traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true>, Lhs, Rhs> >
00148 {};
00149
00150 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
00151 struct traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false> >
00152 : traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false>, Lhs, Rhs> >
00153 {};
00154
00155 }
00156
00157 template<int Mode, typename Lhs, typename Rhs>
00158 struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
00159 : public ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs >
00160 {
00161 EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00162
00163 TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00164
00165 template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00166 {
00167 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00168
00169 const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
00170 const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
00171
00172 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
00173 * RhsBlasTraits::extractScalarFactor(m_rhs);
00174
00175 internal::product_triangular_matrix_vector
00176 <Index,Mode,
00177 typename _ActualLhsType::Scalar, LhsBlasTraits::NeedToConjugate,
00178 typename _ActualRhsType::Scalar, RhsBlasTraits::NeedToConjugate,
00179 (int(internal::traits<Lhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor>
00180 ::run(lhs.rows(),lhs.cols(),lhs.data(),lhs.outerStride(),rhs.data(),rhs.innerStride(),
00181 dst.data(),dst.innerStride(),actualAlpha);
00182 }
00183 };
00184
00185 template<int Mode, typename Lhs, typename Rhs>
00186 struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
00187 : public ProductBase<TriangularProduct<Mode,false,Lhs,true,Rhs,false>, Lhs, Rhs >
00188 {
00189 EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00190
00191 TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00192
00193 template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00194 {
00195
00196 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00197
00198 const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
00199 const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
00200
00201 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
00202 * RhsBlasTraits::extractScalarFactor(m_rhs);
00203
00204 internal::product_triangular_matrix_vector
00205 <Index,(Mode & UnitDiag) | (Mode & Lower) ? Upper : Lower,
00206 typename _ActualRhsType::Scalar, RhsBlasTraits::NeedToConjugate,
00207 typename _ActualLhsType::Scalar, LhsBlasTraits::NeedToConjugate,
00208 (int(internal::traits<Rhs>::Flags)&RowMajorBit) ? ColMajor : RowMajor>
00209 ::run(rhs.rows(),rhs.cols(),rhs.data(),rhs.outerStride(),lhs.data(),lhs.innerStride(),
00210 dst.data(),dst.innerStride(),actualAlpha);
00211 }
00212 };
00213
00214 #endif // EIGEN_TRIANGULARMATRIXVECTOR_H