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>&
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
213 eigen_assert(alpha==Scalar(1) &&
"alpha != 1 is not implemented yet, sorry");
214 typedef typename internal::remove_all<Lhs>::type _Lhs;
215 typedef typename internal::remove_all<Rhs>::type _Rhs;
216 typedef typename _Lhs::InnerIterator LhsInnerIterator;
221 || ( (UpLo&
Upper) && !LhsIsRowMajor)
222 || ( (UpLo&
Lower) && LhsIsRowMajor),
223 ProcessSecondHalf = !ProcessFirstHalf
225 for (Index j=0; j<m_lhs.outerSize(); ++j)
227 LhsInnerIterator i(m_lhs,j);
228 if (ProcessSecondHalf)
230 while (i && i.index()<j) ++i;
231 if(i && i.index()==j)
233 dest.row(j) += i.value() * m_rhs.row(j);
237 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
239 Index a = LhsIsRowMajor ? j : i.index();
240 Index b = LhsIsRowMajor ? i.index() : j;
241 typename Lhs::Scalar v = i.value();
242 dest.row(a) += (v) * m_rhs.row(b);
243 dest.row(b) += internal::conj(v) * m_rhs.row(a);
245 if (ProcessFirstHalf && i && (i.index()==j))
246 dest.row(j) += i.value() * m_rhs.row(j);
251 SparseSelfAdjointTimeDenseProduct&
operator=(
const SparseSelfAdjointTimeDenseProduct&);
255 template<
typename Lhs,
typename Rhs,
int UpLo>
256 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
257 : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
261 template<
typename Lhs,
typename Rhs,
int UpLo>
262 class DenseTimeSparseSelfAdjointProduct
263 :
public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
266 EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
268 DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
271 template<
typename Dest>
void scaleAndAddTo(Dest& , Scalar )
const
277 DenseTimeSparseSelfAdjointProduct&
operator=(
const DenseTimeSparseSelfAdjointProduct&);
285 template<
typename MatrixType,
int UpLo>
286 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
289 template<
int UpLo,
typename MatrixType,
int DestOrder>
290 void permute_symm_to_fullsymm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm)
292 typedef typename MatrixType::Index Index;
293 typedef typename MatrixType::Scalar Scalar;
294 typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
295 typedef Matrix<Index,Dynamic,1> VectorI;
297 Dest& dest(_dest.derived());
299 StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
302 Index size = mat.rows();
306 dest.resize(size,size);
307 for(Index j = 0; j<size; ++j)
309 Index jp = perm ? perm[j] : j;
310 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
312 Index i = it.index();
315 Index ip = perm ? perm[i] : i;
317 count[StorageOrderMatch ? jp : ip]++;
320 else if(( UpLo==
Lower && r>c) || ( UpLo==
Upper && r<c))
327 Index nnz = count.sum();
330 dest.resizeNonZeros(nnz);
331 dest.outerIndexPtr()[0] = 0;
332 for(Index j=0; j<size; ++j)
333 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
334 for(Index j=0; j<size; ++j)
335 count[j] = dest.outerIndexPtr()[j];
338 for(Index j = 0; j<size; ++j)
340 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
342 Index i = it.index();
346 Index jp = perm ? perm[j] : j;
347 Index ip = perm ? perm[i] : i;
351 Index k = count[StorageOrderMatch ? jp : ip]++;
352 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
353 dest.valuePtr()[k] = it.value();
357 Index k = count[ip]++;
358 dest.innerIndexPtr()[k] = ip;
359 dest.valuePtr()[k] = it.value();
361 else if(( (UpLo&
Lower)==Lower && r>c) || ( (UpLo&
Upper)==Upper && r<c))
363 if(!StorageOrderMatch)
365 Index k = count[jp]++;
366 dest.innerIndexPtr()[k] = ip;
367 dest.valuePtr()[k] = it.value();
369 dest.innerIndexPtr()[k] = jp;
370 dest.valuePtr()[k] = internal::conj(it.value());
376 template<
int _SrcUpLo,
int _DstUpLo,
typename MatrixType,
int DstOrder>
377 void permute_symm_to_symm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm)
379 typedef typename MatrixType::Index Index;
380 typedef typename MatrixType::Scalar Scalar;
381 SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived());
382 typedef Matrix<Index,Dynamic,1> VectorI;
385 StorageOrderMatch = int(SrcOrder) == int(DstOrder),
386 DstUpLo = DstOrder==
RowMajor ? (_DstUpLo==Upper ? Lower :
Upper) : _DstUpLo,
387 SrcUpLo = SrcOrder==
RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
390 Index size = mat.rows();
393 dest.resize(size,size);
394 for(Index j = 0; j<size; ++j)
396 Index jp = perm ? perm[j] : j;
397 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
399 Index i = it.index();
400 if((
int(SrcUpLo)==
int(Lower) && i<j) || (
int(SrcUpLo)==
int(Upper) && i>j))
403 Index ip = perm ? perm[i] : i;
404 count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
407 dest.outerIndexPtr()[0] = 0;
408 for(Index j=0; j<size; ++j)
409 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
410 dest.resizeNonZeros(dest.outerIndexPtr()[size]);
411 for(Index j=0; j<size; ++j)
412 count[j] = dest.outerIndexPtr()[j];
414 for(Index j = 0; j<size; ++j)
417 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
419 Index i = it.index();
420 if((
int(SrcUpLo)==
int(Lower) && i<j) || (
int(SrcUpLo)==
int(Upper) && i>j))
423 Index jp = perm ? perm[j] : j;
424 Index ip = perm? perm[i] : i;
426 Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
427 dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
429 if(!StorageOrderMatch) std::swap(ip,jp);
430 if( ((
int(DstUpLo)==
int(Lower) && ip<jp) || (
int(DstUpLo)==
int(Upper) && ip>jp)))
431 dest.valuePtr()[k] = conj(it.value());
433 dest.valuePtr()[k] = it.value();
440 template<
typename MatrixType,
int UpLo>
441 class SparseSymmetricPermutationProduct
442 :
public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
445 typedef typename MatrixType::Scalar Scalar;
446 typedef typename MatrixType::Index Index;
448 typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm;
450 typedef Matrix<Index,Dynamic,1> VectorI;
451 typedef typename MatrixType::Nested MatrixTypeNested;
452 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
454 SparseSymmetricPermutationProduct(
const MatrixType& mat,
const Perm& perm)
455 : m_matrix(mat), m_perm(perm)
458 inline Index rows()
const {
return m_matrix.rows(); }
459 inline Index cols()
const {
return m_matrix.cols(); }
461 template<
typename DestScalar,
int Options,
typename DstIndex>
462 void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest)
const
464 internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
467 template<
typename DestType,
unsigned int DestUpLo>
void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest)
const
469 internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
473 MatrixTypeNested m_matrix;
480 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H