1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2009-2014 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 { 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 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 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \param Mode 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 Kamathnamespace internal { 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, unsigned int Mode> 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct traits<SparseSelfAdjointView<MatrixType,Mode> > : traits<MatrixType> { 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<int SrcMode,int DstMode,typename MatrixType,int DestOrder> 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::StorageIndex>& _dest, const typename MatrixType::StorageIndex* perm = 0); 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<int Mode,typename MatrixType,int DestOrder> 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::StorageIndex>& _dest, const typename MatrixType::StorageIndex* perm = 0); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, unsigned int _Mode> class SparseSelfAdjointView 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public EigenBase<SparseSelfAdjointView<MatrixType,_Mode> > 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Mode = _Mode, 50eda03298de395cf6217486971e6529f92da8da79Miao Wang TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0), 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = internal::traits<SparseSelfAdjointView>::RowsAtCompileTime, 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = internal::traits<SparseSelfAdjointView>::ColsAtCompileTime 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef EigenBase<SparseSelfAdjointView> Base; 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex,Dynamic,1> VectorI; 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::ref_selector<MatrixType>::non_const_type MatrixTypeNested; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit inline SparseSelfAdjointView(MatrixType& matrix) : m_matrix(matrix) 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal \returns a reference to the nested matrix */ 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const _MatrixTypeNested& matrix() const { return m_matrix; } 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename internal::remove_reference<MatrixTypeNested>::type& matrix() { return m_matrix; } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix product between a sparse self-adjoint matrix \c *this and a sparse matrix \a rhs. 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename OtherDerived> 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Product<SparseSelfAdjointView, OtherDerived> 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez operator*(const SparseMatrixBase<OtherDerived>& rhs) const 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Product<SparseSelfAdjointView, OtherDerived>(*this, rhs.derived()); 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix product between a sparse matrix \a lhs and a sparse self-adjoint matrix \a rhs. 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename OtherDerived> friend 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Product<OtherDerived, SparseSelfAdjointView> 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez operator*(const SparseMatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Product<OtherDerived, SparseSelfAdjointView>(lhs.derived(), rhs); 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> 1002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Product<SparseSelfAdjointView,OtherDerived> 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& rhs) const 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Product<SparseSelfAdjointView,OtherDerived>(*this, rhs.derived()); 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived> friend 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Product<OtherDerived,SparseSelfAdjointView> 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Product<OtherDerived,SparseSelfAdjointView>(lhs.derived(), rhs); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to \c *this 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * call this function with u.adjoint(). 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DerivedU> 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of P H P^-1 */ 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // TODO implement twists in a more evaluator friendly fashion 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm) const 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode>(m_matrix, perm); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename SrcMatrixType,int SrcMode> 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcMode>& permutedMatrix) 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::call_assignment_no_alias_no_transpose(*this, permutedMatrix); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> pnull; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename SrcMatrixType,unsigned int SrcMode> 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcMode>& src) 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> pnull; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this = src.twistedBy(pnull); 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void resize(Index rows, Index cols) 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EIGEN_ONLY_USED_FOR_DEBUG(rows); 1552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EIGEN_ONLY_USED_FOR_DEBUG(cols); 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(rows == this->rows() && cols == this->cols() 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang && "SparseSelfadjointView::resize() does not actually allow to resize."); 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixTypeNested m_matrix; 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //mutable VectorI m_countPerRow; 1642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //mutable VectorI m_countPerCol; 1652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang private: 1662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Dest> void evalTo(Dest &) const; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseMatrixBase methods 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 1752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtypename SparseMatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type SparseMatrixBase<Derived>::selfadjointView() const 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return SparseSelfAdjointView<const Derived, UpLo>(derived()); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<unsigned int UpLo> 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtypename SparseMatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type SparseMatrixBase<Derived>::selfadjointView() 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return SparseSelfAdjointView<Derived, UpLo>(derived()); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of SparseSelfAdjointView methods 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, unsigned int Mode> 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename DerivedU> 1932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangSparseSelfAdjointView<MatrixType,Mode>& 1942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangSparseSelfAdjointView<MatrixType,Mode>::rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha) 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<Scalar,(MatrixType::Flags&RowMajorBit)?RowMajor:ColMajor> tmp = u * u.adjoint(); 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(alpha==Scalar(0)) 1982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix = tmp.template triangularView<Mode>(); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 2002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix += alpha * tmp.template triangularView<Mode>(); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 2062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> 2082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// in the future selfadjoint-ness should be defined by the expression traits 2092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work) 2102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, unsigned int Mode> 2112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct evaluator_traits<SparseSelfAdjointView<MatrixType,Mode> > 2122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; 2142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseSelfAdjointShape Shape; 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 2162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct SparseSelfAdjoint2Sparse {}; 2182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<> struct AssignmentKind<SparseShape,SparseSelfAdjointShape> { typedef SparseSelfAdjoint2Sparse Kind; }; 2202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<> struct AssignmentKind<SparseSelfAdjointShape,SparseShape> { typedef Sparse2Sparse Kind; }; 2212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate< typename DstXprType, typename SrcXprType, typename Functor> 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct Assignment<DstXprType, SrcXprType, Functor, SparseSelfAdjoint2Sparse> 2242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename DstXprType::StorageIndex StorageIndex; 2262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> AssignOpType; 2272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestScalar,int StorageOrder> 2292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseMatrix<DestScalar,StorageOrder,StorageIndex> &dst, const SrcXprType &src, const AssignOpType&/*func*/) 2302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::permute_symm_to_fullsymm<SrcXprType::Mode>(src.matrix(), dst); 2322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // FIXME: the handling of += and -= in sparse matrices should be cleanup so that next two overloads could be reduced to: 2352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestScalar,int StorageOrder,typename AssignFunc> 2362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseMatrix<DestScalar,StorageOrder,StorageIndex> &dst, const SrcXprType &src, const AssignFunc& func) 2372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<DestScalar,StorageOrder,StorageIndex> tmp(src.rows(),src.cols()); 2392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang run(tmp, src, AssignOpType()); 2402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang call_assignment_no_alias_no_transpose(dst, tmp, func); 2412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestScalar,int StorageOrder> 2442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseMatrix<DestScalar,StorageOrder,StorageIndex> &dst, const SrcXprType &src, 2452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>& /* func */) 2462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<DestScalar,StorageOrder,StorageIndex> tmp(src.rows(),src.cols()); 2482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang run(tmp, src, AssignOpType()); 2492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dst += tmp; 2502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestScalar,int StorageOrder> 2532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseMatrix<DestScalar,StorageOrder,StorageIndex> &dst, const SrcXprType &src, 2542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>& /* func */) 2552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<DestScalar,StorageOrder,StorageIndex> tmp(src.rows(),src.cols()); 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang run(tmp, src, AssignOpType()); 2582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dst -= tmp; 2592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestScalar> 2622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(DynamicSparseMatrix<DestScalar,ColMajor,StorageIndex>& dst, const SrcXprType &src, const AssignOpType&/*func*/) 2632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // TODO directly evaluate into dst; 2652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<DestScalar,ColMajor,StorageIndex> tmp(dst.rows(),dst.cols()); 2662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::permute_symm_to_fullsymm<SrcXprType::Mode>(src.matrix(), tmp); 2672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dst = tmp; 2682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 2702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 2722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of sparse self-adjoint time dense matrix 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<int Mode, typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType> 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wanginline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha) 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EIGEN_ONLY_USED_FOR_DEBUG(alpha); 2832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::nested_eval<SparseLhsType,DenseRhsType::MaxColsAtCompileTime>::type SparseLhsTypeNested; 2852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::remove_all<SparseLhsTypeNested>::type SparseLhsTypeNestedCleaned; 2862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef evaluator<SparseLhsTypeNestedCleaned> LhsEval; 2872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename LhsEval::InnerIterator LhsIterator; 2882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename SparseLhsType::Scalar LhsScalar; 2892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 2912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsIsRowMajor = (LhsEval::Flags&RowMajorBit)==RowMajorBit, 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ProcessFirstHalf = 2932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ((Mode&(Upper|Lower))==(Upper|Lower)) 2942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang || ( (Mode&Upper) && !LhsIsRowMajor) 2952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang || ( (Mode&Lower) && LhsIsRowMajor), 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ProcessSecondHalf = !ProcessFirstHalf 2972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 2982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseLhsTypeNested lhs_nested(lhs); 3002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsEval lhsEval(lhs_nested); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // work on one column at once 3032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (Index k=0; k<rhs.cols(); ++k) 3042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (Index j=0; j<lhs.outerSize(); ++j) 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsIterator i(lhsEval,j); 3082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // handle diagonal coeff 3092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (ProcessSecondHalf) 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while (i && i.index()<j) ++i; 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(i && i.index()==j) 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res(j,k) += alpha * i.value() * rhs(j,k); 3152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++i; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // premultiplied rhs for scatters 3202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename ScalarBinaryOpTraits<AlphaType, typename DenseRhsType::Scalar>::ReturnType rhs_j(alpha*rhs(j,k)); 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // accumulator for partial scalar product 3222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename DenseResType::Scalar res_j(0); 3232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) 3242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsScalar lhs_ij = i.value(); 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!LhsIsRowMajor) lhs_ij = numext::conj(lhs_ij); 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res_j += lhs_ij * rhs(i.index(),k); 3282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res(i.index(),k) += numext::conj(lhs_ij) * rhs_j; 3292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res(j,k) += alpha * res_j; 3312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // handle diagonal coeff 3332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (ProcessFirstHalf && i && (i.index()==j)) 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res(j,k) += alpha * i.value() * rhs(j,k); 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 3362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename LhsView, typename Rhs, int ProductType> 3412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct generic_product_impl<LhsView, Rhs, SparseSelfAdjointShape, DenseShape, ProductType> 3422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang: generic_product_impl_base<LhsView, Rhs, generic_product_impl<LhsView, Rhs, SparseSelfAdjointShape, DenseShape, ProductType> > 3432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 3442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Dest> 3452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void scaleAndAddTo(Dest& dst, const LhsView& lhsView, const Rhs& rhs, const typename Dest::Scalar& alpha) 3462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename LhsView::_MatrixTypeNested Lhs; 3482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename nested_eval<Lhs,Dynamic>::type LhsNested; 3492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename nested_eval<Rhs,Dynamic>::type RhsNested; 3502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsNested lhsNested(lhsView.matrix()); 3512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RhsNested rhsNested(rhs); 3522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::sparse_selfadjoint_time_dense_product<LhsView::Mode>(lhsNested, rhsNested, dst, alpha); 3542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename Lhs, typename RhsView, int ProductType> 3582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct generic_product_impl<Lhs, RhsView, DenseShape, SparseSelfAdjointShape, ProductType> 3592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang: generic_product_impl_base<Lhs, RhsView, generic_product_impl<Lhs, RhsView, DenseShape, SparseSelfAdjointShape, ProductType> > 3602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 3612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Dest> 3622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const RhsView& rhsView, const typename Dest::Scalar& alpha) 3632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename RhsView::_MatrixTypeNested Rhs; 3652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename nested_eval<Lhs,Dynamic>::type LhsNested; 3662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename nested_eval<Rhs,Dynamic>::type RhsNested; 3672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LhsNested lhsNested(lhs); 3682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RhsNested rhsNested(rhsView.matrix()); 3692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // transpose everything 3712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Transpose<Dest> dstT(dst); 372eda03298de395cf6217486971e6529f92da8da79Miao Wang internal::sparse_selfadjoint_time_dense_product<RhsView::TransposeMode>(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha); 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// NOTE: these two overloads are needed to evaluate the sparse selfadjoint view into a full sparse matrix 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// TODO: maybe the copy could be handled by generic_product_impl so that these overloads would not be needed anymore 3782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename LhsView, typename Rhs, int ProductTag> 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct product_evaluator<Product<LhsView, Rhs, DefaultProduct>, ProductTag, SparseSelfAdjointShape, SparseShape> 3812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public evaluator<typename Product<typename Rhs::PlainObject, Rhs, DefaultProduct>::PlainObject> 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Product<LhsView, Rhs, DefaultProduct> XprType; 3842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename XprType::PlainObject PlainObject; 3852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef evaluator<PlainObject> Base; 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang product_evaluator(const XprType& xpr) 3882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_lhs(xpr.lhs()), m_result(xpr.rows(), xpr.cols()) 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ::new (static_cast<Base*>(this)) Base(m_result); 3912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_product_impl<typename Rhs::PlainObject, Rhs, SparseShape, SparseShape, ProductTag>::evalTo(m_result, m_lhs, xpr.rhs()); 3922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprotected: 3952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename Rhs::PlainObject m_lhs; 3962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PlainObject m_result; 3972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename Lhs, typename RhsView, int ProductTag> 4002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct product_evaluator<Product<Lhs, RhsView, DefaultProduct>, ProductTag, SparseShape, SparseSelfAdjointShape> 4012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public evaluator<typename Product<Lhs, typename Lhs::PlainObject, DefaultProduct>::PlainObject> 4022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 4032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Product<Lhs, RhsView, DefaultProduct> XprType; 4042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename XprType::PlainObject PlainObject; 4052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef evaluator<PlainObject> Base; 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang product_evaluator(const XprType& xpr) 4082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_rhs(xpr.rhs()), m_result(xpr.rows(), xpr.cols()) 4092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ::new (static_cast<Base*>(this)) Base(m_result); 4112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_product_impl<Lhs, typename Lhs::PlainObject, SparseShape, SparseShape, ProductTag>::evalTo(m_result, xpr.lhs(), m_rhs); 4122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprotected: 4152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename Lhs::PlainObject m_rhs; 4162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PlainObject m_result; 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // namespace internal 4202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*************************************************************************** 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of symmetric copies and permutations 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath***************************************************************************/ 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<int Mode,typename MatrixType,int DestOrder> 4272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::StorageIndex>& _dest, const typename MatrixType::StorageIndex* perm) 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 4312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar,DestOrder,StorageIndex> Dest; 4322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex,Dynamic,1> VectorI; 4332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef evaluator<MatrixType> MatEval; 4342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename evaluator<MatrixType>::InnerIterator MatIterator; 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatEval matEval(mat); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Dest& dest(_dest.derived()); 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor) 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count; 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.resize(size); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j = 0; j<size; ++j) 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index jp = perm ? perm[j] : j; 4502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(MatIterator it(matEval,j); it; ++it) 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = it.index(); 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ip = perm ? perm[i] : i; 4562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(Mode==(Upper|Lower)) 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[StorageOrderMatch ? jp : ip]++; 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 4602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else if(( Mode==Lower && r>c) || ( Mode==Upper && r<c)) 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[ip]++; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[jp]++; 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nnz = count.sum(); 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // reserve space 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(nnz); 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // copy data 4782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex j = 0; j<size; ++j) 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(MatIterator it(matEval,j); it; ++it) 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex i = internal::convert_index<StorageIndex>(it.index()); 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index r = it.row(); 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index c = it.col(); 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex jp = perm ? perm[j] : j; 4872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex ip = perm ? perm[i] : i; 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(Mode==(Upper|Lower)) 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[StorageOrderMatch ? jp : ip]++; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp; 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(r==c) 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[ip]++; 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 5012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else if(( (Mode&Lower)==Lower && r>c) || ( (Mode&Upper)==Upper && r<c)) 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(ip,jp); 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k = count[jp]++; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = ip; 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k = count[ip]++; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.innerIndexPtr()[k] = jp; 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.valuePtr()[k] = numext::conj(it.value()); 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<int _SrcMode,int _DstMode,typename MatrixType,int DstOrder> 5172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::StorageIndex>& _dest, const typename MatrixType::StorageIndex* perm) 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 5192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 5212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<Scalar,DstOrder,StorageIndex>& dest(_dest.derived()); 5222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex,Dynamic,1> VectorI; 5232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef evaluator<MatrixType> MatEval; 5242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename evaluator<MatrixType>::InnerIterator MatIterator; 5252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrderMatch = int(SrcOrder) == int(DstOrder), 5292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang DstMode = DstOrder==RowMajor ? (_DstMode==Upper ? Lower : Upper) : _DstMode, 5302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SrcMode = SrcOrder==RowMajor ? (_SrcMode==Upper ? Lower : Upper) : _SrcMode 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 5322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatEval matEval(mat); 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = mat.rows(); 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorI count(size); 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count.setZero(); 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resize(size,size); 5392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex j = 0; j<size; ++j) 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex jp = perm ? perm[j] : j; 5422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(MatIterator it(matEval,j); it; ++it) 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex i = it.index(); 5452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if((int(SrcMode)==int(Lower) && i<j) || (int(SrcMode)==int(Upper) && i>j)) 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex ip = perm ? perm[i] : i; 5492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[0] = 0; 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.resizeNonZeros(dest.outerIndexPtr()[size]); 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index j=0; j<size; ++j) 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath count[j] = dest.outerIndexPtr()[j]; 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex j = 0; j<size; ++j) 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(MatIterator it(matEval,j); it; ++it) 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex i = it.index(); 5652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if((int(SrcMode)==int(Lower) && i<j) || (int(SrcMode)==int(Upper) && i>j)) 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex jp = perm ? perm[j] : j; 5692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex ip = perm? perm[i] : i; 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index k = count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 5722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dest.innerIndexPtr()[k] = int(DstMode)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!StorageOrderMatch) std::swap(ip,jp); 5752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if( ((int(DstMode)==int(Lower) && ip<jp) || (int(DstMode)==int(Upper) && ip>jp))) 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dest.valuePtr()[k] = numext::conj(it.value()); 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.valuePtr()[k] = it.value(); 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// TODO implement twists in a more evaluator friendly fashion 5862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 5882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, int Mode> 5902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct traits<SparseSymmetricPermutationProduct<MatrixType,Mode> > : traits<MatrixType> { 5912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 5922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 5942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType,int Mode> 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SparseSymmetricPermutationProduct 5972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,Mode> > 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 6012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 6022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 6032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = internal::traits<SparseSymmetricPermutationProduct>::RowsAtCompileTime, 6042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = internal::traits<SparseSymmetricPermutationProduct>::ColsAtCompileTime 6052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 6072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> Perm; 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 6092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex,Dynamic,1> VectorI; 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Nested MatrixTypeNested; 6112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::remove_all<MatrixTypeNested>::type NestedExpression; 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(mat), m_perm(perm) 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 6192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const NestedExpression& matrix() const { return m_matrix; } 6212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Perm& perm() const { return m_perm; } 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixTypeNested m_matrix; 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Perm& m_perm; 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 6292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 6302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename DstXprType, typename MatrixType, int Mode, typename Scalar> 6322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct Assignment<DstXprType, SparseSymmetricPermutationProduct<MatrixType,Mode>, internal::assign_op<Scalar,typename MatrixType::Scalar>, Sparse2Sparse> 6332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 6342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseSymmetricPermutationProduct<MatrixType,Mode> SrcXprType; 6352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename DstXprType::StorageIndex DstIndex; 6362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<int Options> 6372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseMatrix<Scalar,Options,DstIndex> &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename MatrixType::Scalar> &) 6382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // internal::permute_symm_to_fullsymm<Mode>(m_matrix,_dest,m_perm.indices().data()); 6402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SparseMatrix<Scalar,(Options&RowMajor)==RowMajor ? ColMajor : RowMajor, DstIndex> tmp; 6412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::permute_symm_to_fullsymm<Mode>(src.matrix(),tmp,src.perm().indices().data()); 6422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dst = tmp; 6432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename DestType,unsigned int DestMode> 6462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(SparseSelfAdjointView<DestType,DestMode>& dst, const SrcXprType &src, const internal::assign_op<Scalar,typename MatrixType::Scalar> &) 6472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::permute_symm_to_symm<Mode,DestMode>(src.matrix(),dst.matrix(),src.perm().indices().data()); 6492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 6512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 6532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H 657