compbio
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,typename OtherDerived::Scalar> >
38  ::run(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::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<>
53  typedef IteratorBased Kind;
54 };
55 
56 template<>
58  typedef SparseShape Shape;
59 };
60 
61 struct Sparse2Sparse {};
62 struct Sparse2Dense {};
63 
64 template<> struct AssignmentKind<SparseShape, SparseShape> { typedef Sparse2Sparse Kind; };
66 template<> struct AssignmentKind<DenseShape, SparseShape> { typedef Sparse2Dense Kind; };
68 
69 
70 template<typename DstXprType, typename SrcXprType>
71 void assign_sparse_to_sparse(DstXprType &dst, const SrcXprType &src)
72 {
73  typedef typename DstXprType::Scalar Scalar;
74  typedef internal::evaluator<DstXprType> DstEvaluatorType;
75  typedef internal::evaluator<SrcXprType> SrcEvaluatorType;
76 
77  SrcEvaluatorType srcEvaluator(src);
78 
79  const bool transpose = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit);
80  const Index outerEvaluationSize = (SrcEvaluatorType::Flags&RowMajorBit) ? src.rows() : src.cols();
81  if ((!transpose) && src.isRValue())
82  {
83  // eval without temporary
84  dst.resize(src.rows(), src.cols());
85  dst.setZero();
86  dst.reserve((std::max)(src.rows(),src.cols())*2);
87  for (Index j=0; j<outerEvaluationSize; ++j)
88  {
89  dst.startVec(j);
90  for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
91  {
92  Scalar v = it.value();
93  dst.insertBackByOuterInner(j,it.index()) = v;
94  }
95  }
96  dst.finalize();
97  }
98  else
99  {
100  // eval through a temporary
101  eigen_assert(( ((internal::traits<DstXprType>::SupportedAccessPatterns & OuterRandomAccessPattern)==OuterRandomAccessPattern) ||
102  (!((DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit)))) &&
103  "the transpose operation is supposed to be handled in SparseMatrix::operator=");
104 
105  enum { Flip = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit) };
106 
107 
108  DstXprType temp(src.rows(), src.cols());
109 
110  temp.reserve((std::max)(src.rows(),src.cols())*2);
111  for (Index j=0; j<outerEvaluationSize; ++j)
112  {
113  temp.startVec(j);
114  for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
115  {
116  Scalar v = it.value();
117  temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v;
118  }
119  }
120  temp.finalize();
121 
122  dst = temp.markAsRValue();
123  }
124 }
125 
126 // Generic Sparse to Sparse assignment
127 template< typename DstXprType, typename SrcXprType, typename Functor>
128 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Sparse>
129 {
130  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
131  {
132  assign_sparse_to_sparse(dst.derived(), src.derived());
133  }
134 };
135 
136 // Generic Sparse to Dense assignment
137 template< typename DstXprType, typename SrcXprType, typename Functor>
138 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense>
139 {
140  static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
141  {
143  dst.setZero();
144 
145  internal::evaluator<SrcXprType> srcEval(src);
146  Index dstRows = src.rows();
147  Index dstCols = src.cols();
148  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
149  dst.resize(dstRows, dstCols);
150  internal::evaluator<DstXprType> dstEval(dst);
151 
152  const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols();
153  for (Index j=0; j<outerEvaluationSize; ++j)
154  for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i)
155  func.assignCoeff(dstEval.coeffRef(i.row(),i.col()), i.value());
156  }
157 };
158 
159 // Specialization for "dst = dec.solve(rhs)"
160 // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error
161 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
162 struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar,Scalar>, Sparse2Sparse>
163 {
165  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
166  {
167  Index dstRows = src.rows();
168  Index dstCols = src.cols();
169  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
170  dst.resize(dstRows, dstCols);
171 
172  src.dec()._solve_impl(src.rhs(), dst);
173  }
174 };
175 
176 struct Diagonal2Sparse {};
177 
179 
180 template< typename DstXprType, typename SrcXprType, typename Functor>
181 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Sparse>
182 {
183  typedef typename DstXprType::StorageIndex StorageIndex;
184  typedef typename DstXprType::Scalar Scalar;
187  template<int Options>
189  {
190  Index dstRows = src.rows();
191  Index dstCols = src.cols();
192  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
193  dst.resize(dstRows, dstCols);
194 
195  Index size = src.diagonal().size();
196  dst.makeCompressed();
197  dst.resizeNonZeros(size);
198  Map<ArrayXI>(dst.innerIndexPtr(), size).setLinSpaced(0,StorageIndex(size)-1);
199  Map<ArrayXI>(dst.outerIndexPtr(), size+1).setLinSpaced(0,StorageIndex(size));
200  Map<ArrayXS>(dst.valuePtr(), size) = src.diagonal();
201  }
202 
203  template<typename DstDerived>
204  static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
205  {
206  dst.diagonal() = src.diagonal();
207  }
208 
209  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
210  { dst.diagonal() += src.diagonal(); }
211 
212  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
213  { dst.diagonal() -= src.diagonal(); }
214 };
215 } // end namespace internal
216 
217 } // end namespace Eigen
218 
219 #endif // EIGEN_SPARSEASSIGN_H
Definition: NonLinearOptimization.cpp:108
const Scalar * valuePtr() const
Definition: SparseMatrix.h:144
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
Definition: AssignmentFunctors.h:67
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:88
Definition: Meta.h:63
Definition: CoreEvaluators.h:90
Namespace containing all symbols from the Eigen library.
Definition: bench_norm.cpp:85
Definition: Constants.h:521
Index cols() const
Definition: SparseMatrix.h:134
const unsigned int RowMajorBit
for a matrix, this means that the storage order is row-major.
Definition: Constants.h:61
Definition: Constants.h:512
Definition: CoreEvaluators.h:29
Definition: Constants.h:529
Definition: AssignmentFunctors.h:21
Definition: AssignEvaluator.h:753
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:281
Definition: Constants.h:515
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Index rows() const
Definition: SparseMatrix.h:132
Definition: SparseAssign.h:62
void resize(Index rows, Index cols)
Resizes the matrix to a rows x cols matrix and initializes it to zero.
Definition: SparseMatrix.h:617
The type used to identify a general sparse storage.
Definition: Constants.h:494
void makeCompressed()
Turns the matrix into the compressed format.
Definition: SparseMatrix.h:460
Definition: SparseUtil.h:137
Definition: benchGeometry.cpp:23
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:153
Definition: BandTriangularSolver.h:13
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:162
Definition: AssignEvaluator.h:740
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:45
Definition: AssignmentFunctors.h:46
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Definition: ForwardDeclarations.h:17
Definition: SparseAssign.h:61
Definition: SparseAssign.h:176