src/SparseCore/SparseAssign.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SPARSEASSIGN_H
0011 #define EIGEN_SPARSEASSIGN_H
0012
0013 namespace Eigen {
0014
0015 template<typename Derived>
0016 template<typename OtherDerived>
0017 Derived& SparseMatrixBase<Derived>::operator=(const EigenBase<OtherDerived> &other)
0018 {
0019   internal::call_assignment_no_alias(derived(), other.derived());
0020   return derived();
0021 }
0022
0023 template<typename Derived>
0024 template<typename OtherDerived>
0025 Derived& SparseMatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other)
0026 {
0027   // TODO use the evaluator mechanism
0028   other.evalTo(derived());
0029   return derived();
0030 }
0031
0032 template<typename Derived>
0033 template<typename OtherDerived>
0034 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseMatrixBase<OtherDerived>& other)
0035 {
0036   // by default sparse evaluation do not alias, so we can safely bypass the generic call_assignment routine
0037   internal::Assignment<Derived,OtherDerived,internal::assign_op<Scalar,typename OtherDerived::Scalar> >
0038           ::run(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
0039   return derived();
0040 }
0041
0042 template<typename Derived>
0043 inline Derived& SparseMatrixBase<Derived>::operator=(const Derived& other)
0044 {
0045   internal::call_assignment_no_alias(derived(), other.derived());
0046   return derived();
0047 }
0048
0049 namespace internal {
0050
0051 template<>
0052 struct storage_kind_to_evaluator_kind<Sparse> {
0053   typedef IteratorBased Kind;
0054 };
0055
0056 template<>
0057 struct storage_kind_to_shape<Sparse> {
0058   typedef SparseShape Shape;
0059 };
0060
0061 struct Sparse2Sparse {};
0062 struct Sparse2Dense  {};
0063
0064 template<> struct AssignmentKind<SparseShape, SparseShape>           { typedef Sparse2Sparse Kind; };
0065 template<> struct AssignmentKind<SparseShape, SparseTriangularShape> { typedef Sparse2Sparse Kind; };
0066 template<> struct AssignmentKind<DenseShape,  SparseShape>           { typedef Sparse2Dense  Kind; };
0067 template<> struct AssignmentKind<DenseShape,  SparseTriangularShape> { typedef Sparse2Dense  Kind; };
0068
0069
0070 template<typename DstXprType, typename SrcXprType>
0071 void assign_sparse_to_sparse(DstXprType &dst, const SrcXprType &src)
0072 {
0073   typedef typename DstXprType::Scalar Scalar;
0074   typedef internal::evaluator<DstXprType> DstEvaluatorType;
0075   typedef internal::evaluator<SrcXprType> SrcEvaluatorType;
0076
0077   SrcEvaluatorType srcEvaluator(src);
0078
0079   const bool transpose = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit);
0080   const Index outerEvaluationSize = (SrcEvaluatorType::Flags&RowMajorBit) ? src.rows() : src.cols();
0081   if ((!transpose) && src.isRValue())
0082   {
0083     // eval without temporary
0084     dst.resize(src.rows(), src.cols());
0085     dst.setZero();
0086     dst.reserve((std::min)(src.rows()*src.cols(), (std::max)(src.rows(),src.cols())*2));
0087     for (Index j=0; j<outerEvaluationSize; ++j)
0088     {
0089       dst.startVec(j);
0090       for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
0091       {
0092         Scalar v = it.value();
0093         dst.insertBackByOuterInner(j,it.index()) = v;
0094       }
0095     }
0096     dst.finalize();
0097   }
0098   else
0099   {
0100     // eval through a temporary
0101     eigen_assert(( ((internal::traits<DstXprType>::SupportedAccessPatterns & OuterRandomAccessPattern)==OuterRandomAccessPattern) ||
0102               (!((DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit)))) &&
0103               "the transpose operation is supposed to be handled in SparseMatrix::operator=");
0104
0105     enum { Flip = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit) };
0106
0107
0108     DstXprType temp(src.rows(), src.cols());
0109
0110     temp.reserve((std::min)(src.rows()*src.cols(), (std::max)(src.rows(),src.cols())*2));
0111     for (Index j=0; j<outerEvaluationSize; ++j)
0112     {
0113       temp.startVec(j);
0114       for (typename SrcEvaluatorType::InnerIterator it(srcEvaluator, j); it; ++it)
0115       {
0116         Scalar v = it.value();
0117         temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v;
0118       }
0119     }
0120     temp.finalize();
0121
0122     dst = temp.markAsRValue();
0123   }
0124 }
0125
0126 // Generic Sparse to Sparse assignment
0127 template< typename DstXprType, typename SrcXprType, typename Functor>
0128 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Sparse>
0129 {
0130   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0131   {
0132     assign_sparse_to_sparse(dst.derived(), src.derived());
0133   }
0134 };
0135
0136 // Generic Sparse to Dense assignment
0137 template< typename DstXprType, typename SrcXprType, typename Functor, typename Weak>
0138 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Weak>
0139 {
0140   static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
0141   {
0142     if(internal::is_same<Functor,internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> >::value)
0143       dst.setZero();
0144
0145     internal::evaluator<SrcXprType> srcEval(src);
0146     resize_if_allowed(dst, src, func);
0147     internal::evaluator<DstXprType> dstEval(dst);
0148
0149     const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols();
0150     for (Index j=0; j<outerEvaluationSize; ++j)
0151       for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i)
0152         func.assignCoeff(dstEval.coeffRef(i.row(),i.col()), i.value());
0153   }
0154 };
0155
0156 // Specialization for dense ?= dense +/- sparse and dense ?= sparse +/- dense
0157 template<typename DstXprType, typename Func1, typename Func2>
0158 struct assignment_from_dense_op_sparse
0159 {
0160   template<typename SrcXprType, typename InitialFunc>
0161   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0162   void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
0163   {
0164     #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN
0165     EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN
0166     #endif
0167
0168     call_assignment_no_alias(dst, src.lhs(), Func1());
0169     call_assignment_no_alias(dst, src.rhs(), Func2());
0170   }
0171
0172   // Specialization for dense1 = sparse + dense2; -> dense1 = dense2; dense1 += sparse;
0173   template<typename Lhs, typename Rhs, typename Scalar>
0174   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0175   typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type
0176   run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_sum_op<Scalar,Scalar>, const Lhs, const Rhs> &src,
0177       const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/)
0178   {
0179     #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN
0180     EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN
0181     #endif
0182
0183     // Apply the dense matrix first, then the sparse one.
0184     call_assignment_no_alias(dst, src.rhs(), Func1());
0185     call_assignment_no_alias(dst, src.lhs(), Func2());
0186   }
0187
0188   // Specialization for dense1 = sparse - dense2; -> dense1 = -dense2; dense1 += sparse;
0189   template<typename Lhs, typename Rhs, typename Scalar>
0190   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0191   typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type
0192   run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_difference_op<Scalar,Scalar>, const Lhs, const Rhs> &src,
0193       const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/)
0194   {
0195     #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN
0196     EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN
0197     #endif
0198
0199     // Apply the dense matrix first, then the sparse one.
0200     call_assignment_no_alias(dst, -src.rhs(), Func1());
0201     call_assignment_no_alias(dst,  src.lhs(), add_assign_op<typename DstXprType::Scalar,typename Lhs::Scalar>());
0202   }
0203 };
0204
0205 #define EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(ASSIGN_OP,BINOP,ASSIGN_OP2) \
0206   template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> \
0207   struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<Scalar,Scalar>, const Lhs, const Rhs>, internal::ASSIGN_OP<typename DstXprType::Scalar,Scalar>, \
0208                     Sparse2Dense, \
0209                     typename internal::enable_if<   internal::is_same<typename internal::evaluator_traits<Lhs>::Shape,DenseShape>::value \
0210                                                  || internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type> \
0211     : assignment_from_dense_op_sparse<DstXprType, internal::ASSIGN_OP<typename DstXprType::Scalar,typename Lhs::Scalar>, internal::ASSIGN_OP2<typename DstXprType::Scalar,typename Rhs::Scalar> > \
0212   {}
0213
0214 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op,    scalar_sum_op,add_assign_op);
0215 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_sum_op,add_assign_op);
0216 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_sum_op,sub_assign_op);
0217
0218 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op,    scalar_difference_op,sub_assign_op);
0219 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_difference_op,sub_assign_op);
0220 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_difference_op,add_assign_op);
0221
0222
0223 // Specialization for "dst = dec.solve(rhs)"
0224 // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error
0225 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
0226 struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar,Scalar>, Sparse2Sparse>
0227 {
0228   typedef Solve<DecType,RhsType> SrcXprType;
0229   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0230   {
0231     Index dstRows = src.rows();
0232     Index dstCols = src.cols();
0233     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0234       dst.resize(dstRows, dstCols);
0235
0236     src.dec()._solve_impl(src.rhs(), dst);
0237   }
0238 };
0239
0240 struct Diagonal2Sparse {};
0241
0242 template<> struct AssignmentKind<SparseShape,DiagonalShape> { typedef Diagonal2Sparse Kind; };
0243
0244 template< typename DstXprType, typename SrcXprType, typename Functor>
0245 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Sparse>
0246 {
0247   typedef typename DstXprType::StorageIndex StorageIndex;
0248   typedef typename DstXprType::Scalar Scalar;
0249
0250   template<int Options, typename AssignFunc>
0251   static void run(SparseMatrix<Scalar,Options,StorageIndex> &dst, const SrcXprType &src, const AssignFunc &func)
0252   { dst.assignDiagonal(src.diagonal(), func); }
0253
0254   template<typename DstDerived>
0255   static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0256   { dst.derived().diagonal() = src.diagonal(); }
0257
0258   template<typename DstDerived>
0259   static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0260   { dst.derived().diagonal() += src.diagonal(); }
0261
0262   template<typename DstDerived>
0263   static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0264   { dst.derived().diagonal() -= src.diagonal(); }
0265 };
0266 } // end namespace internal
0267
0268 } // end namespace Eigen
0269
0270 #endif // EIGEN_SPARSEASSIGN_H