1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SPARSE_SELFADJOINTVIEW_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SparseCore_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SparseSelfAdjointView 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the dense matrix storing the coefficients 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param UpLo can be either \c #Lower or \c #Upper 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is an expression of a sefladjoint matrix from a triangular part of a matrix 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and most of the time this is the only way that it is used. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa SparseMatrixBase::selfadjointView() 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSelfAdjointTimeDenseProduct; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass DenseTimeSparseSelfAdjointProduct; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> { 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder> 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0); 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int UpLo,typename MatrixType,int DestOrder> 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> > 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Nested MatrixTypeNested; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal \returns a reference to the nested matrix */ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const _MatrixTypeNested& matrix() const { return m_matrix; } 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); } 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix product between a sparse self-adjoint matrix \c *this and a sparse matrix \a rhs. 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename OtherDerived> 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived> 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez operator*(const SparseMatrixBase<OtherDerived>& rhs) const 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived>(*this, rhs.derived()); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix product between a sparse matrix \a lhs and a sparse self-adjoint matrix \a rhs. 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename OtherDerived> friend 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject > 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez operator*(const SparseMatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject>(lhs.derived(), rhs); 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo> 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& rhs) const 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived()); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> friend 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo> 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to \c *this 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * call this function with u.adjoint(). 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DerivedU> 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest); 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO directly evaluate into _dest; 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols()); 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _dest = tmp; 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of P H P^-1 */ 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename SrcMatrixType,int SrcUpLo> 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix) 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath permutedMatrix.evalTo(*this); 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic> pnull; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename SrcMatrixType,unsigned int SrcUpLo> 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src) 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic> pnull; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // const SparseLLT<PlainObject, UpLo> llt() const; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // const SparseLDLT<PlainObject, UpLo> ldlt() const; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Nested m_matrix; 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable VectorI m_countPerRow; 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable VectorI m_countPerCol; 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseMatrixBase methods 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseSelfAdjointView methods 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename DerivedU> 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseSelfAdjointView<MatrixType,UpLo>& 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezSparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha) 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint(); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(alpha==Scalar(0)) 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>(); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>(); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of sparse self-adjoint time dense matrix 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> > 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Dense StorageKind; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSelfAdjointTimeDenseProduct 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_ONLY_USED_FOR_DEBUG(alpha); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO use alpha 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry"); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<Lhs>::type _Lhs; 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _Lhs::InnerIterator LhsInnerIterator; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ProcessFirstHalf = 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ((UpLo&(Upper|Lower))==(Upper|Lower)) 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || ( (UpLo&Upper) && !LhsIsRowMajor) 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || ( (UpLo&Lower) && LhsIsRowMajor), 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ProcessSecondHalf = !ProcessFirstHalf 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<m_lhs.outerSize(); ++j) 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LhsInnerIterator i(m_lhs,j); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ProcessSecondHalf) 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (i && i.index()<j) ++i; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(i && i.index()==j) 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(j) += i.value() * m_rhs.row(j); 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index a = LhsIsRowMajor ? j : i.index(); 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index b = LhsIsRowMajor ? i.index() : j; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::Scalar v = i.value(); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(a) += (v) * m_rhs.row(b); 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.row(b) += numext::conj(v) * m_rhs.row(a); 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ProcessFirstHalf && i && (i.index()==j)) 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(j) += i.value() * m_rhs.row(j); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> > 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{}; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass DenseTimeSparseSelfAdjointProduct 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, const Scalar& /*alpha*/) const 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of symmetric copies and permutations 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int UpLo> 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> { 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int UpLo,typename MatrixType,int DestOrder> 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,DestOrder,Index> Dest; 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Dest& dest(_dest.derived()); 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor) 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.resize(size); 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(UpLo==(Upper|Lower)) 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[StorageOrderMatch ? jp : ip]++; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c)) 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[jp]++; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nnz = count.sum(); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // reserve space 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(nnz); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // copy data 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(UpLo==(Upper|Lower)) 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[StorageOrderMatch ? jp : ip]++; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp; 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[ip]++; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c)) 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(ip,jp); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[jp]++; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k = count[ip]++; 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = jp; 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.valuePtr()[k] = numext::conj(it.value()); 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder> 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived()); 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(SrcOrder) == int(DstOrder), 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo, 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count(size); 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(dest.outerIndexPtr()[size]); 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm? perm[i] : i; 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) std::swap(ip,jp); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp))) 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.valuePtr()[k] = numext::conj(it.value()); 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType,int UpLo> 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSymmetricPermutationProduct 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> > 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm; 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Nested MatrixTypeNested; 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(mat), m_perm(perm) 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar, int Options, typename DstIndex> 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data()); 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseMatrix<DestScalar,(Options&RowMajor)==RowMajor ? ColMajor : RowMajor, DstIndex> tmp; 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::permute_symm_to_fullsymm<UpLo>(m_matrix,tmp,m_perm.indices().data()); 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez _dest = tmp; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data()); 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixTypeNested m_matrix; 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Perm& m_perm; 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H 508