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_SKYLINEPRODUCT_H
00026 #define EIGEN_SKYLINEPRODUCT_H
00027
00028 namespace Eigen {
00029
00030 template<typename Lhs, typename Rhs, int ProductMode>
00031 struct SkylineProductReturnType {
00032 typedef const typename internal::nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested;
00033 typedef const typename internal::nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested;
00034
00035 typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type;
00036 };
00037
00038 template<typename LhsNested, typename RhsNested, int ProductMode>
00039 struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > {
00040
00041 typedef typename internal::remove_all<LhsNested>::type _LhsNested;
00042 typedef typename internal::remove_all<RhsNested>::type _RhsNested;
00043 typedef typename _LhsNested::Scalar Scalar;
00044
00045 enum {
00046 LhsCoeffReadCost = _LhsNested::CoeffReadCost,
00047 RhsCoeffReadCost = _RhsNested::CoeffReadCost,
00048 LhsFlags = _LhsNested::Flags,
00049 RhsFlags = _RhsNested::Flags,
00050
00051 RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
00052 ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
00053 InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
00054
00055 MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
00056 MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
00057
00058 EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
00059 ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct,
00060
00061 RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)),
00062
00063 Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
00064 | EvalBeforeAssigningBit
00065 | EvalBeforeNestingBit,
00066
00067 CoeffReadCost = Dynamic
00068 };
00069
00070 typedef typename internal::conditional<ResultIsSkyline,
00071 SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >,
00072 MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::type Base;
00073 };
00074
00075 namespace internal {
00076 template<typename LhsNested, typename RhsNested, int ProductMode>
00077 class SkylineProduct : no_assignment_operator,
00078 public traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base {
00079 public:
00080
00081 EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct)
00082
00083 private:
00084
00085 typedef typename traits<SkylineProduct>::_LhsNested _LhsNested;
00086 typedef typename traits<SkylineProduct>::_RhsNested _RhsNested;
00087
00088 public:
00089
00090 template<typename Lhs, typename Rhs>
00091 EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs)
00092 : m_lhs(lhs), m_rhs(rhs) {
00093 eigen_assert(lhs.cols() == rhs.rows());
00094
00095 enum {
00096 ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic
00097 || _RhsNested::RowsAtCompileTime == Dynamic
00098 || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime),
00099 AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
00100 SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested)
00101 };
00102
00103
00104
00105 EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00106 INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00107 EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00108 INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00109 EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00110 }
00111
00112 EIGEN_STRONG_INLINE Index rows() const {
00113 return m_lhs.rows();
00114 }
00115
00116 EIGEN_STRONG_INLINE Index cols() const {
00117 return m_rhs.cols();
00118 }
00119
00120 EIGEN_STRONG_INLINE const _LhsNested& lhs() const {
00121 return m_lhs;
00122 }
00123
00124 EIGEN_STRONG_INLINE const _RhsNested& rhs() const {
00125 return m_rhs;
00126 }
00127
00128 protected:
00129 LhsNested m_lhs;
00130 RhsNested m_rhs;
00131 };
00132
00133
00134
00135
00136 template<typename Lhs, typename Rhs, typename Dest>
00137 EIGEN_DONT_INLINE void skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
00138 typedef typename remove_all<Lhs>::type _Lhs;
00139 typedef typename remove_all<Rhs>::type _Rhs;
00140 typedef typename traits<Lhs>::Scalar Scalar;
00141
00142 enum {
00143 LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
00144 LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
00145 ProcessFirstHalf = LhsIsSelfAdjoint
00146 && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
00147 || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
00148 || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
00149 ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
00150 };
00151
00152
00153 for (Index col = 0; col < rhs.cols(); col++) {
00154 for (Index row = 0; row < lhs.rows(); row++) {
00155 dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
00156 }
00157 }
00158
00159 for (Index row = 0; row < lhs.rows(); row++) {
00160 typename _Lhs::InnerLowerIterator lIt(lhs, row);
00161 const Index stop = lIt.col() + lIt.size();
00162 for (Index col = 0; col < rhs.cols(); col++) {
00163
00164 Index k = lIt.col();
00165 Scalar tmp = 0;
00166 while (k < stop) {
00167 tmp +=
00168 lIt.value() *
00169 rhs(k++, col);
00170 ++lIt;
00171 }
00172 dst(row, col) += tmp;
00173 lIt += -lIt.size();
00174 }
00175
00176 }
00177
00178
00179 for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
00180 typename _Lhs::InnerUpperIterator uIt(lhs, lhscol);
00181 const Index stop = uIt.size() + uIt.row();
00182 for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
00183
00184
00185 const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
00186 Index k = uIt.row();
00187 while (k < stop) {
00188 dst(k++, rhscol) +=
00189 uIt.value() *
00190 rhsCoeff;
00191 ++uIt;
00192 }
00193 uIt += -uIt.size();
00194 }
00195 }
00196
00197 }
00198
00199 template<typename Lhs, typename Rhs, typename Dest>
00200 EIGEN_DONT_INLINE void skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
00201 typedef typename remove_all<Lhs>::type _Lhs;
00202 typedef typename remove_all<Rhs>::type _Rhs;
00203 typedef typename traits<Lhs>::Scalar Scalar;
00204
00205 enum {
00206 LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
00207 LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
00208 ProcessFirstHalf = LhsIsSelfAdjoint
00209 && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
00210 || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
00211 || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
00212 ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
00213 };
00214
00215
00216 for (Index col = 0; col < rhs.cols(); col++) {
00217 for (Index row = 0; row < lhs.rows(); row++) {
00218 dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
00219 }
00220 }
00221
00222
00223 for (Index row = 0; row < lhs.rows(); row++) {
00224 typename _Lhs::InnerUpperIterator uIt(lhs, row);
00225 const Index stop = uIt.col() + uIt.size();
00226 for (Index col = 0; col < rhs.cols(); col++) {
00227
00228 Index k = uIt.col();
00229 Scalar tmp = 0;
00230 while (k < stop) {
00231 tmp +=
00232 uIt.value() *
00233 rhs(k++, col);
00234 ++uIt;
00235 }
00236
00237
00238 dst(row, col) += tmp;
00239 uIt += -uIt.size();
00240 }
00241 }
00242
00243
00244 for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
00245 typename _Lhs::InnerLowerIterator lIt(lhs, lhscol);
00246 const Index stop = lIt.size() + lIt.row();
00247 for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
00248
00249 const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
00250 Index k = lIt.row();
00251 while (k < stop) {
00252 dst(k++, rhscol) +=
00253 lIt.value() *
00254 rhsCoeff;
00255 ++lIt;
00256 }
00257 lIt += -lIt.size();
00258 }
00259 }
00260
00261 }
00262
00263 template<typename Lhs, typename Rhs, typename ResultType,
00264 int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit>
00265 struct skyline_product_selector;
00266
00267 template<typename Lhs, typename Rhs, typename ResultType>
00268 struct skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> {
00269 typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00270
00271 static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
00272 skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
00273 }
00274 };
00275
00276 template<typename Lhs, typename Rhs, typename ResultType>
00277 struct skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> {
00278 typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00279
00280 static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
00281 skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
00282 }
00283 };
00284
00285 }
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 template<typename Derived>
00301 template<typename OtherDerived >
00302 EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type
00303 SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const {
00304
00305 return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived());
00306 }
00307
00308 }
00309
00310 #endif // EIGEN_SKYLINEPRODUCT_H