1 files changed, 295 insertions, 0 deletions
diff --git a/unsupported/Eigen/src/Skyline/SkylineProduct.h b/unsupported/Eigen/src/Skyline/SkylineProduct.h
new file mode 100644
index 000000000..1ddf455e2
--- /dev/null
+++ b/unsupported/Eigen/src/Skyline/SkylineProduct.h
@@ -0,0 +1,295 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_SKYLINEPRODUCT_H
+#define EIGEN_SKYLINEPRODUCT_H
+
+namespace Eigen { 
+
+template<typename Lhs, typename Rhs, int ProductMode>
+struct SkylineProductReturnType {
+    typedef const typename internal::nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested;
+    typedef const typename internal::nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested;
+
+    typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type;
+};
+
+template<typename LhsNested, typename RhsNested, int ProductMode>
+struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > {
+    // clean the nested types:
+    typedef typename internal::remove_all<LhsNested>::type _LhsNested;
+    typedef typename internal::remove_all<RhsNested>::type _RhsNested;
+    typedef typename _LhsNested::Scalar Scalar;
+
+    enum {
+        LhsCoeffReadCost = _LhsNested::CoeffReadCost,
+        RhsCoeffReadCost = _RhsNested::CoeffReadCost,
+        LhsFlags = _LhsNested::Flags,
+        RhsFlags = _RhsNested::Flags,
+
+        RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
+        ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
+        InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
+
+        MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
+        MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
+
+        EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
+        ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct,
+
+        RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)),
+
+        Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
+        | EvalBeforeAssigningBit
+        | EvalBeforeNestingBit,
+
+        CoeffReadCost = Dynamic
+    };
+
+    typedef typename internal::conditional<ResultIsSkyline,
+            SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >,
+            MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::type Base;
+};
+
+namespace internal {
+template<typename LhsNested, typename RhsNested, int ProductMode>
+class SkylineProduct : no_assignment_operator,
+public traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base {
+public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct)
+
+private:
+
+    typedef typename traits<SkylineProduct>::_LhsNested _LhsNested;
+    typedef typename traits<SkylineProduct>::_RhsNested _RhsNested;
+
+public:
+
+    template<typename Lhs, typename Rhs>
+    EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs)
+    : m_lhs(lhs), m_rhs(rhs) {
+        eigen_assert(lhs.cols() == rhs.rows());
+
+        enum {
+            ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic
+            || _RhsNested::RowsAtCompileTime == Dynamic
+            || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime),
+            AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
+            SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested)
+        };
+        // note to the lost user:
+        //    * for a dot product use: v1.dot(v2)
+        //    * for a coeff-wise product use: v1.cwise()*v2
+        EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
+                INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
+                EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
+                INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
+                EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
+    }
+
+    EIGEN_STRONG_INLINE Index rows() const {
+        return m_lhs.rows();
+    }
+
+    EIGEN_STRONG_INLINE Index cols() const {
+        return m_rhs.cols();
+    }
+
+    EIGEN_STRONG_INLINE const _LhsNested& lhs() const {
+        return m_lhs;
+    }
+
+    EIGEN_STRONG_INLINE const _RhsNested& rhs() const {
+        return m_rhs;
+    }
+
+protected:
+    LhsNested m_lhs;
+    RhsNested m_rhs;
+};
+
+// dense = skyline * dense
+// Note that here we force no inlining and separate the setZero() because GCC messes up otherwise
+
+template<typename Lhs, typename Rhs, typename Dest>
+EIGEN_DONT_INLINE void skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
+    typedef typename remove_all<Lhs>::type _Lhs;
+    typedef typename remove_all<Rhs>::type _Rhs;
+    typedef typename traits<Lhs>::Scalar Scalar;
+
+    enum {
+        LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
+        LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
+        ProcessFirstHalf = LhsIsSelfAdjoint
+        && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
+        || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
+        || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
+        ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
+    };
+
+    //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
+    for (Index col = 0; col < rhs.cols(); col++) {
+        for (Index row = 0; row < lhs.rows(); row++) {
+            dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
+        }
+    }
+    //Use matrix lower triangular part
+    for (Index row = 0; row < lhs.rows(); row++) {
+        typename _Lhs::InnerLowerIterator lIt(lhs, row);
+        const Index stop = lIt.col() + lIt.size();
+        for (Index col = 0; col < rhs.cols(); col++) {
+
+            Index k = lIt.col();
+            Scalar tmp = 0;
+            while (k < stop) {
+                tmp +=
+                        lIt.value() *
+                        rhs(k++, col);
+                ++lIt;
+            }
+            dst(row, col) += tmp;
+            lIt += -lIt.size();
+        }
+
+    }
+
+    //Use matrix upper triangular part
+    for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
+        typename _Lhs::InnerUpperIterator uIt(lhs, lhscol);
+        const Index stop = uIt.size() + uIt.row();
+        for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
+
+
+            const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
+            Index k = uIt.row();
+            while (k < stop) {
+                dst(k++, rhscol) +=
+                        uIt.value() *
+                        rhsCoeff;
+                ++uIt;
+            }
+            uIt += -uIt.size();
+        }
+    }
+
+}
+
+template<typename Lhs, typename Rhs, typename Dest>
+EIGEN_DONT_INLINE void skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
+    typedef typename remove_all<Lhs>::type _Lhs;
+    typedef typename remove_all<Rhs>::type _Rhs;
+    typedef typename traits<Lhs>::Scalar Scalar;
+
+    enum {
+        LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
+        LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
+        ProcessFirstHalf = LhsIsSelfAdjoint
+        && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
+        || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
+        || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
+        ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
+    };
+
+    //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
+    for (Index col = 0; col < rhs.cols(); col++) {
+        for (Index row = 0; row < lhs.rows(); row++) {
+            dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
+        }
+    }
+
+    //Use matrix upper triangular part
+    for (Index row = 0; row < lhs.rows(); row++) {
+        typename _Lhs::InnerUpperIterator uIt(lhs, row);
+        const Index stop = uIt.col() + uIt.size();
+        for (Index col = 0; col < rhs.cols(); col++) {
+
+            Index k = uIt.col();
+            Scalar tmp = 0;
+            while (k < stop) {
+                tmp +=
+                        uIt.value() *
+                        rhs(k++, col);
+                ++uIt;
+            }
+
+
+            dst(row, col) += tmp;
+            uIt += -uIt.size();
+        }
+    }
+
+    //Use matrix lower triangular part
+    for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
+        typename _Lhs::InnerLowerIterator lIt(lhs, lhscol);
+        const Index stop = lIt.size() + lIt.row();
+        for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
+
+            const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
+            Index k = lIt.row();
+            while (k < stop) {
+                dst(k++, rhscol) +=
+                        lIt.value() *
+                        rhsCoeff;
+                ++lIt;
+            }
+            lIt += -lIt.size();
+        }
+    }
+
+}
+
+template<typename Lhs, typename Rhs, typename ResultType,
+        int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit>
+        struct skyline_product_selector;
+
+template<typename Lhs, typename Rhs, typename ResultType>
+struct skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> {
+    typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
+
+    static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
+        skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
+    }
+};
+
+template<typename Lhs, typename Rhs, typename ResultType>
+struct skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> {
+    typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
+
+    static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
+        skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
+    }
+};
+
+} // end namespace internal
+
+// template<typename Derived>
+// template<typename Lhs, typename Rhs >
+// Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) {
+//     typedef typename internal::remove_all<Lhs>::type _Lhs;
+//     internal::skyline_product_selector<typename internal::remove_all<Lhs>::type,
+//             typename internal::remove_all<Rhs>::type,
+//             Derived>::run(product.lhs(), product.rhs(), derived());
+// 
+//     return derived();
+// }
+
+// skyline * dense
+
+template<typename Derived>
+template<typename OtherDerived >
+EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type
+SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const {
+
+    return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived());
+}
+
+} // end namespace Eigen
+
+#endif // EIGEN_SKYLINEPRODUCT_H