New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
View | Details | Raw Unified | Return to bug 838
Collapse All | Expand All

(-)a/Eigen/src/SparseCore/SparseDenseProduct.h (-26 / +75 lines)
Lines 14-139 namespace Eigen { Link Here
14
14
15
template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductReturnType
15
template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductReturnType
16
{
16
{
17
  typedef SparseTimeDenseProduct<Lhs,Rhs> Type;
17
  typedef SparseTimeDenseProduct<Lhs,Rhs> Type;
18
};
18
};
19
19
20
template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
20
template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
21
{
21
{
22
  typedef SparseDenseOuterProduct<Lhs,Rhs,false> Type;
22
  typedef SparseDenseOuterProduct<Lhs,Rhs,true> Type;
23
};
23
};
24
24
25
template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
25
template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
26
{
26
{
27
  typedef DenseTimeSparseProduct<Lhs,Rhs> Type;
27
  typedef DenseTimeSparseProduct<Lhs,Rhs> Type;
28
};
28
};
29
29
30
template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
30
template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
31
{
31
{
32
  typedef SparseDenseOuterProduct<Rhs,Lhs,true> Type;
32
  typedef SparseDenseOuterProduct<Rhs,Lhs,false> Type;
33
};
33
};
34
34
35
namespace internal {
35
namespace internal {
36
36
37
template<typename Lhs, typename Rhs, bool Tr>
37
template<typename Lhs, typename Rhs, bool SparseFirst>
38
struct traits<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
38
struct traits<SparseDenseOuterProduct<Lhs,Rhs,SparseFirst> >
39
{
39
{
40
  typedef Sparse StorageKind;
40
  typedef Sparse StorageKind;
41
  typedef typename scalar_product_traits<typename traits<Lhs>::Scalar,
41
  typedef typename scalar_product_traits<typename traits<Lhs>::Scalar,
42
                                         typename traits<Rhs>::Scalar>::ReturnType Scalar;
42
                                         typename traits<Rhs>::Scalar>::ReturnType Scalar;
43
  typedef typename Lhs::Index Index;
43
  typedef typename Lhs::Index Index; 
44
  typedef typename Lhs::Nested LhsNested;
44
  typedef typename Lhs::Nested LhsNested;
45
  typedef typename Rhs::Nested RhsNested;
45
  typedef typename Rhs::Nested RhsNested;
46
46
  typedef typename remove_all<LhsNested>::type _LhsNested;
47
  typedef typename remove_all<LhsNested>::type _LhsNested;
47
  typedef typename remove_all<RhsNested>::type _RhsNested;
48
  typedef typename remove_all<RhsNested>::type _RhsNested;
48
49
49
  enum {
50
  enum {
50
    LhsCoeffReadCost = traits<_LhsNested>::CoeffReadCost,
51
    LhsCoeffReadCost = traits<_LhsNested>::CoeffReadCost,
51
    RhsCoeffReadCost = traits<_RhsNested>::CoeffReadCost,
52
    RhsCoeffReadCost = traits<_RhsNested>::CoeffReadCost,
52
53
53
    RowsAtCompileTime    = Tr ? int(traits<Rhs>::RowsAtCompileTime)     : int(traits<Lhs>::RowsAtCompileTime),
54
    RowsAtCompileTime    = SparseFirst ? int(traits<Rhs>::RowsAtCompileTime)     : int(traits<Lhs>::RowsAtCompileTime),
54
    ColsAtCompileTime    = Tr ? int(traits<Lhs>::ColsAtCompileTime)     : int(traits<Rhs>::ColsAtCompileTime),
55
    ColsAtCompileTime    = SparseFirst ? int(traits<Lhs>::ColsAtCompileTime)     : int(traits<Rhs>::ColsAtCompileTime),
55
    MaxRowsAtCompileTime = Tr ? int(traits<Rhs>::MaxRowsAtCompileTime)  : int(traits<Lhs>::MaxRowsAtCompileTime),
56
    MaxRowsAtCompileTime = SparseFirst ? int(traits<Rhs>::MaxRowsAtCompileTime)  : int(traits<Lhs>::MaxRowsAtCompileTime),
56
    MaxColsAtCompileTime = Tr ? int(traits<Lhs>::MaxColsAtCompileTime)  : int(traits<Rhs>::MaxColsAtCompileTime),
57
    MaxColsAtCompileTime = SparseFirst ? int(traits<Lhs>::MaxColsAtCompileTime)  : int(traits<Rhs>::MaxColsAtCompileTime),
57
58
58
    Flags = Tr ? RowMajorBit : 0,
59
    Flags = traits<Lhs>::Flags & RowMajorBit,
59
60
60
    CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + NumTraits<Scalar>::MulCost
61
    CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + NumTraits<Scalar>::MulCost
61
  };
62
  };
62
};
63
};
63
64
64
} // end namespace internal
65
} // end namespace internal
65
66
66
template<typename Lhs, typename Rhs, bool Tr>
67
68
template<typename Lhs, typename Rhs, bool UseSparseIterator, bool SparseFirst>
69
class sdop_inner_iterator;
70
71
template<typename Lhs, typename Rhs, bool SparseFirst>
67
class SparseDenseOuterProduct
72
class SparseDenseOuterProduct
68
 : public SparseMatrixBase<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
73
 : public SparseMatrixBase<SparseDenseOuterProduct<Lhs,Rhs,SparseFirst> >
69
{
74
{
70
  public:
75
  public:
71
76
72
    typedef SparseMatrixBase<SparseDenseOuterProduct> Base;
77
    typedef SparseMatrixBase<SparseDenseOuterProduct> Base;
73
    EIGEN_DENSE_PUBLIC_INTERFACE(SparseDenseOuterProduct)
78
    EIGEN_DENSE_PUBLIC_INTERFACE(SparseDenseOuterProduct)
74
    typedef internal::traits<SparseDenseOuterProduct> Traits;
79
    typedef internal::traits<SparseDenseOuterProduct> Traits;
75
80
76
  private:
81
  private:
77
82
78
    typedef typename Traits::LhsNested LhsNested;
83
    typedef typename Traits::LhsNested LhsNested;
79
    typedef typename Traits::RhsNested RhsNested;
84
    typedef typename Traits::RhsNested RhsNested;
80
    typedef typename Traits::_LhsNested _LhsNested;
85
    typedef typename Traits::_LhsNested _LhsNested;
81
    typedef typename Traits::_RhsNested _RhsNested;
86
    typedef typename Traits::_RhsNested _RhsNested;
82
87
83
  public:
88
  public:
84
89
85
    class InnerIterator;
90
    typedef sdop_inner_iterator<Lhs, Rhs, SparseFirst == ((Traits::Flags&RowMajorBit)==0), SparseFirst> InnerIterator;
86
91
87
    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Lhs& lhs, const Rhs& rhs)
92
    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Lhs& lhs, const Rhs& rhs)
88
      : m_lhs(lhs), m_rhs(rhs)
93
      : m_lhs(lhs), m_rhs(rhs)
89
    {
94
    {
90
      EIGEN_STATIC_ASSERT(!Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
95
      EIGEN_STATIC_ASSERT(SparseFirst,YOU_MADE_A_PROGRAMMING_MISTAKE);
91
    }
96
    }
92
97
93
    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Rhs& rhs, const Lhs& lhs)
98
    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Rhs& rhs, const Lhs& lhs)
94
      : m_lhs(lhs), m_rhs(rhs)
99
      : m_lhs(lhs), m_rhs(rhs)
95
    {
100
    {
96
      EIGEN_STATIC_ASSERT(Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
101
      EIGEN_STATIC_ASSERT(!SparseFirst,YOU_MADE_A_PROGRAMMING_MISTAKE);
97
    }
102
    }
98
103
99
    EIGEN_STRONG_INLINE Index rows() const { return Tr ? m_rhs.rows() : m_lhs.rows(); }
104
    EIGEN_STRONG_INLINE Index rows() const { return SparseFirst ? m_lhs.rows() : m_rhs.rows(); }
100
    EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : m_rhs.cols(); }
105
    EIGEN_STRONG_INLINE Index cols() const { return SparseFirst ? m_rhs.cols() : m_lhs.cols(); }
101
106
    
102
    EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
107
    EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
103
    EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
108
    EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
104
109
105
  protected:
110
  protected:
106
    LhsNested m_lhs;
111
    LhsNested m_lhs;
107
    RhsNested m_rhs;
112
    RhsNested m_rhs;
108
};
113
};
109
114
110
template<typename Lhs, typename Rhs, bool Transpose>
115
111
class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNested::InnerIterator
116
template<typename Lhs, typename Rhs, bool SparseFirst>
117
class sdop_inner_iterator<Lhs,Rhs,true,SparseFirst>
118
  : public SparseDenseOuterProduct<Lhs, Rhs, SparseFirst>::Traits::_LhsNested::InnerIterator
112
{
119
{
113
    typedef typename _LhsNested::InnerIterator Base;
120
    typedef SparseDenseOuterProduct<Lhs, Rhs, SparseFirst> sdop;
114
    typedef typename SparseDenseOuterProduct::Index Index;
121
    enum { ColMajor = ((sdop::Traits::Flags&RowMajorBit)==0) };
122
    typedef typename sdop::Traits::_LhsNested::InnerIterator Base;
123
    typedef typename sdop::Index Index;
124
    typedef typename sdop::Scalar Scalar;
125
115
  public:
126
  public:
116
    EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer)
127
    EIGEN_STRONG_INLINE sdop_inner_iterator(const sdop& prod, Index outer) 
117
      : Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer))
128
      : Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer))
118
    {
129
    {
119
    }
130
    }
120
131
    
121
    inline Index outer() const { return m_outer; }
132
    inline Index outer() const { return m_outer; }
122
    inline Index row() const { return Transpose ? Base::row() : m_outer; }
133
    inline Index row() const { return ColMajor ? Base::index() : m_outer; }
123
    inline Index col() const { return Transpose ? m_outer : Base::row(); }
134
    inline Index col() const { return ColMajor ? m_outer : Base::index(); }
124
135
125
    inline Scalar value() const { return Base::value() * m_factor; }
136
    inline Scalar value() const { return Base::value() * m_factor; }
126
137
127
  protected:
138
  protected:
128
    Index m_outer;
139
    Index m_outer;
129
    Scalar m_factor;
140
    Scalar m_factor;
130
};
141
};
131
142
143
template<typename Lhs, typename Rhs, bool SparseFirst>
144
class sdop_inner_iterator<Lhs,Rhs,false,SparseFirst> 
145
{
146
    typedef SparseDenseOuterProduct<Lhs, Rhs, SparseFirst> sdop;
147
    enum { ColMajor = ((sdop::Traits::Flags&RowMajorBit)==0) };
148
    typedef typename sdop::Index Index;
149
    typedef typename sdop::Scalar Scalar;
150
151
  public:
152
    EIGEN_STRONG_INLINE sdop_inner_iterator(const sdop& prod, Index outer) 
153
      : m_outer(outer), m_prod(prod), m_id(0), m_end(0)
154
    {
155
      typename sdop::Traits::_LhsNested::InnerIterator i(prod.lhs(), outer);
156
      if (i) {
157
	m_end = SparseFirst ? prod.rhs().cols() : prod.rhs().rows();
158
	m_factor = i.value();
159
      }
160
    }
161
  
162
    inline sdop_inner_iterator& operator++() { m_id++; return *this; }
163
   
164
    inline Scalar value() const { return m_factor * m_prod.rhs().coeff(m_id); }
165
    
166
    inline Index index() const { return m_id; }
167
    inline Index outer() const { return m_outer; }
168
    inline Index row() const { return ColMajor ? m_id : m_outer; }
169
    inline Index col() const { return ColMajor ? m_outer : m_id; }
170
    
171
    inline operator bool() const { return (m_id < m_end); }
172
   
173
  protected:
174
    const Index m_outer;
175
    const sdop& m_prod;
176
    Index m_id;
177
    Index m_end;
178
    Scalar m_factor;
179
};
180
132
namespace internal {
181
namespace internal {
133
template<typename Lhs, typename Rhs>
182
template<typename Lhs, typename Rhs>
134
struct traits<SparseTimeDenseProduct<Lhs,Rhs> >
183
struct traits<SparseTimeDenseProduct<Lhs,Rhs> >
135
 : traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
184
 : traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
136
{
185
{
137
  typedef Dense StorageKind;
186
  typedef Dense StorageKind;
138
  typedef MatrixXpr XprKind;
187
  typedef MatrixXpr XprKind;
139
};
188
};
(-)a/Eigen/src/SparseCore/SparseUtil.h (-2 / +2 lines)
Lines 76-96 template<typename _Scalar, int _Flags = Link Here
76
template<typename MatrixType, int Mode>           class SparseTriangularView;
76
template<typename MatrixType, int Mode>           class SparseTriangularView;
77
template<typename MatrixType, unsigned int UpLo>  class SparseSelfAdjointView;
77
template<typename MatrixType, unsigned int UpLo>  class SparseSelfAdjointView;
78
template<typename Lhs, typename Rhs>              class SparseDiagonalProduct;
78
template<typename Lhs, typename Rhs>              class SparseDiagonalProduct;
79
template<typename MatrixType> class SparseView;
79
template<typename MatrixType> class SparseView;
80
80
81
template<typename Lhs, typename Rhs>        class SparseSparseProduct;
81
template<typename Lhs, typename Rhs>        class SparseSparseProduct;
82
template<typename Lhs, typename Rhs>        class SparseTimeDenseProduct;
82
template<typename Lhs, typename Rhs>        class SparseTimeDenseProduct;
83
template<typename Lhs, typename Rhs>        class DenseTimeSparseProduct;
83
template<typename Lhs, typename Rhs>        class DenseTimeSparseProduct;
84
template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct;
84
template<typename Lhs, typename Rhs, bool SparseFirst> class SparseDenseOuterProduct;
85
85
86
template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType;
86
template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType;
87
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType;
87
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType;
88
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct SparseDenseProductReturnType;
88
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Rhs>::RowsAtCompileTime> struct SparseDenseProductReturnType;
89
template<typename MatrixType,int UpLo> class SparseSymmetricPermutationProduct;
89
template<typename MatrixType,int UpLo> class SparseSymmetricPermutationProduct;
90
90
91
namespace internal {
91
namespace internal {
92
92
93
template<typename T,int Rows,int Cols> struct sparse_eval;
93
template<typename T,int Rows,int Cols> struct sparse_eval;
94
94
95
template<typename T> struct eval<T,Sparse>
95
template<typename T> struct eval<T,Sparse>
96
  : public sparse_eval<T, traits<T>::RowsAtCompileTime,traits<T>::ColsAtCompileTime>
96
  : public sparse_eval<T, traits<T>::RowsAtCompileTime,traits<T>::ColsAtCompileTime>
(-)a/test/CMakeLists.txt (+2 lines)
Lines 246-261 ei_add_test(rvalue_types) Link Here
246
ei_add_test(dense_storage)
246
ei_add_test(dense_storage)
247
247
248
ei_add_test(simplicial_cholesky)
248
ei_add_test(simplicial_cholesky)
249
ei_add_test(conjugate_gradient)
249
ei_add_test(conjugate_gradient)
250
ei_add_test(bicgstab)
250
ei_add_test(bicgstab)
251
ei_add_test(sparselu)
251
ei_add_test(sparselu)
252
ei_add_test(sparseqr)
252
ei_add_test(sparseqr)
253
253
254
ei_add_test(sparse_dense_outer_product)
255
254
# ei_add_test(denseLM)
256
# ei_add_test(denseLM)
255
257
256
if(QT4_FOUND)
258
if(QT4_FOUND)
257
  ei_add_test(qtvector "" "${QT_QTCORE_LIBRARY}")
259
  ei_add_test(qtvector "" "${QT_QTCORE_LIBRARY}")
258
endif(QT4_FOUND)
260
endif(QT4_FOUND)
259
261
260
if(UMFPACK_FOUND)
262
if(UMFPACK_FOUND)
261
  ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}")
263
  ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}")
(-)a/test/sparse_dense_outer_product.cpp (+87 lines)
Line 0 Link Here
1
#include "sparse.h"
2
3
typedef typename SparseMatrix<double>::Index Index;
4
5
void testit(Index rows, Index cols) 
6
{
7
  double density = (rows*cols==1) ? 1 : ((double)internal::random<int>(1,4)/4.0);
8
9
  std::cout << "Testing " << rows << "x" << cols << ", density = " << density << std::endl;
10
11
  Matrix<double, Dynamic, Dynamic> D1, D2;
12
  SparseMatrix<double> S1;
13
  SparseMatrix<double, RowMajorBit> S2;
14
15
  {
16
    Matrix<double, Dynamic, 1> A = Matrix<double, Dynamic, 1>::Random(rows);
17
    Matrix<double, Dynamic, Dynamic> R(1, cols);
18
    {
19
      SparseMatrix<double> B(1,cols);
20
21
      initSparse((cols==1) ? 1 : density, R, B);
22
      D2 = A * B.toDense();
23
      VERIFY_IS_APPROX(D1 = A * B, D2);
24
      VERIFY_IS_APPROX(S1 = A * B, D2);
25
      VERIFY_IS_APPROX(S2 = A * B, D2);
26
    }
27
    {
28
      SparseMatrix<double, RowMajorBit> B(1,cols);  
29
      initSparse((cols==1) ? 1 : density, R, B);
30
      D2 = A * B.toDense();
31
      VERIFY_IS_APPROX(D1 = A * B, D2);
32
      VERIFY_IS_APPROX(S1 = A * B, D2);
33
      VERIFY_IS_APPROX(S2 = A * B, D2);
34
    }
35
  }
36
  {
37
    Matrix<double, 1, Dynamic> A = Matrix<double, 1, Dynamic>::Random(cols);
38
    Matrix<double, Dynamic, Dynamic> R(rows, 1);
39
    {
40
      SparseMatrix<double> B(rows,1);
41
      initSparse((rows==1) ? 1 : density, R, B);
42
      D2 = B.toDense() * A;
43
      VERIFY_IS_APPROX(D1 = B * A, D2);
44
      VERIFY_IS_APPROX(S1 = B * A, D2);
45
      VERIFY_IS_APPROX(S2 = B * A, D2);
46
    }
47
    {
48
      SparseMatrix<double, RowMajorBit> B(rows, 1);
49
      initSparse((rows==1) ? 1 : density, R, B);
50
      D2 = B.toDense() * A;
51
      VERIFY_IS_APPROX(D1 = B * A, D2);
52
      VERIFY_IS_APPROX(S1 = B * A, D2);
53
      VERIFY_IS_APPROX(S2 = B * A, D2);
54
    }
55
  }
56
}
57
58
void test_sparse_dense_outer_product() {
59
#ifdef EIGEN_TEST_PART_1
60
  for(int i = 0; i < g_repeat; i++) testit(1, 1);
61
#endif
62
  
63
#ifdef EIGEN_TEST_PART_2
64
  Index n = 100;
65
  for(int i = 0; i < g_repeat; i++) {
66
      Index cols = internal::random<Index>(2,n);
67
      testit(1, cols);
68
  }
69
#endif
70
71
#ifdef EIGEN_TEST_PART_3
72
  Index n = 100;
73
  for(int i = 0; i < g_repeat; i++) {
74
    Index rows = internal::random<Index>(2,n);
75
    testit(rows, 1);
76
  }
77
#endif
78
79
#ifdef EIGEN_TEST_PART_4
80
  Index n = 100;
81
  for(int i = 0; i < g_repeat; i++) {
82
    Index rows  = internal::random<Index>(2,n);
83
    Index cols  = internal::random<Index>(2,n);
84
    testit(rows, cols);
85
  }
86
#endif
87
}

Return to bug 838