10 #ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H
11 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
29 template<
typename Lhs,
typename Rhs,
int UpLo>
30 class SparseSelfAdjointTimeDenseProduct;
32 template<
typename Lhs,
typename Rhs,
int UpLo>
33 class DenseTimeSparseSelfAdjointProduct;
37 template<
typename MatrixType,
unsigned int UpLo>
38 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
41 template<
int SrcUpLo,
int DstUpLo,
typename MatrixType,
int DestOrder>
42 void permute_symm_to_symm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm = 0);
44 template<
int UpLo,
typename MatrixType,
int DestOrder>
45 void permute_symm_to_fullsymm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm = 0);
50 :
public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
54 typedef typename MatrixType::Scalar Scalar;
55 typedef typename MatrixType::Index Index;
57 typedef typename MatrixType::Nested MatrixTypeNested;
58 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
62 eigen_assert(rows()==cols() &&
"SelfAdjointView is only for squared matrices");
65 inline Index rows()
const {
return m_matrix.rows(); }
66 inline Index cols()
const {
return m_matrix.cols(); }
69 const _MatrixTypeNested& matrix()
const {
return m_matrix; }
70 _MatrixTypeNested& matrix() {
return m_matrix.const_cast_derived(); }
73 template<
typename OtherDerived>
74 SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
77 return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
81 template<
typename OtherDerived>
friend
82 DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
85 return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
96 template<
typename DerivedU>
102 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
105 template<
typename DestScalar>
void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest)
const
108 SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols());
109 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
116 return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);
119 template<
typename SrcMatrixType,
int SrcUpLo>
120 SparseSelfAdjointView& operator=(
const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix)
122 permutedMatrix.evalTo(*
this);
127 SparseSelfAdjointView& operator=(
const SparseSelfAdjointView& src)
129 PermutationMatrix<Dynamic> pnull;
130 return *
this = src.twistedBy(pnull);
133 template<
typename SrcMatrixType,
unsigned int SrcUpLo>
134 SparseSelfAdjointView& operator=(
const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src)
136 PermutationMatrix<Dynamic> pnull;
137 return *
this = src.twistedBy(pnull);
146 typename MatrixType::Nested m_matrix;
147 mutable VectorI m_countPerRow;
148 mutable VectorI m_countPerCol;
155 template<
typename Derived>
156 template<
unsigned int UpLo>
157 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
const
162 template<
typename Derived>
163 template<
unsigned int UpLo>
164 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
173 template<
typename MatrixType,
unsigned int UpLo>
174 template<
typename DerivedU>
175 SparseSelfAdjointView<MatrixType,UpLo>&
178 SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
180 m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
182 m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
192 template<
typename Lhs,
typename Rhs,
int UpLo>
193 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
194 : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
196 typedef Dense StorageKind;
200 template<
typename Lhs,
typename Rhs,
int UpLo>
201 class SparseSelfAdjointTimeDenseProduct
202 :
public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
205 EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
207 SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
210 template<
typename Dest>
void scaleAndAddTo(Dest& dest, Scalar alpha)
const
212 EIGEN_ONLY_USED_FOR_DEBUG(alpha);
214 eigen_assert(alpha==Scalar(1) &&
"alpha != 1 is not implemented yet, sorry");
215 typedef typename internal::remove_all<Lhs>::type _Lhs;
216 typedef typename internal::remove_all<Rhs>::type _Rhs;
217 typedef typename _Lhs::InnerIterator LhsInnerIterator;
222 || ( (UpLo&
Upper) && !LhsIsRowMajor)
223 || ( (UpLo&
Lower) && LhsIsRowMajor),
224 ProcessSecondHalf = !ProcessFirstHalf
226 for (Index j=0; j<m_lhs.outerSize(); ++j)
228 LhsInnerIterator i(m_lhs,j);
229 if (ProcessSecondHalf)
231 while (i && i.index()<j) ++i;
232 if(i && i.index()==j)
234 dest.row(j) += i.value() * m_rhs.row(j);
238 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
240 Index a = LhsIsRowMajor ? j : i.index();
241 Index b = LhsIsRowMajor ? i.index() : j;
242 typename Lhs::Scalar v = i.value();
243 dest.row(a) += (v) * m_rhs.row(b);
244 dest.row(b) += internal::conj(v) * m_rhs.row(a);
246 if (ProcessFirstHalf && i && (i.index()==j))
247 dest.row(j) += i.value() * m_rhs.row(j);
252 SparseSelfAdjointTimeDenseProduct&
operator=(
const SparseSelfAdjointTimeDenseProduct&);
256 template<
typename Lhs,
typename Rhs,
int UpLo>
257 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
258 : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
262 template<
typename Lhs,
typename Rhs,
int UpLo>
263 class DenseTimeSparseSelfAdjointProduct
264 :
public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
267 EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
269 DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
272 template<
typename Dest>
void scaleAndAddTo(Dest& , Scalar )
const
278 DenseTimeSparseSelfAdjointProduct&
operator=(
const DenseTimeSparseSelfAdjointProduct&);
286 template<
typename MatrixType,
int UpLo>
287 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
290 template<
int UpLo,
typename MatrixType,
int DestOrder>
291 void permute_symm_to_fullsymm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm)
293 typedef typename MatrixType::Index Index;
294 typedef typename MatrixType::Scalar Scalar;
295 typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
296 typedef Matrix<Index,Dynamic,1> VectorI;
298 Dest& dest(_dest.derived());
300 StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
303 Index size = mat.rows();
307 dest.resize(size,size);
308 for(Index j = 0; j<size; ++j)
310 Index jp = perm ? perm[j] : j;
311 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
313 Index i = it.index();
316 Index ip = perm ? perm[i] : i;
318 count[StorageOrderMatch ? jp : ip]++;
321 else if(( UpLo==
Lower && r>c) || ( UpLo==
Upper && r<c))
328 Index nnz = count.sum();
331 dest.resizeNonZeros(nnz);
332 dest.outerIndexPtr()[0] = 0;
333 for(Index j=0; j<size; ++j)
334 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
335 for(Index j=0; j<size; ++j)
336 count[j] = dest.outerIndexPtr()[j];
339 for(Index j = 0; j<size; ++j)
341 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
343 Index i = it.index();
347 Index jp = perm ? perm[j] : j;
348 Index ip = perm ? perm[i] : i;
352 Index k = count[StorageOrderMatch ? jp : ip]++;
353 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
354 dest.valuePtr()[k] = it.value();
358 Index k = count[ip]++;
359 dest.innerIndexPtr()[k] = ip;
360 dest.valuePtr()[k] = it.value();
362 else if(( (UpLo&
Lower)==Lower && r>c) || ( (UpLo&
Upper)==Upper && r<c))
364 if(!StorageOrderMatch)
366 Index k = count[jp]++;
367 dest.innerIndexPtr()[k] = ip;
368 dest.valuePtr()[k] = it.value();
370 dest.innerIndexPtr()[k] = jp;
371 dest.valuePtr()[k] = internal::conj(it.value());
377 template<
int _SrcUpLo,
int _DstUpLo,
typename MatrixType,
int DstOrder>
378 void permute_symm_to_symm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm)
380 typedef typename MatrixType::Index Index;
381 typedef typename MatrixType::Scalar Scalar;
382 SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived());
383 typedef Matrix<Index,Dynamic,1> VectorI;
386 StorageOrderMatch = int(SrcOrder) == int(DstOrder),
387 DstUpLo = DstOrder==
RowMajor ? (_DstUpLo==Upper ? Lower :
Upper) : _DstUpLo,
388 SrcUpLo = SrcOrder==
RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
391 Index size = mat.rows();
394 dest.resize(size,size);
395 for(Index j = 0; j<size; ++j)
397 Index jp = perm ? perm[j] : j;
398 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
400 Index i = it.index();
401 if((
int(SrcUpLo)==
int(Lower) && i<j) || (
int(SrcUpLo)==
int(Upper) && i>j))
404 Index ip = perm ? perm[i] : i;
405 count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
408 dest.outerIndexPtr()[0] = 0;
409 for(Index j=0; j<size; ++j)
410 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
411 dest.resizeNonZeros(dest.outerIndexPtr()[size]);
412 for(Index j=0; j<size; ++j)
413 count[j] = dest.outerIndexPtr()[j];
415 for(Index j = 0; j<size; ++j)
418 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
420 Index i = it.index();
421 if((
int(SrcUpLo)==
int(Lower) && i<j) || (
int(SrcUpLo)==
int(Upper) && i>j))
424 Index jp = perm ? perm[j] : j;
425 Index ip = perm? perm[i] : i;
427 Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
428 dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
430 if(!StorageOrderMatch) std::swap(ip,jp);
431 if( ((
int(DstUpLo)==
int(Lower) && ip<jp) || (
int(DstUpLo)==
int(Upper) && ip>jp)))
432 dest.valuePtr()[k] = conj(it.value());
434 dest.valuePtr()[k] = it.value();
441 template<
typename MatrixType,
int UpLo>
442 class SparseSymmetricPermutationProduct
443 :
public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
446 typedef typename MatrixType::Scalar Scalar;
447 typedef typename MatrixType::Index Index;
449 typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm;
451 typedef Matrix<Index,Dynamic,1> VectorI;
452 typedef typename MatrixType::Nested MatrixTypeNested;
453 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
455 SparseSymmetricPermutationProduct(
const MatrixType& mat,
const Perm& perm)
456 : m_matrix(mat), m_perm(perm)
459 inline Index rows()
const {
return m_matrix.rows(); }
460 inline Index cols()
const {
return m_matrix.cols(); }
462 template<
typename DestScalar,
int Options,
typename DstIndex>
463 void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest)
const
465 internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
468 template<
typename DestType,
unsigned int DestUpLo>
void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest)
const
470 internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
474 MatrixTypeNested m_matrix;
481 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H