SparseSelfAdjointView.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
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 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo> 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& rhs) const 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived()); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> friend 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo> 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to \c *this 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * call this function with u.adjoint(). 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DerivedU> 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO directly evaluate into _dest; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols()); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _dest = tmp; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of P H P^-1 */ 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename SrcMatrixType,int SrcUpLo> 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix) 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath permutedMatrix.evalTo(*this); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic> pnull; 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename SrcMatrixType,unsigned int SrcUpLo> 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src) 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic> pnull; 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // const SparseLLT<PlainObject, UpLo> llt() const; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // const SparseLDLT<PlainObject, UpLo> ldlt() const; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Nested m_matrix; 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable VectorI m_countPerRow; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable VectorI m_countPerCol; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseMatrixBase methods 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseSelfAdjointView methods 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename DerivedU> 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseSelfAdjointView<MatrixType,UpLo>& 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha) 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint(); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(alpha==Scalar(0)) 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>(); 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>(); 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of sparse self-adjoint time dense matrix 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> > 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Dense StorageKind; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSelfAdjointTimeDenseProduct 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO use alpha 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry"); 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<Lhs>::type _Lhs; 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<Rhs>::type _Rhs; 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _Lhs::InnerIterator LhsInnerIterator; 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ProcessFirstHalf = 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ((UpLo&(Upper|Lower))==(Upper|Lower)) 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || ( (UpLo&Upper) && !LhsIsRowMajor) 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || ( (UpLo&Lower) && LhsIsRowMajor), 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ProcessSecondHalf = !ProcessFirstHalf 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<m_lhs.outerSize(); ++j) 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LhsInnerIterator i(m_lhs,j); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ProcessSecondHalf) 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (i && i.index()<j) ++i; 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(i && i.index()==j) 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(j) += i.value() * m_rhs.row(j); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index a = LhsIsRowMajor ? j : i.index(); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index b = LhsIsRowMajor ? i.index() : j; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::Scalar v = i.value(); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(a) += (v) * m_rhs.row(b); 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(b) += internal::conj(v) * m_rhs.row(a); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ProcessFirstHalf && i && (i.index()==j)) 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.row(j) += i.value() * m_rhs.row(j); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> > 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{}; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int UpLo> 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass DenseTimeSparseSelfAdjointProduct 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of symmetric copies and permutations 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int UpLo> 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> { 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int UpLo,typename MatrixType,int DestOrder> 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,DestOrder,Index> Dest; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Dest& dest(_dest.derived()); 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor) 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.resize(size); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(UpLo==(Upper|Lower)) 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[StorageOrderMatch ? jp : ip]++; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c)) 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[jp]++; 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nnz = count.sum(); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // reserve space 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(nnz); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // copy data 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(UpLo==(Upper|Lower)) 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[StorageOrderMatch ? jp : ip]++; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[ip]++; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c)) 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(ip,jp); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[jp]++; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k = count[ip]++; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = jp; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = internal::conj(it.value()); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder> 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived()); 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(SrcOrder) == int(DstOrder), 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo, 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count(size); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(dest.outerIndexPtr()[size]); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm? perm[i] : i; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) std::swap(ip,jp); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp))) 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = conj(it.value()); 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType,int UpLo> 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSymmetricPermutationProduct 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> > 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm; 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Index,Dynamic,1> VectorI; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Nested MatrixTypeNested; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(mat), m_perm(perm) 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestScalar, int Options, typename DstIndex> 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data()); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data()); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixTypeNested m_matrix; 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Perm& m_perm; 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H 481