Eigen  3.2.91
SparseAssign.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SPARSEASSIGN_H
11 #define EIGEN_SPARSEASSIGN_H
12 
13 namespace Eigen {
14 
15 template<typename Derived>
16 template<typename OtherDerived>
17 Derived& SparseMatrixBase<Derived>::operator=(const EigenBase<OtherDerived> &other)
18 {
19  internal::call_assignment_no_alias(derived(), other.derived());
20  return derived();
21 }
22 
23 template<typename Derived>
24 template<typename OtherDerived>
25 Derived& SparseMatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other)
26 {
27  // TODO use the evaluator mechanism
28  other.evalTo(derived());
29  return derived();
30 }
31 
32 template<typename Derived>
33 template<typename OtherDerived>
34 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseMatrixBase<OtherDerived>& other)
35 {
36  // by default sparse evaluation do not alias, so we can safely bypass the generic call_assignment routine
37  internal::Assignment<Derived,OtherDerived,internal::assign_op<Scalar> >
38  ::run(derived(), other.derived(), internal::assign_op<Scalar>());
39  return derived();
40 }
41 
42 template<typename Derived>
43 inline Derived& SparseMatrixBase<Derived>::operator=(const Derived& other)
44 {
45  internal::call_assignment_no_alias(derived(), other.derived());
46  return derived();
47 }
48 
49 namespace internal {
50 
51 template<>
52 struct storage_kind_to_evaluator_kind<Sparse> {
53  typedef IteratorBased Kind;
54 };
55 
56 template<>
57 struct storage_kind_to_shape<Sparse> {
58  typedef SparseShape Shape;
59 };
60 
61 struct Sparse2Sparse {};
62 struct Sparse2Dense {};
63 
64 template<> struct AssignmentKind<SparseShape, SparseShape> { typedef Sparse2Sparse Kind; };
65 template<> struct AssignmentKind<SparseShape, SparseTriangularShape> { typedef Sparse2Sparse Kind; };
66 template<> struct AssignmentKind<DenseShape, SparseShape> { typedef Sparse2Dense Kind; };
67 
68 
69 template<typename DstXprType, typename SrcXprType>
70 void assign_sparse_to_sparse(DstXprType &dst, const SrcXprType &src)
71 {
72  typedef typename DstXprType::Scalar Scalar;
73  typedef internal::evaluator<DstXprType> DstEvaluatorType;
74  typedef internal::evaluator<SrcXprType> SrcEvaluatorType;
75 
76  SrcEvaluatorType srcEvaluator(src);
77 
78  const bool transpose = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit);
79  const Index outerEvaluationSize = (SrcEvaluatorType::Flags&RowMajorBit) ? src.rows() : src.cols();
80  if ((!transpose) && src.isRValue())
81  {
82  // eval without temporary
83  dst.resize(src.rows(), src.cols());
84  dst.setZero();
85  dst.reserve((std::max)(src.rows(),src.cols())*2);
86  for (Index j=0; j<outerEvaluationSize; ++j)
87  {
88  dst.startVec(j);
89  for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
90  {
91  Scalar v = it.value();
92  dst.insertBackByOuterInner(j,it.index()) = v;
93  }
94  }
95  dst.finalize();
96  }
97  else
98  {
99  // eval through a temporary
100  eigen_assert(( ((internal::traits<DstXprType>::SupportedAccessPatterns & OuterRandomAccessPattern)==OuterRandomAccessPattern) ||
101  (!((DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit)))) &&
102  "the transpose operation is supposed to be handled in SparseMatrix::operator=");
103 
104  enum { Flip = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit) };
105 
106 
107  DstXprType temp(src.rows(), src.cols());
108 
109  temp.reserve((std::max)(src.rows(),src.cols())*2);
110  for (Index j=0; j<outerEvaluationSize; ++j)
111  {
112  temp.startVec(j);
113  for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
114  {
115  Scalar v = it.value();
116  temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v;
117  }
118  }
119  temp.finalize();
120 
121  dst = temp.markAsRValue();
122  }
123 }
124 
125 // Generic Sparse to Sparse assignment
126 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
127 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Sparse, Scalar>
128 {
129  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &/*func*/)
130  {
131  assign_sparse_to_sparse(dst.derived(), src.derived());
132  }
133 };
134 
135 // Sparse to Dense assignment
136 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
137 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Scalar>
138 {
139  static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
140  {
141  eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
142 
143  internal::evaluator<SrcXprType> srcEval(src);
144  internal::evaluator<DstXprType> dstEval(dst);
145  const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols();
146  for (Index j=0; j<outerEvaluationSize; ++j)
147  for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i)
148  func.assignCoeff(dstEval.coeffRef(i.row(),i.col()), i.value());
149  }
150 };
151 
152 template< typename DstXprType, typename SrcXprType, typename Scalar>
153 struct Assignment<DstXprType, SrcXprType, internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense, Scalar>
154 {
155  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &)
156  {
157  eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
158 
159  dst.setZero();
160  internal::evaluator<SrcXprType> srcEval(src);
161  internal::evaluator<DstXprType> dstEval(dst);
162  const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols();
163  for (Index j=0; j<outerEvaluationSize; ++j)
164  for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i)
165  dstEval.coeffRef(i.row(),i.col()) = i.value();
166  }
167 };
168 
169 // Specialization for "dst = dec.solve(rhs)"
170 // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error
171 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
172 struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar>, Sparse2Sparse, Scalar>
173 {
174  typedef Solve<DecType,RhsType> SrcXprType;
175  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
176  {
177  src.dec()._solve_impl(src.rhs(), dst);
178  }
179 };
180 
181 struct Diagonal2Sparse {};
182 
183 template<> struct AssignmentKind<SparseShape,DiagonalShape> { typedef Diagonal2Sparse Kind; };
184 
185 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
186 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Sparse, Scalar>
187 {
188  typedef typename DstXprType::StorageIndex StorageIndex;
189  typedef Array<StorageIndex,Dynamic,1> ArrayXI;
190  typedef Array<Scalar,Dynamic,1> ArrayXS;
191  template<int Options>
192  static void run(SparseMatrix<Scalar,Options,StorageIndex> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &/*func*/)
193  {
194  Index size = src.diagonal().size();
195  dst.makeCompressed();
196  dst.resizeNonZeros(size);
197  Map<ArrayXI>(dst.innerIndexPtr(), size).setLinSpaced(0,StorageIndex(size)-1);
198  Map<ArrayXI>(dst.outerIndexPtr(), size+1).setLinSpaced(0,StorageIndex(size));
199  Map<ArrayXS>(dst.valuePtr(), size) = src.diagonal();
200  }
201 
202  template<typename DstDerived>
203  static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &/*func*/)
204  {
205  dst.diagonal() = src.diagonal();
206  }
207 
208  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar> &/*func*/)
209  { dst.diagonal() += src.diagonal(); }
210 
211  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar> &/*func*/)
212  { dst.diagonal() -= src.diagonal(); }
213 };
214 } // end namespace internal
215 
216 } // end namespace Eigen
217 
218 #endif // EIGEN_SPARSEASSIGN_H
Definition: LDLT.h:16
const unsigned int RowMajorBit
Definition: Constants.h:53
Definition: Eigen_Colamd.h:54