10 #ifndef EIGEN_TRIANGULARMATRIXVECTOR_H 11 #define EIGEN_TRIANGULARMATRIXVECTOR_H 17 template<
typename Index,
int Mode,
typename LhsScalar,
bool ConjLhs,
typename RhsScalar,
bool ConjRhs,
int StorageOrder,
int Version=Specialized>
18 struct triangular_matrix_vector_product;
20 template<
typename Index,
int Mode,
typename LhsScalar,
bool ConjLhs,
typename RhsScalar,
bool ConjRhs,
int Version>
21 struct triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,
ColMajor,Version>
23 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
27 HasZeroDiag = (Mode &
ZeroDiag)==ZeroDiag
29 static EIGEN_DONT_INLINE
void run(Index _rows, Index _cols,
const LhsScalar* _lhs, Index lhsStride,
30 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr,
const ResScalar& alpha);
33 template<
typename Index,
int Mode,
typename LhsScalar,
bool ConjLhs,
typename RhsScalar,
bool ConjRhs,
int Version>
34 EIGEN_DONT_INLINE
void triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor,Version>
35 ::run(Index _rows, Index _cols,
const LhsScalar* _lhs, Index lhsStride,
36 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr,
const ResScalar& alpha)
38 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
39 Index size = (std::min)(_rows,_cols);
40 Index rows = IsLower ? _rows : (std::min)(_rows,_cols);
41 Index cols = IsLower ? (std::min)(_rows,_cols) : _cols;
43 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
44 const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
45 typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
47 typedef Map<const Matrix<RhsScalar,Dynamic,1>, 0, InnerStride<> > RhsMap;
48 const RhsMap rhs(_rhs,cols,InnerStride<>(rhsIncr));
49 typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
51 typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap;
52 ResMap res(_res,rows);
54 typedef const_blas_data_mapper<LhsScalar,Index,ColMajor> LhsMapper;
55 typedef const_blas_data_mapper<RhsScalar,Index,RowMajor> RhsMapper;
57 for (Index pi=0; pi<size; pi+=PanelWidth)
59 Index actualPanelWidth = (std::min)(PanelWidth, size-pi);
60 for (Index k=0; k<actualPanelWidth; ++k)
63 Index s = IsLower ? ((HasUnitDiag||HasZeroDiag) ? i+1 : i ) : pi;
64 Index r = IsLower ? actualPanelWidth-k : k+1;
65 if ((!(HasUnitDiag||HasZeroDiag)) || (--r)>0)
66 res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
68 res.coeffRef(i) += alpha * cjRhs.coeff(i);
70 Index r = IsLower ? rows - pi - actualPanelWidth : pi;
73 Index s = IsLower ? pi+actualPanelWidth : 0;
74 general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjLhs,RhsScalar,RhsMapper,ConjRhs,BuiltIn>::run(
76 LhsMapper(&lhs.coeffRef(s,pi), lhsStride),
77 RhsMapper(&rhs.coeffRef(pi), rhsIncr),
78 &res.coeffRef(s), resIncr, alpha);
81 if((!IsLower) && cols>size)
83 general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjLhs,RhsScalar,RhsMapper,ConjRhs>::run(
85 LhsMapper(&lhs.coeffRef(0,size), lhsStride),
86 RhsMapper(&rhs.coeffRef(size), rhsIncr),
87 _res, resIncr, alpha);
91 template<
typename Index,
int Mode,
typename LhsScalar,
bool ConjLhs,
typename RhsScalar,
bool ConjRhs,
int Version>
92 struct triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,
RowMajor,Version>
94 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
98 HasZeroDiag = (Mode &
ZeroDiag)==ZeroDiag
100 static EIGEN_DONT_INLINE
void run(Index _rows, Index _cols,
const LhsScalar* _lhs, Index lhsStride,
101 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr,
const ResScalar& alpha);
104 template<
typename Index,
int Mode,
typename LhsScalar,
bool ConjLhs,
typename RhsScalar,
bool ConjRhs,
int Version>
105 EIGEN_DONT_INLINE
void triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor,Version>
106 ::run(Index _rows, Index _cols,
const LhsScalar* _lhs, Index lhsStride,
107 const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr,
const ResScalar& alpha)
109 static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
110 Index diagSize = (std::min)(_rows,_cols);
111 Index rows = IsLower ? _rows : diagSize;
112 Index cols = IsLower ? diagSize : _cols;
114 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
115 const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
116 typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
118 typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap;
119 const RhsMap rhs(_rhs,cols);
120 typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
122 typedef Map<Matrix<ResScalar,Dynamic,1>, 0, InnerStride<> > ResMap;
123 ResMap res(_res,rows,InnerStride<>(resIncr));
125 typedef const_blas_data_mapper<LhsScalar,Index,RowMajor> LhsMapper;
126 typedef const_blas_data_mapper<RhsScalar,Index,RowMajor> RhsMapper;
128 for (Index pi=0; pi<diagSize; pi+=PanelWidth)
130 Index actualPanelWidth = (std::min)(PanelWidth, diagSize-pi);
131 for (Index k=0; k<actualPanelWidth; ++k)
134 Index s = IsLower ? pi : ((HasUnitDiag||HasZeroDiag) ? i+1 : i);
135 Index r = IsLower ? k+1 : actualPanelWidth-k;
136 if ((!(HasUnitDiag||HasZeroDiag)) || (--r)>0)
137 res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum();
139 res.coeffRef(i) += alpha * cjRhs.coeff(i);
141 Index r = IsLower ? pi : cols - pi - actualPanelWidth;
144 Index s = IsLower ? 0 : pi + actualPanelWidth;
145 general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjLhs,RhsScalar,RhsMapper,ConjRhs,BuiltIn>::run(
147 LhsMapper(&lhs.coeffRef(pi,s), lhsStride),
148 RhsMapper(&rhs.coeffRef(s), rhsIncr),
149 &res.coeffRef(pi), resIncr, alpha);
152 if(IsLower && rows>diagSize)
154 general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjLhs,RhsScalar,RhsMapper,ConjRhs>::run(
156 LhsMapper(&lhs.coeffRef(diagSize,0), lhsStride),
157 RhsMapper(&rhs.coeffRef(0), rhsIncr),
158 &res.coeffRef(diagSize), resIncr, alpha);
166 template<
int Mode,
int StorageOrder>
167 struct trmv_selector;
173 template<
int Mode,
typename Lhs,
typename Rhs>
174 struct triangular_product_impl<Mode,true,Lhs,false,Rhs,true>
176 template<
typename Dest>
static void run(Dest& dst,
const Lhs &lhs,
const Rhs &rhs,
const typename Dest::Scalar& alpha)
178 eigen_assert(dst.rows()==lhs.rows() && dst.cols()==rhs.cols());
184 template<
int Mode,
typename Lhs,
typename Rhs>
185 struct triangular_product_impl<Mode,false,Lhs,true,Rhs,false>
187 template<
typename Dest>
static void run(Dest& dst,
const Lhs &lhs,
const Rhs &rhs,
const typename Dest::Scalar& alpha)
189 eigen_assert(dst.rows()==lhs.rows() && dst.cols()==rhs.cols());
191 Transpose<Dest> dstT(dst);
194 ::run(rhs.transpose(),lhs.transpose(), dstT, alpha);
204 template<
int Mode>
struct trmv_selector<Mode,
ColMajor>
206 template<
typename Lhs,
typename Rhs,
typename Dest>
207 static void run(
const Lhs &lhs,
const Rhs &rhs, Dest& dest,
const typename Dest::Scalar& alpha)
209 typedef typename Lhs::Scalar LhsScalar;
210 typedef typename Rhs::Scalar RhsScalar;
211 typedef typename Dest::Scalar ResScalar;
212 typedef typename Dest::RealScalar RealScalar;
214 typedef internal::blas_traits<Lhs> LhsBlasTraits;
215 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
216 typedef internal::blas_traits<Rhs> RhsBlasTraits;
217 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
219 typedef Map<Matrix<ResScalar,Dynamic,1>,
Aligned> MappedDest;
221 typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(lhs);
222 typename internal::add_const_on_value_type<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(rhs);
224 ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs)
225 * RhsBlasTraits::extractScalarFactor(rhs);
230 EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
231 ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
232 MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
235 gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
237 bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
238 bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
240 RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
242 ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
243 evalToDest ? dest.data() : static_dest.data());
247 #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN 248 Index size = dest.size();
249 EIGEN_DENSE_STORAGE_CTOR_PLUGIN
251 if(!alphaIsCompatible)
253 MappedDest(actualDestPtr, dest.size()).setZero();
254 compatibleAlpha = RhsScalar(1);
257 MappedDest(actualDestPtr, dest.size()) = dest;
260 internal::triangular_matrix_vector_product
262 LhsScalar, LhsBlasTraits::NeedToConjugate,
263 RhsScalar, RhsBlasTraits::NeedToConjugate,
265 ::run(actualLhs.rows(),actualLhs.cols(),
266 actualLhs.data(),actualLhs.outerStride(),
267 actualRhs.data(),actualRhs.innerStride(),
268 actualDestPtr,1,compatibleAlpha);
272 if(!alphaIsCompatible)
273 dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
275 dest = MappedDest(actualDestPtr, dest.size());
280 template<
int Mode>
struct trmv_selector<Mode,
RowMajor>
282 template<
typename Lhs,
typename Rhs,
typename Dest>
283 static void run(
const Lhs &lhs,
const Rhs &rhs, Dest& dest,
const typename Dest::Scalar& alpha)
285 typedef typename Lhs::Scalar LhsScalar;
286 typedef typename Rhs::Scalar RhsScalar;
287 typedef typename Dest::Scalar ResScalar;
289 typedef internal::blas_traits<Lhs> LhsBlasTraits;
290 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
291 typedef internal::blas_traits<Rhs> RhsBlasTraits;
292 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
293 typedef typename internal::remove_all<ActualRhsType>::type ActualRhsTypeCleaned;
295 typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(lhs);
296 typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(rhs);
298 ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs)
299 * RhsBlasTraits::extractScalarFactor(rhs);
302 DirectlyUseRhs = ActualRhsTypeCleaned::InnerStrideAtCompileTime==1
305 gemv_static_vector_if<RhsScalar,ActualRhsTypeCleaned::SizeAtCompileTime,ActualRhsTypeCleaned::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
307 ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
308 DirectlyUseRhs ?
const_cast<RhsScalar*
>(actualRhs.data()) : static_rhs.data());
312 #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN 313 Index size = actualRhs.size();
314 EIGEN_DENSE_STORAGE_CTOR_PLUGIN
316 Map<typename ActualRhsTypeCleaned::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
319 internal::triangular_matrix_vector_product
321 LhsScalar, LhsBlasTraits::NeedToConjugate,
322 RhsScalar, RhsBlasTraits::NeedToConjugate,
324 ::run(actualLhs.rows(),actualLhs.cols(),
325 actualLhs.data(),actualLhs.outerStride(),
327 dest.data(),dest.innerStride(),
336 #endif // EIGEN_TRIANGULARMATRIXVECTOR_H Definition: Constants.h:210
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: Constants.h:320
Definition: Constants.h:204
Definition: Constants.h:322
Definition: Eigen_Colamd.h:54
Definition: Constants.h:208
Definition: Constants.h:235
Definition: Constants.h:206