SkylineProduct.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_SKYLINEPRODUCT_H
00011 #define EIGEN_SKYLINEPRODUCT_H
00012 
00013 namespace Eigen { 
00014 
00015 template<typename Lhs, typename Rhs, int ProductMode>
00016 struct SkylineProductReturnType {
00017     typedef const typename internal::nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested;
00018     typedef const typename internal::nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested;
00019 
00020     typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type;
00021 };
00022 
00023 template<typename LhsNested, typename RhsNested, int ProductMode>
00024 struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > {
00025     // clean the nested types:
00026     typedef typename internal::remove_all<LhsNested>::type _LhsNested;
00027     typedef typename internal::remove_all<RhsNested>::type _RhsNested;
00028     typedef typename _LhsNested::Scalar Scalar;
00029 
00030     enum {
00031         LhsCoeffReadCost = _LhsNested::CoeffReadCost,
00032         RhsCoeffReadCost = _RhsNested::CoeffReadCost,
00033         LhsFlags = _LhsNested::Flags,
00034         RhsFlags = _RhsNested::Flags,
00035 
00036         RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
00037         ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
00038         InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
00039 
00040         MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
00041         MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
00042 
00043         EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
00044         ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct,
00045 
00046         RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)),
00047 
00048         Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
00049         | EvalBeforeAssigningBit
00050         | EvalBeforeNestingBit,
00051 
00052         CoeffReadCost = Dynamic
00053     };
00054 
00055     typedef typename internal::conditional<ResultIsSkyline,
00056             SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >,
00057             MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::type Base;
00058 };
00059 
00060 namespace internal {
00061 template<typename LhsNested, typename RhsNested, int ProductMode>
00062 class SkylineProduct : no_assignment_operator,
00063 public traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base {
00064 public:
00065 
00066     EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct)
00067 
00068 private:
00069 
00070     typedef typename traits<SkylineProduct>::_LhsNested _LhsNested;
00071     typedef typename traits<SkylineProduct>::_RhsNested _RhsNested;
00072 
00073 public:
00074 
00075     template<typename Lhs, typename Rhs>
00076     EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs)
00077     : m_lhs(lhs), m_rhs(rhs) {
00078         eigen_assert(lhs.cols() == rhs.rows());
00079 
00080         enum {
00081             ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic
00082             || _RhsNested::RowsAtCompileTime == Dynamic
00083             || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime),
00084             AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
00085             SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested)
00086         };
00087         // note to the lost user:
00088         //    * for a dot product use: v1.dot(v2)
00089         //    * for a coeff-wise product use: v1.cwise()*v2
00090         EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00091                 INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00092                 EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00093                 INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00094                 EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00095     }
00096 
00097     EIGEN_STRONG_INLINE Index rows() const {
00098         return m_lhs.rows();
00099     }
00100 
00101     EIGEN_STRONG_INLINE Index cols() const {
00102         return m_rhs.cols();
00103     }
00104 
00105     EIGEN_STRONG_INLINE const _LhsNested& lhs() const {
00106         return m_lhs;
00107     }
00108 
00109     EIGEN_STRONG_INLINE const _RhsNested& rhs() const {
00110         return m_rhs;
00111     }
00112 
00113 protected:
00114     LhsNested m_lhs;
00115     RhsNested m_rhs;
00116 };
00117 
00118 // dense = skyline * dense
00119 // Note that here we force no inlining and separate the setZero() because GCC messes up otherwise
00120 
00121 template<typename Lhs, typename Rhs, typename Dest>
00122 EIGEN_DONT_INLINE void skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
00123     typedef typename remove_all<Lhs>::type _Lhs;
00124     typedef typename remove_all<Rhs>::type _Rhs;
00125     typedef typename traits<Lhs>::Scalar Scalar;
00126 
00127     enum {
00128         LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
00129         LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
00130         ProcessFirstHalf = LhsIsSelfAdjoint
00131         && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
00132         || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
00133         || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
00134         ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
00135     };
00136 
00137     //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
00138     for (Index col = 0; col < rhs.cols(); col++) {
00139         for (Index row = 0; row < lhs.rows(); row++) {
00140             dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
00141         }
00142     }
00143     //Use matrix lower triangular part
00144     for (Index row = 0; row < lhs.rows(); row++) {
00145         typename _Lhs::InnerLowerIterator lIt(lhs, row);
00146         const Index stop = lIt.col() + lIt.size();
00147         for (Index col = 0; col < rhs.cols(); col++) {
00148 
00149             Index k = lIt.col();
00150             Scalar tmp = 0;
00151             while (k < stop) {
00152                 tmp +=
00153                         lIt.value() *
00154                         rhs(k++, col);
00155                 ++lIt;
00156             }
00157             dst(row, col) += tmp;
00158             lIt += -lIt.size();
00159         }
00160 
00161     }
00162 
00163     //Use matrix upper triangular part
00164     for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
00165         typename _Lhs::InnerUpperIterator uIt(lhs, lhscol);
00166         const Index stop = uIt.size() + uIt.row();
00167         for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
00168 
00169 
00170             const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
00171             Index k = uIt.row();
00172             while (k < stop) {
00173                 dst(k++, rhscol) +=
00174                         uIt.value() *
00175                         rhsCoeff;
00176                 ++uIt;
00177             }
00178             uIt += -uIt.size();
00179         }
00180     }
00181 
00182 }
00183 
00184 template<typename Lhs, typename Rhs, typename Dest>
00185 EIGEN_DONT_INLINE void skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
00186     typedef typename remove_all<Lhs>::type _Lhs;
00187     typedef typename remove_all<Rhs>::type _Rhs;
00188     typedef typename traits<Lhs>::Scalar Scalar;
00189 
00190     enum {
00191         LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
00192         LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
00193         ProcessFirstHalf = LhsIsSelfAdjoint
00194         && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
00195         || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
00196         || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
00197         ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
00198     };
00199 
00200     //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
00201     for (Index col = 0; col < rhs.cols(); col++) {
00202         for (Index row = 0; row < lhs.rows(); row++) {
00203             dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
00204         }
00205     }
00206 
00207     //Use matrix upper triangular part
00208     for (Index row = 0; row < lhs.rows(); row++) {
00209         typename _Lhs::InnerUpperIterator uIt(lhs, row);
00210         const Index stop = uIt.col() + uIt.size();
00211         for (Index col = 0; col < rhs.cols(); col++) {
00212 
00213             Index k = uIt.col();
00214             Scalar tmp = 0;
00215             while (k < stop) {
00216                 tmp +=
00217                         uIt.value() *
00218                         rhs(k++, col);
00219                 ++uIt;
00220             }
00221 
00222 
00223             dst(row, col) += tmp;
00224             uIt += -uIt.size();
00225         }
00226     }
00227 
00228     //Use matrix lower triangular part
00229     for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
00230         typename _Lhs::InnerLowerIterator lIt(lhs, lhscol);
00231         const Index stop = lIt.size() + lIt.row();
00232         for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
00233 
00234             const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
00235             Index k = lIt.row();
00236             while (k < stop) {
00237                 dst(k++, rhscol) +=
00238                         lIt.value() *
00239                         rhsCoeff;
00240                 ++lIt;
00241             }
00242             lIt += -lIt.size();
00243         }
00244     }
00245 
00246 }
00247 
00248 template<typename Lhs, typename Rhs, typename ResultType,
00249         int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit>
00250         struct skyline_product_selector;
00251 
00252 template<typename Lhs, typename Rhs, typename ResultType>
00253 struct skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> {
00254     typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00255 
00256     static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
00257         skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
00258     }
00259 };
00260 
00261 template<typename Lhs, typename Rhs, typename ResultType>
00262 struct skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> {
00263     typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00264 
00265     static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
00266         skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
00267     }
00268 };
00269 
00270 } // end namespace internal
00271 
00272 // template<typename Derived>
00273 // template<typename Lhs, typename Rhs >
00274 // Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) {
00275 //     typedef typename internal::remove_all<Lhs>::type _Lhs;
00276 //     internal::skyline_product_selector<typename internal::remove_all<Lhs>::type,
00277 //             typename internal::remove_all<Rhs>::type,
00278 //             Derived>::run(product.lhs(), product.rhs(), derived());
00279 // 
00280 //     return derived();
00281 // }
00282 
00283 // skyline * dense
00284 
00285 template<typename Derived>
00286 template<typename OtherDerived >
00287 EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type
00288 SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const {
00289 
00290     return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived());
00291 }
00292 
00293 } // end namespace Eigen
00294 
00295 #endif // EIGEN_SKYLINEPRODUCT_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:00:17