SparseAssign.h (11368B)
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<> 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 template<> struct AssignmentKind<DenseShape, SparseTriangularShape> { 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::min)(src.rows()*src.cols(), (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::min)(src.rows()*src.cols(), (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, typename Weak> 138 struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Weak> 139 { 140 static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) 141 { 142 if(internal::is_same<Functor,internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> >::value) 143 dst.setZero(); 144 145 internal::evaluator<SrcXprType> srcEval(src); 146 resize_if_allowed(dst, src, func); 147 internal::evaluator<DstXprType> dstEval(dst); 148 149 const Index outerEvaluationSize = (internal::evaluator<SrcXprType>::Flags&RowMajorBit) ? src.rows() : src.cols(); 150 for (Index j=0; j<outerEvaluationSize; ++j) 151 for (typename internal::evaluator<SrcXprType>::InnerIterator i(srcEval,j); i; ++i) 152 func.assignCoeff(dstEval.coeffRef(i.row(),i.col()), i.value()); 153 } 154 }; 155 156 // Specialization for dense ?= dense +/- sparse and dense ?= sparse +/- dense 157 template<typename DstXprType, typename Func1, typename Func2> 158 struct assignment_from_dense_op_sparse 159 { 160 template<typename SrcXprType, typename InitialFunc> 161 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 162 void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/) 163 { 164 #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN 165 EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN 166 #endif 167 168 call_assignment_no_alias(dst, src.lhs(), Func1()); 169 call_assignment_no_alias(dst, src.rhs(), Func2()); 170 } 171 172 // Specialization for dense1 = sparse + dense2; -> dense1 = dense2; dense1 += sparse; 173 template<typename Lhs, typename Rhs, typename Scalar> 174 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 175 typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type 176 run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_sum_op<Scalar,Scalar>, const Lhs, const Rhs> &src, 177 const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/) 178 { 179 #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN 180 EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN 181 #endif 182 183 // Apply the dense matrix first, then the sparse one. 184 call_assignment_no_alias(dst, src.rhs(), Func1()); 185 call_assignment_no_alias(dst, src.lhs(), Func2()); 186 } 187 188 // Specialization for dense1 = sparse - dense2; -> dense1 = -dense2; dense1 += sparse; 189 template<typename Lhs, typename Rhs, typename Scalar> 190 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 191 typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type 192 run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_difference_op<Scalar,Scalar>, const Lhs, const Rhs> &src, 193 const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/) 194 { 195 #ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN 196 EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN 197 #endif 198 199 // Apply the dense matrix first, then the sparse one. 200 call_assignment_no_alias(dst, -src.rhs(), Func1()); 201 call_assignment_no_alias(dst, src.lhs(), add_assign_op<typename DstXprType::Scalar,typename Lhs::Scalar>()); 202 } 203 }; 204 205 #define EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(ASSIGN_OP,BINOP,ASSIGN_OP2) \ 206 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> \ 207 struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<Scalar,Scalar>, const Lhs, const Rhs>, internal::ASSIGN_OP<typename DstXprType::Scalar,Scalar>, \ 208 Sparse2Dense, \ 209 typename internal::enable_if< internal::is_same<typename internal::evaluator_traits<Lhs>::Shape,DenseShape>::value \ 210 || internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type> \ 211 : assignment_from_dense_op_sparse<DstXprType, internal::ASSIGN_OP<typename DstXprType::Scalar,typename Lhs::Scalar>, internal::ASSIGN_OP2<typename DstXprType::Scalar,typename Rhs::Scalar> > \ 212 {} 213 214 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op, scalar_sum_op,add_assign_op); 215 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_sum_op,add_assign_op); 216 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_sum_op,sub_assign_op); 217 218 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op, scalar_difference_op,sub_assign_op); 219 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_difference_op,sub_assign_op); 220 EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_difference_op,add_assign_op); 221 222 223 // Specialization for "dst = dec.solve(rhs)" 224 // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error 225 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar> 226 struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar,Scalar>, Sparse2Sparse> 227 { 228 typedef Solve<DecType,RhsType> SrcXprType; 229 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &) 230 { 231 Index dstRows = src.rows(); 232 Index dstCols = src.cols(); 233 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) 234 dst.resize(dstRows, dstCols); 235 236 src.dec()._solve_impl(src.rhs(), dst); 237 } 238 }; 239 240 struct Diagonal2Sparse {}; 241 242 template<> struct AssignmentKind<SparseShape,DiagonalShape> { typedef Diagonal2Sparse Kind; }; 243 244 template< typename DstXprType, typename SrcXprType, typename Functor> 245 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Sparse> 246 { 247 typedef typename DstXprType::StorageIndex StorageIndex; 248 typedef typename DstXprType::Scalar Scalar; 249 250 template<int Options, typename AssignFunc> 251 static void run(SparseMatrix<Scalar,Options,StorageIndex> &dst, const SrcXprType &src, const AssignFunc &func) 252 { dst.assignDiagonal(src.diagonal(), func); } 253 254 template<typename DstDerived> 255 static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 256 { dst.derived().diagonal() = src.diagonal(); } 257 258 template<typename DstDerived> 259 static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 260 { dst.derived().diagonal() += src.diagonal(); } 261 262 template<typename DstDerived> 263 static void run(SparseMatrixBase<DstDerived> &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 264 { dst.derived().diagonal() -= src.diagonal(); } 265 }; 266 } // end namespace internal 267 268 } // end namespace Eigen 269 270 #endif // EIGEN_SPARSEASSIGN_H