SkylineProduct.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SKYLINEPRODUCT_H
26 #define EIGEN_SKYLINEPRODUCT_H
27 
28 namespace Eigen {
29 
30 template<typename Lhs, typename Rhs, int ProductMode>
31 struct SkylineProductReturnType {
32  typedef const typename internal::nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested;
33  typedef const typename internal::nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested;
34 
35  typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type;
36 };
37 
38 template<typename LhsNested, typename RhsNested, int ProductMode>
39 struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > {
40  // clean the nested types:
41  typedef typename internal::remove_all<LhsNested>::type _LhsNested;
42  typedef typename internal::remove_all<RhsNested>::type _RhsNested;
43  typedef typename _LhsNested::Scalar Scalar;
44 
45  enum {
46  LhsCoeffReadCost = _LhsNested::CoeffReadCost,
47  RhsCoeffReadCost = _RhsNested::CoeffReadCost,
48  LhsFlags = _LhsNested::Flags,
49  RhsFlags = _RhsNested::Flags,
50 
51  RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
52  ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
53  InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
54 
55  MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
56  MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
57 
58  EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
59  ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct,
60 
61  RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)),
62 
63  Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
64  | EvalBeforeAssigningBit
65  | EvalBeforeNestingBit,
66 
67  CoeffReadCost = Dynamic
68  };
69 
70  typedef typename internal::conditional<ResultIsSkyline,
71  SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >,
72  MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::type Base;
73 };
74 
75 namespace internal {
76 template<typename LhsNested, typename RhsNested, int ProductMode>
77 class SkylineProduct : no_assignment_operator,
78 public traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base {
79 public:
80 
81  EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct)
82 
83 private:
84 
85  typedef typename traits<SkylineProduct>::_LhsNested _LhsNested;
86  typedef typename traits<SkylineProduct>::_RhsNested _RhsNested;
87 
88 public:
89 
90  template<typename Lhs, typename Rhs>
91  EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs)
92  : m_lhs(lhs), m_rhs(rhs) {
93  eigen_assert(lhs.cols() == rhs.rows());
94 
95  enum {
96  ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic
97  || _RhsNested::RowsAtCompileTime == Dynamic
98  || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime),
99  AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
100  SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested)
101  };
102  // note to the lost user:
103  // * for a dot product use: v1.dot(v2)
104  // * for a coeff-wise product use: v1.cwise()*v2
105  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
106  INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
107  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
108  INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
109  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
110  }
111 
112  EIGEN_STRONG_INLINE Index rows() const {
113  return m_lhs.rows();
114  }
115 
116  EIGEN_STRONG_INLINE Index cols() const {
117  return m_rhs.cols();
118  }
119 
120  EIGEN_STRONG_INLINE const _LhsNested& lhs() const {
121  return m_lhs;
122  }
123 
124  EIGEN_STRONG_INLINE const _RhsNested& rhs() const {
125  return m_rhs;
126  }
127 
128 protected:
129  LhsNested m_lhs;
130  RhsNested m_rhs;
131 };
132 
133 // dense = skyline * dense
134 // Note that here we force no inlining and separate the setZero() because GCC messes up otherwise
135 
136 template<typename Lhs, typename Rhs, typename Dest>
137 EIGEN_DONT_INLINE void skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
138  typedef typename remove_all<Lhs>::type _Lhs;
139  typedef typename remove_all<Rhs>::type _Rhs;
140  typedef typename traits<Lhs>::Scalar Scalar;
141 
142  enum {
143  LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
144  LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
145  ProcessFirstHalf = LhsIsSelfAdjoint
146  && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
147  || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
148  || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
149  ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
150  };
151 
152  //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
153  for (Index col = 0; col < rhs.cols(); col++) {
154  for (Index row = 0; row < lhs.rows(); row++) {
155  dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
156  }
157  }
158  //Use matrix lower triangular part
159  for (Index row = 0; row < lhs.rows(); row++) {
160  typename _Lhs::InnerLowerIterator lIt(lhs, row);
161  const Index stop = lIt.col() + lIt.size();
162  for (Index col = 0; col < rhs.cols(); col++) {
163 
164  Index k = lIt.col();
165  Scalar tmp = 0;
166  while (k < stop) {
167  tmp +=
168  lIt.value() *
169  rhs(k++, col);
170  ++lIt;
171  }
172  dst(row, col) += tmp;
173  lIt += -lIt.size();
174  }
175 
176  }
177 
178  //Use matrix upper triangular part
179  for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
180  typename _Lhs::InnerUpperIterator uIt(lhs, lhscol);
181  const Index stop = uIt.size() + uIt.row();
182  for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
183 
184 
185  const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
186  Index k = uIt.row();
187  while (k < stop) {
188  dst(k++, rhscol) +=
189  uIt.value() *
190  rhsCoeff;
191  ++uIt;
192  }
193  uIt += -uIt.size();
194  }
195  }
196 
197 }
198 
199 template<typename Lhs, typename Rhs, typename Dest>
200 EIGEN_DONT_INLINE void skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) {
201  typedef typename remove_all<Lhs>::type _Lhs;
202  typedef typename remove_all<Rhs>::type _Rhs;
203  typedef typename traits<Lhs>::Scalar Scalar;
204 
205  enum {
206  LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit,
207  LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit,
208  ProcessFirstHalf = LhsIsSelfAdjoint
209  && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0)
210  || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor)
211  || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)),
212  ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
213  };
214 
215  //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix.
216  for (Index col = 0; col < rhs.cols(); col++) {
217  for (Index row = 0; row < lhs.rows(); row++) {
218  dst(row, col) = lhs.coeffDiag(row) * rhs(row, col);
219  }
220  }
221 
222  //Use matrix upper triangular part
223  for (Index row = 0; row < lhs.rows(); row++) {
224  typename _Lhs::InnerUpperIterator uIt(lhs, row);
225  const Index stop = uIt.col() + uIt.size();
226  for (Index col = 0; col < rhs.cols(); col++) {
227 
228  Index k = uIt.col();
229  Scalar tmp = 0;
230  while (k < stop) {
231  tmp +=
232  uIt.value() *
233  rhs(k++, col);
234  ++uIt;
235  }
236 
237 
238  dst(row, col) += tmp;
239  uIt += -uIt.size();
240  }
241  }
242 
243  //Use matrix lower triangular part
244  for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) {
245  typename _Lhs::InnerLowerIterator lIt(lhs, lhscol);
246  const Index stop = lIt.size() + lIt.row();
247  for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) {
248 
249  const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol);
250  Index k = lIt.row();
251  while (k < stop) {
252  dst(k++, rhscol) +=
253  lIt.value() *
254  rhsCoeff;
255  ++lIt;
256  }
257  lIt += -lIt.size();
258  }
259  }
260 
261 }
262 
263 template<typename Lhs, typename Rhs, typename ResultType,
264  int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit>
265  struct skyline_product_selector;
266 
267 template<typename Lhs, typename Rhs, typename ResultType>
268 struct skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> {
269  typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
270 
271  static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
272  skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
273  }
274 };
275 
276 template<typename Lhs, typename Rhs, typename ResultType>
277 struct skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> {
278  typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
279 
280  static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) {
281  skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res);
282  }
283 };
284 
285 } // end namespace internal
286 
287 // template<typename Derived>
288 // template<typename Lhs, typename Rhs >
289 // Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) {
290 // typedef typename internal::remove_all<Lhs>::type _Lhs;
291 // internal::skyline_product_selector<typename internal::remove_all<Lhs>::type,
292 // typename internal::remove_all<Rhs>::type,
293 // Derived>::run(product.lhs(), product.rhs(), derived());
294 //
295 // return derived();
296 // }
297 
298 // skyline * dense
299 
300 template<typename Derived>
301 template<typename OtherDerived >
302 EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type
303 SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const {
304 
305  return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived());
306 }
307 
308 } // end namespace Eigen
309 
310 #endif // EIGEN_SKYLINEPRODUCT_H