Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
alt_matrix.hpp
1 /*
2  This file is part of the Kernel Quantum Probability library (KQP).
3 
4  KQP is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  KQP is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with KQP. If not, see <http://www.gnu.org/licenses/>.
16  */
17 #ifndef _KQP_ALTMATRIX_H_
18 #define _KQP_ALTMATRIX_H_
19 
20 #include <boost/shared_ptr.hpp>
21 #include <algorithm>
22 
23 #include <kqp/kqp.hpp>
24 #include <kqp/eigen_identity.hpp>
25 #include <Eigen/Core>
26 #include <Eigen/Sparse>
27 
28 namespace kqp {
29 
30  // Should use template template aliases?
31 
32 
33  // Forward declarations
34  template<typename T1, typename T2> class AltMatrix;
35  template<typename _AltMatrix> class AltBlock;
36  template<typename _AltMatrix> class RowWise;
37  template<typename Derived> struct scalar;
38  template<typename Derived> struct DiagonalBlockWrapper;
39  template<typename Derived> class AltArrayWrapper;
40 
41  // Returns the associated scalar
42  template<typename Derived> struct scalar { typedef typename Eigen::internal::traits<Derived>::Scalar type; };
43  template<> struct scalar<double> { typedef double type; };
44  template<> struct scalar<float> { typedef float type; };
45  template<> struct scalar<std::complex<double>> { typedef std::complex<double> type; };
46  template<> struct scalar<std::complex<float>> { typedef std::complex<float> type; };
47 
48 
49  // ---- Predefined Alt-based matrices
50 
52  template<typename Scalar> struct AltDiagonal {
53  typedef Eigen::DiagonalWrapper<Eigen::Matrix<Scalar,Dynamic,1>> Diagonal;
54  typedef typename Eigen::MatrixBase<Eigen::Matrix<Scalar,Dynamic,Dynamic> >::IdentityReturnType Identity;
55 
57  };
58 
60  template<typename Scalar> struct AltVector {
61  typedef Eigen::Matrix<Scalar,Dynamic,1> VectorType;
62  typedef typename Eigen::Matrix<Scalar,Dynamic,1>::ConstantReturnType ConstantVectorType;
63 
65  };
66 
67 
69  template<typename Scalar> struct AltDense {
70  typedef Eigen::Matrix<Scalar,Dynamic,Dynamic> DenseType;
73 
74  static inline type Identity(Index n) {
75  return type(Eigen::Matrix<Scalar,Dynamic,Dynamic>::Identity(n,n));
76  }
77  };
78 
79 
81  struct AltMatrixStorage {};
82 
83  template<typename Operator, typename Derived> class AltCwiseUnaryOp;
84  template<typename Derived> class AltAsDiagonal;
85 
87  template<typename Derived>
88  class AltMatrixBase : public Eigen::EigenBase<Derived> {
89  public:
90  typedef typename Eigen::internal::traits<Derived>::Scalar Scalar;
91 
92  AltCwiseUnaryOp<Eigen::internal::scalar_sqrt_op<Scalar>, const Derived> cwiseSqrt() const {
93  return this->derived();
94  }
95 
96  AltCwiseUnaryOp<Eigen::internal::scalar_abs_op<Scalar>, const Derived> cwiseAbs() const {
97  return this->derived();
98  }
99 
100  AltCwiseUnaryOp<Eigen::internal::scalar_abs2_op<Scalar>, const Derived> cwiseAbs2() const {
101  return this->derived();
102  }
103 
104  AltCwiseUnaryOp<Eigen::internal::scalar_inverse_op<Scalar>, const Derived> cwiseInverse() const {
105  return this->derived();
106  }
107 
108  AltAsDiagonal<const Derived> asDiagonal() const {
109  return this->derived();
110  }
111 
112  AltArrayWrapper<const Derived> array() const {
113  return this->derived();
114  }
115 
116  };
117 
118 
120  template<typename MatrixType, unsigned int UpLo, typename Derived>
121  void rankUpdate(Eigen::SelfAdjointView<MatrixType, UpLo> &&matrix, const Eigen::MatrixBase<Derived> &mA, const typename MatrixType::Scalar alpha) {
122  matrix.rankUpdate(mA, alpha);
123  }
124 
131  template<typename MatrixType, unsigned int UpLo, typename Derived>
132  void rankUpdate2(Eigen::SelfAdjointView<MatrixType, UpLo> &&matrix, const AltMatrixBase<Derived> &mA, const typename MatrixType::Scalar alpha) {
133  if (mA.derived().isT1())
134  rankUpdate(std::move(matrix), mA.derived().t1(), alpha);
135  else
136  rankUpdate(std::move(matrix), mA.derived().t2(), alpha);
137  }
138 
139 
140  // --- As diagonal
141  template<typename XprType> class AltAsDiagonal : Eigen::internal::no_assignment_operator, public AltMatrixBase<AltAsDiagonal<XprType>> {
142  protected:
143  const typename XprType::Nested m_xpr;
144  public:
145 
146  inline AltAsDiagonal(const XprType& xpr) : m_xpr(xpr) {
147  }
148 
149  EIGEN_STRONG_INLINE Index rows() const { return m_xpr.rows(); }
150  EIGEN_STRONG_INLINE Index cols() const { return m_xpr.rows(); }
151 
152  bool isT1() const { return m_xpr.derived().isT1(); }
153 
154  auto t1() const -> decltype(m_xpr.derived().t1().asDiagonal()) { return m_xpr.derived().t1().asDiagonal(); }
155  auto t2() const -> decltype(m_xpr.derived().t2().asDiagonal()) { return m_xpr.derived().t2().asDiagonal(); }
156 
157  template<typename Dest> inline void evalTo(Dest& dst) const
158  { if (isT1()) dst = t1(); else dst = t2(); }
159 
160 
161  };
162 
163  // --- Unary operator
164  template<typename UnaryOp, typename XprType> class AltCwiseUnaryOp :
165  Eigen::internal::no_assignment_operator, public AltMatrixBase<AltCwiseUnaryOp<UnaryOp,XprType>> {
166  protected:
167  const typename XprType::Nested m_xpr;
168  const UnaryOp m_functor;
169  public:
171  typedef typename Eigen::internal::traits<Nested>::Scalar Scalar;
172 
173  inline AltCwiseUnaryOp(const XprType& xpr, const UnaryOp& func = UnaryOp())
174  : m_xpr(xpr), m_functor(func) {
175  }
176 
177  EIGEN_STRONG_INLINE Index rows() const { return m_xpr.rows(); }
178  EIGEN_STRONG_INLINE Index cols() const { return m_xpr.cols(); }
179 
180  bool isT1() const { return m_xpr.derived().isT1(); }
181 
182  typedef typename Eigen::internal::remove_all<decltype(m_xpr.derived().t1())>::type T1;
183  Eigen::CwiseUnaryOp<UnaryOp,T1> t1() const {
184  return Eigen::CwiseUnaryOp<UnaryOp,T1>(m_xpr.derived().t1(), m_functor);
185  }
186 
187  typedef typename Eigen::internal::remove_all<decltype(m_xpr.derived().t2())>::type T2;
188  Eigen::CwiseUnaryOp<UnaryOp,T2> t2() const {
189  return Eigen::CwiseUnaryOp<UnaryOp,T2>(m_xpr.derived().t2(), m_functor);
190  }
191 
192  template<typename Dest> inline void evalTo(Dest& dst) const {
193  if (isT1()) dst = t1(); else dst = t2();
194  }
195 
196  Scalar sum() const {
197  if (isT1()) return t1().sum();
198  else return t2().sum();
199  }
200  };
201 
202 
203 
204  // --- Helper functions
205  template<typename Derived>
206  const typename Eigen::MatrixBase<Derived>::AdjointReturnType adjoint(const Eigen::MatrixBase<Derived>& x) { return x.adjoint(); }
207 
208  template<typename Derived>
209  Eigen::DiagonalWrapper<const Derived> adjoint(const Eigen::DiagonalWrapper<const Derived> & x) { return x; }
210 
211  template<typename Scalar>
212  Eigen::Identity<Scalar> adjoint(const Eigen::Identity<Scalar> & x) { return x; }
213 
214  template<typename Derived>
215  auto blockSquaredNorm(const Eigen::MatrixBase<Derived>& x, Index row, Index col, Index rowSize, Index colSize) -> decltype(x.squaredNorm()) {
216  return x.block(row,col,rowSize,colSize).squaredNorm();
217  }
218 
219  template<typename Derived>
220  auto blockSquaredNorm(const Eigen::DiagonalWrapper<Derived> & x, Index row, Index col, Index rowSize, Index colSize) -> decltype(x.diagonal().squaredNorm()) {
221  return x.diagonal().segment(std::max(row,col), std::min(row+rowSize, col+colSize)).squaredNorm();
222  }
223 
224  template<typename Scalar>
225  Scalar blockSquaredNorm(const Eigen::Identity<Scalar> &, Index row, Index col, Index rowSize, Index colSize) {
226  return std::max(row,col) - std::min(row+rowSize, col+colSize);
227  }
228 
229  template<typename Derived>
230  auto squaredNorm(const Eigen::MatrixBase<Derived>& x) -> decltype(x.derived().squaredNorm()) { return x.derived().squaredNorm(); }
231 
232  template<typename Derived>
233  auto squaredNorm(const Eigen::DiagonalWrapper<Derived> & x) -> decltype(x.diagonal().squaredNorm()) { return x.diagonal().squaredNorm(); }
234 
235  // ---- Transpose
236 
237 
238 
240  template<typename Derived>
241  struct Adjoint : public AltMatrixBase< Adjoint<Derived> > {
242  typename Eigen::internal::ref_selector<Derived>::type nested;
243  public:
244 
245  Adjoint(Derived &nested) : nested(nested) {}
246 
247  Index rows() const { return nested.cols(); }
248  Index cols() const { return nested.rows(); }
249 
250  bool isT1() const { return nested.isT1(); }
251 
252  typedef decltype(kqp::adjoint(nested.t1())) T1;
253  typedef decltype(kqp::adjoint(nested.t2())) T2;
254 
255  inline T1 t1() const { return adjoint(nested.t1()); }
256  inline T2 t2() const { return adjoint(nested.t2()); }
257 
258 
259  void printExpression(std::ostream &out) const {
260  out << "Adjoint(";
261  nested.printExpression(out);
262  out << ")";
263  }
264 
265 
266  };
267 
268 
269 
270  template <typename Derived>
271  void printExpression(std::ostream &out, const kqp::AltMatrixBase<Derived> &o) {
272  (static_cast<const Derived&>(o)).printExpression(out);
273  }
274 
275  template <typename Derived>
276  void printExpression(std::ostream &out, const Eigen::EigenBase<Derived> &x) {
277  out << KQP_DEMANGLE(x) << "[" << x.rows() << " x " << x.cols() << "]";
278  }
279 
280  template <typename Scalar>
281  void printExpression(std::ostream &out, const Eigen::Identity<Scalar> &x) {
282  out << "Id[" << KQP_DEMANGLE(Scalar) << "; " << x.rows() << " x " << x.cols() << "]";
283  }
284 
285  template <typename Derived>
286  void printExpression(std::ostream &out, const Eigen::Transpose<Derived> &x) {
287  out << "adjoint(";
288  printExpression(out, static_cast<const Derived&>(x.derived()));
289  out << ")";
290  }
291 
292  template <typename Derived>
293  void printExpression(std::ostream &out, const Eigen::DiagonalWrapper<Derived> &x) {
294  out << "diag[" << x.rows() << " x " << x.cols() << "]";
295  }
296 
297  template <typename Scalar, int a, int b>
298  void printExpression(std::ostream &out, const Eigen::Matrix<Scalar, a, b> &x) {
299  out << "dense/" << &x << "[" << x.rows() << " x " << x.cols() << "]";
300  }
301 
302  template <typename Derived, typename OtherDerived>
303  void printExpression(std::ostream &out, const Eigen::DiagonalProduct<Derived, OtherDerived, Eigen::OnTheLeft> &x) {
304  out << "dproduct[";
305  kqp::printExpression(out, x.m_diagonal);
306  out << " x ";
307  kqp::printExpression(out, x.m_matrix);
308  out << "]";
309  }
310  template <typename Derived, typename OtherDerived>
311  void printExpression(std::ostream &out, const Eigen::DiagonalProduct<Derived, OtherDerived, Eigen::OnTheRight> &x) {
312  out << "dproduct[";
313  kqp::printExpression(out, x.m_matrix);
314  out << " x ";
315  kqp::printExpression(out, x.m_diagonal);
316  out << "]";
317  }
318 
319 
320  template <typename Derived, typename OtherDerived>
321  void printExpression(std::ostream &out, const Eigen::GeneralProduct<Derived, OtherDerived> &x) {
322  out << "product[";
323  kqp::printExpression(out, x.m_lhs);
324  out << " x ";
325  kqp::printExpression(out, x.m_rhs);
326  out << "]";
327  }
328 
329 
330 
331  // ---- Resizing
332 
334  template<class Derived>
335  typename boost::enable_if_c<(Eigen::internal::traits<Derived>::RowsAtCompileTime == Dynamic) && (Eigen::internal::traits<Derived>::ColsAtCompileTime == Dynamic), void>::type
336  resize(Eigen::EigenBase<Derived> &matrix, bool conservative, Index rows, Index cols) {
337  if (conservative) matrix.derived().conservativeResize(rows, cols); else matrix.derived().resize(rows, cols);
338  }
339 
340 
342  template<class Derived>
343  typename boost::enable_if_c<(Eigen::internal::traits<Derived>::RowsAtCompileTime == Dynamic) && (Eigen::internal::traits<Derived>::ColsAtCompileTime != Dynamic), void>::type
344  resize(Eigen::EigenBase<Derived> &matrix, bool conservative, Index rows, Index cols) {
345  if (cols != matrix.cols())
346  KQP_THROW_EXCEPTION_F(out_of_bound_exception, "Cannot change the number of columns in a fixed column-sized matrix (%d to %d)", %matrix.cols() %cols);
347  if (conservative) matrix.derived().conservativeResize(rows); else matrix.derived().resize(rows);
348  }
349 
351  template<class Derived>
352  typename boost::enable_if_c<(Eigen::internal::traits<Derived>::RowsAtCompileTime != Dynamic) && (Eigen::internal::traits<Derived>::ColsAtCompileTime == Dynamic), void>::type
353  resize(Eigen::EigenBase<Derived> &matrix, bool conservative, Index rows, Index cols) {
354  if (rows != matrix.rows())
355  KQP_THROW_EXCEPTION_F(out_of_bound_exception, "Cannot change the number of rows in a fixed row-sized matrix (%d to %d)", %matrix.rows() %rows);
356  if (conservative) matrix.derived().conservativeResize(cols); else matrix.derived().resize(cols);
357  }
358 
360  template<typename Scalar>
361  void resize(Eigen::Identity<Scalar> &matrix, bool conservative, Index rows, Index cols) {
362  if (conservative)
363  matrix.conservativeResize(rows, cols);
364  else
365  matrix.resize(rows, cols);
366  }
367 
368 
369 
370  // --- AltMatrix inner storage of Eigen matrices
371 
373  template<typename Derived>
374  struct storage {
375  typedef Derived & ReturnType;
376  typedef const Derived & ConstReturnType;
377  typedef typename Eigen::internal::traits<Derived>::Scalar Scalar;
378 
379  Derived m_value;
380 
381  storage() {}
382  storage(const Derived &value) : m_value(value) {}
383  ConstReturnType get() const { return m_value; }
384  ReturnType get() { return m_value; }
385 
386  void swap(storage &other) { m_value.swap(other.m_value); }
387  void swap(Derived &value) { m_value.swap(value); }
388  Index rows() const { return m_value.rows(); }
389  Index cols() const { return m_value.cols(); }
390 
391  Scalar trace() const { return m_value.trace(); }
392  void resize(Index rows, Index cols) {
393  kqp::resize(m_value, false, rows, cols);
394  }
395 
396  void conservativeResize(Index rows, Index cols) {
397  kqp::resize(m_value, true, rows, cols);
398  }
399 
400  typename Eigen::internal::traits<Derived>::Scalar operator()(Index i, Index j) const {
401  return m_value(i,j);
402  }
403 
404  template<typename CwiseUnaryOp>
405  void unaryExprInPlace(const CwiseUnaryOp &op) {
406  m_value = m_value.unaryExpr(op);
407  }
408 
409  typedef typename Eigen::NumTraits<typename Eigen::internal::traits<Derived>::Scalar>::Real Real;
410  Real squaredNorm() const { return m_value.squaredNorm(); }
411  Scalar sum() const { return m_value.sum(); }
412  const std::type_info &getTypeId() const { return typeid(Derived); }
413 
414  auto block(Index startRow, Index startCol, Index blockRows, Index blockCols) -> decltype(m_value.block(0,0,0,0)) {
415  return m_value.block(startRow, startCol, blockRows, blockCols);
416  }
417 
418  auto block(Index startRow, Index startCol, Index blockRows, Index blockCols) const
419  -> decltype(const_cast<ConstReturnType>(m_value).block(0,0,0,0)) {
420  return m_value.block(startRow, startCol, blockRows, blockCols);
421  }
422 
423  };
424 
426  template<typename Scalar, typename Derived>
427  struct storage< typename Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Derived> > {
428  typedef typename Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Derived> ReturnType;
429  typedef const ReturnType ConstReturnType;
430  typedef typename Eigen::NumTraits<Scalar>::Real Real;
431 
432  Scalar m_value;
433  Index m_rows;
434  Index m_cols;
435 
436  storage() : m_value(0), m_rows(0), m_cols(0) {}
437  storage(const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>,Derived> &value)
438  : m_value(value.coeff(0,0)), m_rows(value.rows()), m_cols(value.cols()) {}
439 
440  ReturnType get() const { return ReturnType(m_rows, m_cols, m_value); }
441 
442 
443  void swap(storage &other) {
444  std::swap(m_value,other.m_value);
445  std::swap(m_rows,other.m_rows);
446  std::swap(m_cols,other.m_cols);
447  }
448  void swap(ReturnType &) { KQP_THROW_EXCEPTION(illegal_argument_exception, "Cannot swap a constant matrix"); }
449  Index rows() const { return m_rows; }
450  Index cols() const { return m_cols; }
451 
452  void resize(Index rows, Index cols) {
453  m_rows = rows;
454  m_cols = cols;
455  }
456 
457  void conservativeResize(Index rows, Index cols) {
458  if (rows > m_rows || cols > m_cols)
459  KQP_THROW_EXCEPTION_F(illegal_argument_exception, "Cannot resize to a non smaller size (%d x %d to %d x %d)",
460  %m_rows%m_cols%rows%cols);
461  this->resize(rows,cols);
462  }
463 
464  Scalar operator()(Index, Index) const {
465  return m_value;
466  }
467 
468  template<typename CwiseUnaryOp>
469  void unaryExprInPlace(const CwiseUnaryOp &op) {
470  m_value = op(m_value);
471  }
472 
473  Real squaredNorm() const { return std::abs(m_value)*std::abs(m_value) * (Real)m_rows * (Real)m_cols; }
474  Scalar sum() const { return (Scalar)m_rows * (Scalar)m_cols; }
475 
476  const std::type_info &getTypeId() const { return typeid(ReturnType); }
477 
478 
479  ReturnType block(Index, Index, Index blockRows, Index blockCols) const {
480  return ReturnType(blockRows, blockCols, m_value);
481  }
482 
483  };
484 
485 
486 
488  template<typename Derived>
489  struct storage< Eigen::DiagonalWrapper<Derived> > {
490  typedef Eigen::DiagonalWrapper<const Derived> ConstReturnType;
491  typedef Eigen::DiagonalWrapper<const Derived> ReturnType;
492  typedef typename Eigen::internal::traits<Derived>::Scalar Scalar;
493  typedef typename Eigen::NumTraits<Scalar>::Real Real;
494 
495  Derived m_value;
496 
497  storage() {}
498  storage(const Eigen::DiagonalWrapper<Derived> &value) : m_value(value.diagonal()) {}
499  ConstReturnType get() const { return m_value.asDiagonal(); }
500  ReturnType get() { return static_cast<ReturnType>(m_value.asDiagonal()); }
501 
502  void swap(storage &other) {
503  m_value.swap(other.m_value);
504  }
505  void swap(ReturnType &value) { m_value.swap(value); }
506  Index rows() const { return m_value.rows(); }
507  Index cols() const { return m_value.rows(); }
508 
509  void resize(Index rows, Index cols) {
510  if (rows != cols) KQP_THROW_EXCEPTION_F(illegal_argument_exception, "Cannot resize to a non diagonal size (%d x%d)", %rows%cols);
511  m_value.resize(rows);
512  }
513 
514  void conservativeResize(Index rows, Index cols) {
515  if (rows != cols) KQP_THROW_EXCEPTION_F(illegal_argument_exception, "Cannot resize to a non diagonal size (%d x%d)", %rows%cols);
516  m_value.conservativeResize(rows);
517  }
518 
519  typename Eigen::internal::traits<Derived>::Scalar operator()(Index i, Index j) const {
520  return i == j ? m_value(i) : 0;
521  }
522 
523  template<typename CwiseUnaryOp>
524  void unaryExprInPlace(const CwiseUnaryOp &op) {
525  m_value = m_value.unaryExpr(op);
526  }
527 
528  Real squaredNorm() const { return m_value.squaredNorm(); }
529  Scalar sum() const { return m_value.sum(); }
530  const std::type_info &getTypeId() const { return typeid(ReturnType); }
531 
532  DiagonalBlockWrapper<Derived> block(Index startRow, Index startCol, Index blockRows, Index blockCols) const {
533  return DiagonalBlockWrapper<Derived>(m_value, startRow, startCol, blockRows, blockCols);
534  }
535 
536  };
537 
538 
540  template<typename Scalar>
541  struct storage< Eigen::CwiseNullaryOp<Eigen::internal::scalar_identity_op<Scalar>, Eigen::Matrix<Scalar,Dynamic,Dynamic> > > {
542  Index m_rows, m_cols;
543 
544  typedef Eigen::CwiseNullaryOp<Eigen::internal::scalar_identity_op<Scalar>, Eigen::Matrix<Scalar,Dynamic,Dynamic> > Type;
545  typedef Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Eigen::Matrix<Scalar,Dynamic,1> > VectorType;
546  typedef Type ReturnType;
547  typedef const Type ConstReturnType;
548 
549  storage() {}
550  storage(const Type &value) : m_rows(value.rows()), m_cols(value.cols()) {}
551  ReturnType get() const { return Eigen::Matrix<Scalar,Dynamic,Dynamic>::Identity(m_rows, m_cols); }
552 
553  void swap(storage &other) {
554  std::swap(m_rows, other.m_rows);
555  std::swap(m_cols, other.m_cols);
556  }
557  void swap(Type &) {
558  KQP_THROW_EXCEPTION(not_implemented_exception, "Not sure what to do");
559  }
560 
561 
562  void resize(Index rows, Index cols) {
563  m_rows = rows;
564  m_cols = cols;
565  }
566 
567  void conservativeResize(Index rows, Index cols) {
568  resize(rows, cols);
569  }
570 
571  Index rows() const { return m_rows; }
572  Index cols() const { return m_cols; }
573 
574  Scalar operator()(Index i, Index j) const {
575  return i == j ? 1 : 0;
576  }
577 
578  template<typename CwiseUnaryOp>
579  void unaryExprInPlace(const CwiseUnaryOp &) {
580  KQP_THROW_EXCEPTION(illegal_operation_exception, "Cannot modify the identity");
581  }
582 
583  typedef typename Eigen::NumTraits<Scalar>::Real Real;
584  Real squaredNorm() const { return std::min(m_rows, m_cols); }
585  Scalar sum() const { return std::min(m_rows, m_cols); }
586  const std::type_info &getTypeId() const { return typeid(Type); }
587 
588 
589 
590  DiagonalBlockWrapper<VectorType> block(Index startCol, Index startRow, Index blockRows, Index blockCols) const {
591  return DiagonalBlockWrapper<VectorType>(VectorType(std::min(m_rows,m_cols),1,1), startCol, startRow, blockRows, blockCols);
592  }
593 
594  };
595 
596 
597 
598 
600  template<typename T1, typename T2>
601  class AltMatrix : public AltMatrixBase< AltMatrix<T1, T2> > {
602  storage<T1> m_t1;
603  storage<T2> m_t2;
604 
605  bool m_isT1;
606 
607  public:
608  typedef typename Eigen::internal::traits<T1>::Scalar Scalar;
609  typedef typename Eigen::NumTraits<Scalar>::Real Real;
610  typedef const AltMatrix& Nested;
611 
612  AltMatrix() : m_isT1(true) {}
613  AltMatrix(const T1 &t1) : m_t1(t1), m_t2(), m_isT1(true) {}
614  AltMatrix(const T2 &t2) : m_t1(), m_t2(t2), m_isT1(false) {}
615 
616 
617  Index rows() const { return m_isT1 ? m_t1.rows() : m_t2.rows(); }
618  Index cols() const { return m_isT1 ? m_t1.cols() : m_t2.cols(); }
619 
620  Real squaredNorm() const { return m_isT1 ? m_t1.squaredNorm() : m_t2.squaredNorm(); }
621  Scalar trace() const { return m_isT1 ? m_t1.get().trace() : m_t2.get().trace(); }
622  Scalar sum() const { return m_isT1 ? m_t1.get().sum() : m_t2.get().sum(); }
623 
624  template<typename CwiseUnaryOp>
625  void unaryExprInPlace(const CwiseUnaryOp &op) {
626  if (isT1()) m_t1.unaryExprInPlace(op);
627  else m_t2.unaryExprInPlace(op);
628  }
629 
630  const std::type_info &getTypeId() const {
631  if (isT1()) return m_t1.getTypeId();
632  return m_t2.getTypeId();
633  }
634 
635  bool isT1() const { return m_isT1; }
636  bool isT2() const { return !m_isT1; }
637 
638  inline typename storage<T1>::ConstReturnType t1() const { return m_t1.get(); }
639  inline typename storage<T2>::ConstReturnType t2() const { return m_t2.get(); }
640 
641  inline typename storage<T1>::ReturnType t1() { return m_t1.get(); }
642  inline typename storage<T2>::ReturnType t2() { return m_t2.get(); }
643 
644  const storage<T1> &getStorage1() const { return m_t1; }
645  const storage<T2> &getStorage2() const { return m_t2; }
646 
647 
648  void swap(AltMatrix &other) {
649  m_t1.swap(other.m_t1);
650  m_t2.swap(other.m_t2);
651  std::swap(m_isT1, other.m_isT1);
652  }
653 
654  void swap(T1 &t1) { m_isT1 = true; m_t1.swap(t1); }
655  void swap(T2 &t2) { m_isT1 = true; m_t2.swap(t2); }
656 
658  Adjoint<AltMatrix> adjoint() const { return Adjoint<AltMatrix>(const_cast<AltMatrix&>(*this)); }
659 
660  Scalar operator()(Index i, Index j) const {
661  return m_isT1 ? m_t1(i,j) : m_t2(i,j);
662  }
663 
664  inline Scalar coeff(Index i, Index j) const {
665  return operator()(i,j);
666  }
667 
668  template<typename Dest> inline void evalTo(Dest& dst) const {
669  if (m_isT1)
670  dst = m_t1.get();
671  else
672  dst = m_t2.get();
673  }
674 
675 
676  void printExpression(std::ostream &out) const {
677  out << "AltMatrix/" << this << "[";
678  if (m_isT1) kqp::printExpression(out, m_t1.get()); else kqp::printExpression(out, m_t2.get());
679  out << "]";
680  }
681 
682  // ---- Resizing ----
683 
684 
685  void conservativeResize(Index rows, Index cols) {
686  if (m_isT1)
687  m_t1.conservativeResize(rows, cols);
688  else
689  m_t2.conservativeResize(rows, cols);
690  }
691 
692  void resize(Index rows, Index cols) {
693  if (m_isT1)
694  m_t1.resize(rows, cols);
695  else
696  m_t2.resize(rows, cols);
697  }
698 
699  // ---- Blocks ----
700 
701  typedef AltMatrix<T1,T2> Self;
702 
703  typedef AltBlock< Self > Block;
704  typedef AltBlock< const Self > ConstBlock;
705 
706  friend class AltBlock<Self> ;
707  friend class AltBlock<const Self>;
708 
709  Block row(Index i) { return Block(*this, i, 0, 1, cols()); }
710  ConstBlock row(Index i) const { return ConstBlock(const_cast<Self&>(*this), i, 0, 1, cols()); }
711 
712  Block col(Index j) { return Block(*this, 0, j, rows(), 1); }
713  ConstBlock col(Index j) const { return ConstBlock(const_cast<Self&>(*this), 0, j, rows(), 1); }
714 
715  Block block(Index i, Index j, Index height, Index width) { return Block(*this, i, j, height, width); }
716  ConstBlock block(Index i, Index j, Index height, Index width) const { return ConstBlock(const_cast<Self&>(*this), i, j, height, width); }
717 
718  ConstBlock topRows(Index h) const { return ConstBlock(*this, 0, 0, h, cols()); }
719  Block topRows(Index h) { return Block(*this, 0, 0, h, cols()); }
720 
721  ConstBlock bottomRows(Index h) const { return ConstBlock(*this, rows() - h, 0, h, cols()); }
722  Block bottomRows(Index h) { return Block(*this, rows() - h, 0, h, cols()); }
723 
724  const RowWise<AltMatrix> rowwise() const {
725  return RowWise<AltMatrix>(const_cast<Self&>(*this));
726  }
727 
728  };
729 
730 
731  template<typename T1, typename T2>
732  std::ostream &operator<<(std::ostream &out, const AltMatrix<T1,T2> &alt_matrix) {
733  return out << "Alt/" << (alt_matrix.isT1() ? kqp::demangle(typeid(T1)) : kqp::demangle(typeid(T2)));
734  }
735 
736 
737  template <typename Derived> auto rowwise(const Derived &x) -> decltype(x.rowwise()) {
738  return x.rowwise();
739  }
740 
741  template <typename Derived> auto rowwise(const Eigen::DiagonalWrapper<Derived> &x) -> decltype(x.diagonal().rowwise()) {
742  return x.diagonal().rowwise();
743  }
744 
745  template <typename Scalar> typename Eigen::Identity<Scalar>::RowWise rowwise(const Eigen::Identity<Scalar> &x) {
746  return x.rowwise();
747  }
748 
749 
751  template<typename AltMatrix> class RowWise {
752  AltMatrix &alt_matrix;
753  public:
754  typedef typename AltMatrix::Scalar Scalar;
755  typedef typename Eigen::NumTraits< Scalar >::Real Real;
756  typedef Eigen::Matrix<Real,Dynamic,1> RealVector;
757 
758  RowWise(AltMatrix &alt_matrix) : alt_matrix(alt_matrix) {
759  }
760 
761  RealVector squaredNorm() const {
762  if (alt_matrix.isT1()) return rowwise(alt_matrix.t1()).squaredNorm();
763  else return rowwise(alt_matrix.t2()).squaredNorm();
764  }
765  };
766 
767 
768 
769 
770 
771 
772  // ---- Block view of an AltMatrix
773 
774 
775  template<typename AltMatrix> class AltBlock : public AltMatrixBase<AltBlock<AltMatrix>> {
776  public:
777  typedef typename AltMatrix::Scalar Scalar;
778  typedef typename Eigen::NumTraits<Scalar>::Real Real;
779  typedef const AltBlock<AltBlock<AltMatrix>> & Nested;
780 
781  AltBlock(AltMatrix &alt_matrix, Index row, Index col, Index height, Index width) :
782  alt_matrix(alt_matrix), row(row), col(col), height(height), width(width),
783  range(std::max(row, col), std::min(row+width, col+height) - std::max(row, col) + 1)
784  {
785  }
786 
787 
788  bool isT1() const { return alt_matrix.isT1(); }
789 
790  Real squaredNorm() const {
791  if (alt_matrix.isT1())
792  return kqp::blockSquaredNorm(alt_matrix.t1(),row,col,height,width);
793  else
794  return kqp::blockSquaredNorm(alt_matrix.t2(),row,col,height,width);
795  }
796 
797  Index rows() const { return width; }
798  Index cols() const { return height; }
799 
800  // Assignement
801 
802  template<typename Scalar, typename Derived, typename OtherDerived>
803  void assign(const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Derived> &op1,
804  const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Derived> &op2, Index , Index ) {
805  if (op1(0,0) != op2(0,0))
806  KQP_THROW_EXCEPTION_F(not_implemented_exception, "Cannot assign a constant matrix to a constant matrix with different values (%g vs %g)", %op1(0,0) %op2(0,0));
807  }
808 
809 
810  template<typename Op, typename Derived, typename OtherDerived>
811  void assign(const Eigen::CwiseNullaryOp<Op, Derived> &, const OtherDerived &, Index, Index) {
812  KQP_THROW_EXCEPTION_F(not_implemented_exception, "Cannot assign a constant matrix to [%s]", % KQP_DEMANGLE(OtherDerived));
813  }
814 
815  template<typename Scalar, typename OtherDerived>
816  void assign(const Eigen::Identity<Scalar> &, const OtherDerived &, Index, Index) {
817  KQP_THROW_EXCEPTION_F(not_implemented_exception, "Cannot assign an Id matrix to [%s]", % KQP_DEMANGLE(OtherDerived));
818  }
819 
820  template<typename Derived, int Rows, int Cols, typename OtherDerived>
821  void assign(Eigen::Matrix<Derived, Rows, Cols> &mTo, const Eigen::MatrixBase<OtherDerived> &mFrom, Index fromRow, Index fromCol) {
822  mTo.block(row,col,height,width) = mFrom.derived().block(fromRow, fromCol, height, width);
823  }
824 
825  template<typename Scalar, int Rows, int Cols, typename Derived>
826  void assign(Eigen::Matrix<Derived, Rows, Cols> &, const Eigen::Identity<Scalar> &, Index , Index ) {
827  KQP_THROW_EXCEPTION(not_implemented_exception, "Not implemented Matrix = Identity");
828  }
829 
830  template<typename T> friend class AltBlock;
831 
832  template<typename Derived>
833  AltBlock<AltMatrix>& assignTo(const AltBlock<Derived> &from) {
834  if (from.height != this->height || from.width != this->width)
835  KQP_THROW_EXCEPTION_F(out_of_bound_exception, "Block sizes differ in assignement (%d x %d != %d x %d)", %height %width %from.height %from.width);
836  if (alt_matrix.isT1()) {
837  if (from.alt_matrix.isT1()) this->assign(alt_matrix.t1(), from.alt_matrix.t1(), from.row, from.col);
838  else this->assign(alt_matrix.t1(), from.alt_matrix.t2(), from.row, from.col);
839  } else {
840  if (from.alt_matrix.isT1()) this->assign(alt_matrix.t2(), from.alt_matrix.t1(), from.row, from.col);
841  else this->assign(alt_matrix.t2(), from.alt_matrix.t2(), from.row, from.col);
842  }
843 
844  return *this;
845  }
846 
847 
848  private:
849  AltMatrix &alt_matrix;
850  Index row, col, height, width;
851  std::pair<Index, Index> range;
852  public:
853  // Type dependent return
854 
855  auto t1() -> decltype(alt_matrix.m_t1.block(0,0,0,0)) { return alt_matrix.m_t1.block(row, col, height, width); };
856  auto t2() -> decltype(alt_matrix.m_t2.block(0,0,0,0)) { return alt_matrix.m_t2.block(row, col, height, width); };
857  auto t1() const -> decltype(alt_matrix.m_t1.block(0,0,0,0)) { return alt_matrix.m_t1.block(row, col, height, width); };
858  auto t2() const -> decltype(alt_matrix.m_t2.block(0,0,0,0)) { return alt_matrix.m_t2.block(row, col, height, width); };
859 
860 
861  };
862 
863 
864  template <typename Derived, typename OtherDerived>
865  void copy(const Eigen::MatrixBase<Derived> &from, const Eigen::MatrixBase<OtherDerived> &to) {
866  const_cast<Eigen::MatrixBase<OtherDerived>&>(to).derived() = from;
867  }
868 
869 
870  template <typename Derived, typename OtherDerived>
871  void copy(const AltBlock<Derived> &from, const AltBlock<OtherDerived> &to) {
872  const_cast<AltBlock<OtherDerived>&>(to).template assignTo<Derived>(from);
873  }
874 
875 
876 
877  // --- Lazy evaluation
878  template<typename Derived>
879  struct NoAlias {
880  Derived &m;
881 
882  NoAlias(Derived &m) : m(m) {}
883 
884  template<typename OtherDerived>
885  Derived &operator=(const AltMatrixBase<OtherDerived> &expr) {
886  expr.derived().lazyAssign(m);
887  return m;
888  }
889 
890  template<typename OtherDerived>
891  Derived &operator=(const Eigen::MatrixBase<OtherDerived> &expr) {
892  return m.noalias() = expr.derived();
893  }
894 
895  };
896 
897  template<typename Derived>
898  NoAlias<Derived> noalias(Eigen::MatrixBase<Derived> &m) {
899  return NoAlias<Derived>(m.derived());
900  }
901 
902 
903  // ---
904  // --- Expressions
905  // ----
906 
907 
908  // --- Forward declarations
909 
910  template<typename Op, typename Lhs, typename Rhs, int ProductOrder> class AltMatrixOp;
911  template<typename Op, typename Lhs, typename Rhs, int ProductOrder> class AltArrayOp;
912 
914  template<typename Op, typename Lhs, typename Rhs, int Side, bool isT1> struct ExprType;
915 
916 
918  struct MultOp {
919  template<typename T1, typename T2> static inline auto apply(const T1 &t1, const T2 &t2) -> decltype(t1 * t2) { return t1 * t2; }
920  };
922  struct MinusOp {
923  template<typename T1, typename T2> static inline auto apply(const T1 &t1, const T2 &t2) -> decltype(t1 - t2) { return t1 - t2; }
924  };
925 
926 
928  template<typename Op, typename _Lhs, typename _Rhs, bool isT1> struct ExprType<Op, _Lhs, _Rhs, Eigen::OnTheLeft, isT1> {
929  const _Lhs &_lhs;
930  typedef decltype(_lhs.t1()) LhsT1;
931  typedef decltype(_lhs.t2()) LhsT2;
932  typedef typename Eigen::internal::conditional<isT1, LhsT1, LhsT2>::type Lhs;
933  const Lhs &lhs;
934 
935  typedef _Rhs Rhs;
936  const Rhs &rhs;
937 
938  typedef decltype(Op::apply(lhs, rhs)) Type;
939  };
940 
942  template<typename Op, typename _Lhs, typename _Rhs, bool isT1> struct ExprType<Op, _Lhs, _Rhs, Eigen::OnTheRight, isT1> {
943  typedef _Lhs Lhs;
944  const Lhs &lhs;
945 
946  const _Rhs &_rhs;
947  typedef decltype(_rhs.t1()) RhsT1;
948  typedef decltype(_rhs.t2()) RhsT2;
949  typedef typename Eigen::internal::conditional<isT1, RhsT1, RhsT2>::type Rhs;
950  const Rhs &rhs;
951 
952  typedef decltype(Op::apply(lhs, rhs)) Type;
953  };
954 
955 
957  template<typename Op, typename _Lhs, typename _Rhs, int Side, bool isT1>
958  struct AltExpression {
960 
961  typedef typename Types::Type Expression;
962  typedef typename Types::Lhs Lhs;
963  typedef typename Types::Rhs Rhs;
964 
965  Expression expression;
966  AltExpression (const Lhs &lhs, const Rhs &rhs) : expression(Op::apply(lhs, rhs)) {}
967  };
968 
969 
970 
971  // --- Array wrapper
972 
973  template<typename Derived>
974  class AltArrayBase {
975  public:
976  const Derived& derived() const { return static_cast<const Derived&>(*this); }
977  typename Eigen::internal::traits<Derived>::MatrixExpr matrix() const { return derived().matrix(); }
978  };
979 
980  template<typename Derived>
981  class AltArrayWrapper: public AltArrayBase<AltArrayWrapper<Derived>> {
982  typename Eigen::internal::ref_selector<Derived>::type m_alt;
983  public:
984  AltArrayWrapper(const Derived &alt) : m_alt(alt) {}
985  bool isT1() const { return m_alt.isT1(); }
986  auto t1() const -> decltype(m_alt.t1().array()) { return m_alt.t1().array(); }
987  auto t2() const -> decltype(m_alt.t2().array()) { return m_alt.t2().array(); }
988 
989  auto matrix() const -> decltype(m_alt) { return m_alt; }
990 
991  Index rows() const {
992  return m_alt.rows();
993  }
994  Index cols() const { return m_alt.cols(); }
995 
996  };
997 
999  template<typename Derived>
1000  class AltMatrixWrapper: public AltMatrixBase<AltMatrixWrapper<Derived>> {
1001  typename Eigen::internal::ref_selector<Derived>::type m_expr;
1002  public:
1003  AltMatrixWrapper(const Derived &expr) : m_expr(expr) {}
1004  Index rows() const {
1005  return m_expr.rows();
1006  }
1007  Index cols() const {
1008  return m_expr.cols();
1009  }
1010 
1011  bool isT1() const { return m_expr.isT1(); }
1012  auto t1() const -> decltype(m_expr.t1().matrix()) { return m_expr.t1().matrix(); }
1013  auto t2() const -> decltype(m_expr.t2().matrix()) { return m_expr.t2().matrix(); }
1014 
1015  template<typename Dest> void evalTo(Dest& dest) const {
1016  if (isT1())
1017  dest = t1();
1018  else
1019  dest = t2();
1020  }
1021 
1022  };
1023 
1024 #define KQP_ALT_OP_LEFT(op, opname, Result, AltType, OtherType)\
1025  template<typename Derived, typename OtherDerived>\
1026  Result<opname, Derived, OtherDerived, Eigen::OnTheLeft>\
1027  operator op(const AltType<Derived> &lhs, const OtherType<OtherDerived> &rhs) {\
1028  return Result<opname, Derived, OtherDerived, Eigen::OnTheLeft>(lhs.derived(), rhs.derived());\
1029  };
1030 #define KQP_ALT_OP_RIGHT(op, opname, Result, AltType, OtherType)\
1031  template<typename Derived, typename OtherDerived>\
1032  Result<opname, Derived, OtherDerived, Eigen::OnTheRight>\
1033  operator op(const OtherType<Derived> &lhs, const AltType<OtherDerived> &rhs) {\
1034  return Result<opname, Derived, OtherDerived, Eigen::OnTheRight>(lhs.derived(), rhs.derived());\
1035  };
1036 
1037 #define KQP_ALT_OP(op, opname, Result, AltType, OtherType)\
1038  KQP_ALT_OP_LEFT(op, opname, Result, AltType, OtherType)\
1039  KQP_ALT_OP_RIGHT(op, opname, Result, AltType, OtherType)\
1040  KQP_ALT_OP_LEFT(op, opname, Result, AltType, AltType)
1041 
1042 KQP_ALT_OP(-, MinusOp, AltArrayOp, AltArrayBase, Eigen::ArrayBase)
1043 KQP_ALT_OP(*, MultOp, AltArrayOp, AltArrayBase, Eigen::ArrayBase)
1044 
1045  template<typename Op, typename Lhs, typename Rhs, int Side>
1046  class AltArrayOp : Eigen::internal::no_assignment_operator, public AltArrayBase< AltArrayOp<Op, Lhs,Rhs,Side> >
1047  {
1050 
1051  typedef typename HolderT1::Expression ExprIfT1;
1052  typedef typename HolderT2::Expression ExprIfT2;
1053 
1054  boost::shared_ptr<HolderT1> _t1;
1055  boost::shared_ptr<HolderT2> _t2;
1056 
1057  bool m_isT1;
1058 
1059  public:
1060  typedef AltArrayOp Nested;
1061 
1062  typedef typename kqp::scalar<Lhs>::type Scalar;
1063  typedef typename Eigen::NumTraits<Scalar>::Real Real;
1064 
1065  // Initialisation when the Alt is on the left
1066  friend void initAltArrayOp(AltArrayOp<Op, Lhs, Rhs, Eigen::OnTheLeft> &op, const Lhs &lhs, const Rhs &rhs) {
1069  if ((op.m_isT1 = lhs.isT1())) op._t1.reset(new HolderT1(lhs.t1(), rhs));
1070  else op._t2.reset(new HolderT2(lhs.t2(), rhs));
1071  }
1072 
1073  // Initialisation when the Alt is on the right
1074  friend void initAltArrayOp(AltArrayOp<Op, Lhs, Rhs, Eigen::OnTheRight> &op, const Lhs &lhs, const Rhs &rhs) {
1077  if ((op.m_isT1 = rhs.isT1())) op._t1.reset(new HolderT1(lhs, rhs.t1()));
1078  else op._t2.reset(new HolderT2(lhs, rhs.t2()));
1079  }
1080 
1081 
1082  AltArrayOp(const Lhs& lhs, const Rhs& rhs) {
1083  initAltArrayOp(*this, lhs, rhs);
1084  }
1085 
1086  bool isT1() const { return m_isT1; }
1087 
1088  inline const ExprIfT1 & t1() const { return _t1->expression; }
1089  inline const ExprIfT2 & t2() const { return _t2->expression; }
1090 
1091  AltMatrixWrapper<AltArrayOp> matrix() const { return AltMatrixWrapper<AltArrayOp>(*this); }
1092 
1093  Index rows() const {
1094  return m_isT1 ? _t1->expression.rows() : _t2->expression.rows();
1095  }
1096  Index cols() const {return m_isT1 ? _t1->expression.cols() : _t2->expression.cols(); }
1097  };
1098 
1099  // --- Multiplication between an Alt matrix and anything
1100 
1101  template<typename Op, typename Lhs, typename Rhs, int Side>
1102  class AltMatrixOp : Eigen::internal::no_assignment_operator, public AltMatrixBase< AltMatrixOp<Op, Lhs,Rhs,Side> >
1103  {
1106 
1107  typedef typename HolderT1::Expression ExprIfT1;
1108  typedef typename HolderT2::Expression ExprIfT2;
1109 
1110  boost::shared_ptr<HolderT1> _t1;
1111  boost::shared_ptr<HolderT2> _t2;
1112 
1113  bool m_isT1;
1114 
1115  public:
1116  typedef typename kqp::scalar<Lhs>::type Scalar;
1117  typedef typename Eigen::NumTraits<Scalar>::Real Real;
1118 
1119  // Initialisation when the Alt is on the left
1120  friend void initAltMatrixOp(AltMatrixOp<Op, Lhs, Rhs, Eigen::OnTheLeft> &op, const Lhs &lhs, const Rhs &rhs) {
1123  if ((op.m_isT1 = lhs.isT1())) op._t1.reset(new HolderT1(lhs.t1(), rhs));
1124  else op._t2.reset(new HolderT2(lhs.t2(), rhs));
1125  }
1126 
1127  // Initialisation when the Alt is on the right
1128  friend void initAltMatrixOp(AltMatrixOp<Op, Lhs, Rhs, Eigen::OnTheRight> &op, const Lhs &lhs, const Rhs &rhs) {
1131  if ((op.m_isT1 = rhs.isT1())) op._t1.reset(new HolderT1(lhs, rhs.t1()));
1132  else op._t2.reset(new HolderT2(lhs, rhs.t2()));
1133  }
1134 
1135 
1136  AltMatrixOp(const Lhs& lhs, const Rhs& rhs) {
1137  initAltMatrixOp(*this, lhs, rhs);
1138  }
1139 
1140 
1141  Index rows() const { return m_isT1 ? _t1->expression.rows() : _t2->expression.rows(); }
1142  Index cols() const {return m_isT1 ? _t1->expression.cols() : _t2->expression.cols(); }
1143 
1144  Real squaredNorm() const { return m_isT1 ? _t1->expression.squaredNorm() : _t2->expression.squaredNorm(); }
1145  Scalar trace() const { return m_isT1 ? _t1->expression.trace() : _t2->expression.trace(); }
1146 
1147 
1148  bool isT1() const { return m_isT1; }
1149 
1150 
1151  inline const ExprIfT1 & t1() const { return _t1->expression; }
1152  inline const ExprIfT2 & t2() const { return _t2->expression; }
1153 
1154  template<typename Dest> void evalTo(Dest& dest) const {
1155  if (m_isT1)
1156  dest = t1();
1157  else
1158  dest = t2();
1159  }
1160 
1161 
1162  template<typename Dest> void lazyAssign(Dest& dest) const {
1163  if (m_isT1)
1164  noalias(dest) = t1();
1165  else
1166  noalias(dest) = t2();
1167  }
1168 
1169 
1170  void printExpression(std::ostream &out) const {
1171  out << "altx[";
1172  if (m_isT1)
1173  kqp::printExpression(out, t1());
1174  else
1175  kqp::printExpression(out, t2());
1176  out << "]";
1177  }
1178  };
1179 
1180 
1181 
1182  KQP_ALT_OP(*, MultOp, kqp::AltMatrixOp, kqp::AltMatrixBase, Eigen::EigenBase);
1183 
1184  // --- Scalar * AltMatrix
1185  template<class Derived>
1187  operator*(typename Eigen::internal::traits<Derived>::Scalar alpha, const kqp::AltMatrixBase<Derived> &a) {
1188  return AltMatrixOp<MultOp,typename Eigen::internal::traits<Derived>::Scalar, Derived, Eigen::OnTheRight>(alpha, a.derived());
1189  }
1190 
1191 
1192 
1193 
1194  // ---- Diagonal Wrapper
1195  template<typename Derived> struct DiagonalBlockWrapper {
1196  typename Derived::Nested m_value;
1197  Index startRow, startCol;
1198  Index rows, cols;
1199  DiagonalBlockWrapper(const Derived &value, Index startRow, Index startCol, Index blockRows, Index blockCols)
1200  : m_value(value), startRow(startRow), startCol(startCol), rows(blockRows), cols(blockCols) {}
1201  const Derived &derived() const { return m_value; }
1202  };
1203 
1204 
1205 
1206  template<typename Lhs, typename Rhs, int side> struct DiagonalBlockWrapperDenseMult;
1207 
1208  template<typename Lhs, typename Rhs> struct DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheRight>
1209  : public Eigen::MatrixBase<DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheRight>> {
1211  typename Rhs::Nested m_rhs;
1212  typedef typename Eigen::internal::traits<DiagonalBlockWrapperDenseMult>::Scalar Scalar;
1213 
1214  Index zerosAbove, zerosLeft, first, size;
1215 
1216  DiagonalBlockWrapperDenseMult(const DiagonalBlockWrapper<Lhs>& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) {
1217  zerosAbove = std::max(0l, m_lhs.startCol - m_lhs.startRow);
1218  zerosLeft = std::max(0l, m_lhs.startRow - m_lhs.startCol);
1219  // First index of the diagonal
1220  first = std::max(m_lhs.startCol, m_lhs.startRow);
1221  // Number of values
1222  size = std::min(m_lhs.rows - zerosAbove, m_lhs.cols - zerosLeft);
1223  }
1224 
1225 
1226  Scalar coeff(Index row, Index col) const {
1227  if (row < zerosAbove || row >= zerosAbove + size) return 0;
1228  return m_lhs.m_value[first+row-zerosAbove] * m_rhs(row+zerosLeft-zerosAbove, col);
1229 
1230  }
1231 
1232  template<typename Dest> inline void evalTo(Dest&) const {
1233  KQP_THROW_EXCEPTION(not_implemented_exception, "evalTo");
1234  }
1235 
1236 
1237  Index rows() const { return m_lhs.rows; }
1238  Index cols() const { return m_rhs.cols(); }
1239  };
1240 
1241  template<typename Lhs, typename Rhs> struct DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheLeft>
1242  : public Eigen::MatrixBase<DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheLeft>> {
1243  typename Lhs::Nested m_lhs;
1245  typedef typename Eigen::internal::traits<DiagonalBlockWrapperDenseMult>::Scalar Scalar;
1246 
1247  Index zerosAbove, zerosLeft, first, size;
1248 
1249  DiagonalBlockWrapperDenseMult(const Lhs& lhs, const DiagonalBlockWrapper<Rhs>& rhs) : m_lhs(lhs), m_rhs(rhs) {
1250  zerosAbove = std::max(0l, m_rhs.startCol - m_rhs.startRow);
1251  zerosLeft = std::max(0l, m_rhs.startRow - m_rhs.startCol);
1252  // First index of the diagonal
1253  first = std::max(m_rhs.startCol, m_rhs.startRow);
1254  // Number of values
1255  size = std::min(m_rhs.rows - zerosAbove, m_rhs.cols - zerosLeft);
1256  }
1257 
1258 
1259  Scalar coeff(Index row, Index col) const {
1260  if (col < zerosLeft || col >= zerosLeft + size) return 0;
1261  return m_lhs(row, col+zerosAbove-zerosLeft) * m_rhs.m_value[first+col-zerosLeft];
1262 
1263  }
1264 
1265  template<typename Dest> inline void evalTo(Dest&) const {
1266  KQP_THROW_EXCEPTION(not_implemented_exception, "evalTo");
1267  }
1268 
1269 
1270  Index rows() const { return m_lhs.rows(); }
1271  Index cols() const { return m_rhs.cols; }
1272  };
1273 
1274 
1275  // --- DiagonalBlockWrapper * Dense
1276  template<class Lhs,class Rhs>
1277  DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheRight> operator*(const kqp::DiagonalBlockWrapper<Lhs> &a, const Eigen::MatrixBase<Rhs> &b)
1278  {
1280  }
1281 
1282  template<class Lhs,class Rhs>
1283  DiagonalBlockWrapperDenseMult<Lhs,Rhs,Eigen::OnTheLeft> operator*(const Eigen::MatrixBase<Lhs> &a, const kqp::DiagonalBlockWrapper<Rhs> &b)
1284  {
1286  }
1287 
1288 
1289 }
1290 
1291 namespace Eigen {
1292 
1293  // --- Multiplication of two diagonal wrappers
1294  template<typename Derived, typename OtherDerived>
1295  auto operator*(const DiagonalWrapper<Derived> &a, const DiagonalWrapper<OtherDerived> &b)
1296  -> decltype(a.diagonal().cwiseProduct(b.diagonal()).asDiagonal()) {
1297  return a.diagonal().cwiseProduct(b.diagonal()).asDiagonal();
1298  }
1299 
1300  // --- Multiplication scalar * diagonal
1301  template<typename Derived>
1302  auto operator*(typename Eigen::internal::traits<Derived>::Scalar alpha, const DiagonalWrapper<Derived> &a) -> decltype((alpha * a.diagonal()).asDiagonal()) {
1303  return (alpha * a.diagonal()).asDiagonal();
1304  }
1305 
1306 
1307  // --- Traits
1308  namespace internal {
1309  // AlMatrixBase
1310  template<typename Derived>
1311  struct traits< kqp::AltMatrixBase<Derived> > {
1313  typedef typename MatrixXd::Index Index;
1314  typedef typename traits<Derived>::Scalar Scalar;
1315 
1316  enum {
1317  Flags = 0,
1318  RowsAtCompileTime = Dynamic,
1319  ColsAtCompileTime = Dynamic
1320  };
1321 
1322  };
1323 
1324  // AltMatrix
1325  template<typename T1, typename T2>
1326  struct traits< kqp::AltMatrix<T1,T2> > {
1328  typedef typename MatrixXd::Index Index;
1329  typedef typename traits<T1>::Scalar Scalar;
1330  enum {
1331  Flags = NestByRefBit,
1332  RowsAtCompileTime = Dynamic,
1333  ColsAtCompileTime = Dynamic
1334  };
1335  };
1336 
1337  // Adjoint of AltMatrix
1338  template<typename Derived>
1339  struct traits< kqp::Adjoint< kqp::AltMatrixBase<Derived> > > {
1341  typedef typename MatrixXd::Index Index;
1342  typedef typename traits<Derived>::Scalar Scalar;
1343  enum {
1344  Flags = 0,
1345  RowsAtCompileTime = Dynamic,
1346  ColsAtCompileTime = Dynamic
1347  };
1348  };
1349 
1350  template<typename BaseT1, typename BaseT2>
1351  struct traits< kqp::Adjoint< kqp::AltMatrix<BaseT1, BaseT2> > > {
1353  typedef typename MatrixXd::Index Index;
1354  typedef typename traits< kqp::AltMatrix<BaseT1, BaseT2> >::Scalar Scalar;
1355  enum {
1356  Flags = 0,
1357  RowsAtCompileTime = Dynamic,
1358  ColsAtCompileTime = Dynamic
1359  };
1360  };
1361  template<typename UnaryOp, typename Derived>
1362  struct traits<kqp::AltCwiseUnaryOp<UnaryOp, Derived>> {
1364  typedef typename MatrixXd::Index Index;
1365  typedef typename traits<Derived>::Scalar Scalar;
1366  enum {
1367  Flags = 0,
1368  RowsAtCompileTime = Dynamic,
1369  ColsAtCompileTime = Dynamic
1370  };
1371  };
1372 
1373  template<typename Derived>
1374  struct traits<kqp::AltAsDiagonal<Derived>> {
1376  typedef typename MatrixXd::Index Index;
1377  typedef typename traits<Derived>::Scalar Scalar;
1378  enum {
1379  Flags = 0,
1380  RowsAtCompileTime = Dynamic,
1381  ColsAtCompileTime = Dynamic
1382  };
1383  };
1384 
1385  template<typename AltMatrix>
1386  struct traits<kqp::AltBlock<AltMatrix>> {
1388  typedef typename MatrixXd::Index Index;
1389  typedef typename traits<AltMatrix>::Scalar Scalar;
1390  enum {
1391  Flags = 0,
1392  RowsAtCompileTime = Dynamic,
1393  ColsAtCompileTime = Dynamic
1394  };
1395  };
1396 
1397  template<typename Derived>
1398  struct traits<kqp::AltArrayWrapper<Derived>> {
1399  typedef typename traits<Derived>::Scalar Scalar;
1400  typedef Derived MatrixExpr;
1401  };
1402 
1403  template<typename Op, typename Lhs, typename Rhs, int Side>
1404  struct traits<kqp::AltArrayOp<Op, Lhs, Rhs, Side>> {
1405  typedef typename scalar_product_traits<typename kqp::scalar<Lhs>::type, typename kqp::scalar<Rhs>::type>::ReturnType Scalar;
1406  typedef MatrixWrapper<kqp::AltArrayOp<Op, Lhs, Rhs, Side>> MatrixExpr;
1407 
1408  typedef kqp::Index Index;
1410 
1411  enum {
1412  Flags = 0,
1413  };
1414 
1415  };
1416 
1417 
1418  template<typename Derived>
1419  struct traits<kqp::AltMatrixWrapper<Derived>> {
1420  typedef typename traits<Derived>::Scalar Scalar;
1422  typedef typename traits<Derived>::Index Index;
1423  };
1424 
1425 
1426  template<typename Lhs, typename Rhs, int Side>
1427  struct traits<kqp::DiagonalBlockWrapperDenseMult<Lhs,Rhs,Side>> {
1428  typedef Dense StorageKind;
1429  typedef typename MatrixXd::Index Index;
1430  typedef typename scalar_product_traits<typename kqp::scalar<Lhs>::type, typename kqp::scalar<Rhs>::type>::ReturnType Scalar;
1431  typedef MatrixXpr XprKind;
1432  enum {
1433  Flags = 0,
1434  RowsAtCompileTime = Dynamic,
1435  ColsAtCompileTime = Dynamic,
1436  MaxRowsAtCompileTime = Dynamic,
1437  MaxColsAtCompileTime = Dynamic,
1438  CoeffReadCost = 1,
1439  };
1440  };
1441 
1442 
1443  // Alt * Eigen
1444  template<typename Op, typename Lhs, typename Rhs, int Side>
1445  struct traits< kqp::AltMatrixOp<Op,Lhs, Rhs, Side> > {
1447  typedef typename MatrixXd::Index Index;
1448 
1449  typedef typename scalar_product_traits<typename kqp::scalar<Lhs>::type, typename kqp::scalar<Rhs>::type>::ReturnType Scalar;
1450  enum {
1451  Flags = 0,
1452  RowsAtCompileTime = Dynamic,
1453  ColsAtCompileTime = Dynamic
1454  };
1455  };
1456 
1457 
1458  }
1459 
1460 
1461 }
1462 
1463 #endif