SkylineProduct.h
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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     // clean the nested types:
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         // note to the lost user:
00103         //    * for a dot product use: v1.dot(v2)
00104         //    * for a coeff-wise product use: v1.cwise()*v2
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 // dense = skyline * dense
00134 // Note that here we force no inlining and separate the setZero() because GCC messes up otherwise
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     //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
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     //Use matrix lower triangular part
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     //Use matrix upper triangular part
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     //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
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     //Use matrix upper triangular part
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     //Use matrix lower triangular part
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 } // end namespace internal
00286 
00287 // template<typename Derived>
00288 // template<typename Lhs, typename Rhs >
00289 // Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) {
00290 //     typedef typename internal::remove_all<Lhs>::type _Lhs;
00291 //     internal::skyline_product_selector<typename internal::remove_all<Lhs>::type,
00292 //             typename internal::remove_all<Rhs>::type,
00293 //             Derived>::run(product.lhs(), product.rhs(), derived());
00294 // 
00295 //     return derived();
00296 // }
00297 
00298 // skyline * dense
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 } // end namespace Eigen
00309 
00310 #endif // EIGEN_SKYLINEPRODUCT_H