00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H
00026 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
00027
00028 namespace Eigen {
00029
00044 template<typename Lhs, typename Rhs, int UpLo>
00045 class SparseSelfAdjointTimeDenseProduct;
00046
00047 template<typename Lhs, typename Rhs, int UpLo>
00048 class DenseTimeSparseSelfAdjointProduct;
00049
00050 namespace internal {
00051
00052 template<typename MatrixType, unsigned int UpLo>
00053 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
00054 };
00055
00056 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
00057 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00058
00059 template<int UpLo,typename MatrixType,int DestOrder>
00060 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00061
00062 }
00063
00064 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
00065 : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
00066 {
00067 public:
00068
00069 typedef typename MatrixType::Scalar Scalar;
00070 typedef typename MatrixType::Index Index;
00071 typedef Matrix<Index,Dynamic,1> VectorI;
00072 typedef typename MatrixType::Nested MatrixTypeNested;
00073 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00074
00075 inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
00076 {
00077 eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
00078 }
00079
00080 inline Index rows() const { return m_matrix.rows(); }
00081 inline Index cols() const { return m_matrix.cols(); }
00082
00084 const _MatrixTypeNested& matrix() const { return m_matrix; }
00085 _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
00086
00088 template<typename OtherDerived>
00089 SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
00090 operator*(const MatrixBase<OtherDerived>& rhs) const
00091 {
00092 return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
00093 }
00094
00096 template<typename OtherDerived> friend
00097 DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
00098 operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
00099 {
00100 return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
00101 }
00102
00111 template<typename DerivedU>
00112 SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00113
00115 template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const
00116 {
00117 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
00118 }
00119
00120 template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const
00121 {
00122
00123 SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols());
00124 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
00125 _dest = tmp;
00126 }
00127
00129 SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
00130 {
00131 return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);
00132 }
00133
00134 template<typename SrcMatrixType,int SrcUpLo>
00135 SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix)
00136 {
00137 permutedMatrix.evalTo(*this);
00138 return *this;
00139 }
00140
00141
00142 SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src)
00143 {
00144 PermutationMatrix<Dynamic> pnull;
00145 return *this = src.twistedBy(pnull);
00146 }
00147
00148 template<typename SrcMatrixType,unsigned int SrcUpLo>
00149 SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src)
00150 {
00151 PermutationMatrix<Dynamic> pnull;
00152 return *this = src.twistedBy(pnull);
00153 }
00154
00155
00156
00157
00158
00159 protected:
00160
00161 typename MatrixType::Nested m_matrix;
00162 mutable VectorI m_countPerRow;
00163 mutable VectorI m_countPerCol;
00164 };
00165
00166
00167
00168
00169
00170 template<typename Derived>
00171 template<unsigned int UpLo>
00172 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
00173 {
00174 return derived();
00175 }
00176
00177 template<typename Derived>
00178 template<unsigned int UpLo>
00179 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
00180 {
00181 return derived();
00182 }
00183
00184
00185
00186
00187
00188 template<typename MatrixType, unsigned int UpLo>
00189 template<typename DerivedU>
00190 SparseSelfAdjointView<MatrixType,UpLo>&
00191 SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha)
00192 {
00193 SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
00194 if(alpha==Scalar(0))
00195 m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
00196 else
00197 m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
00198
00199 return *this;
00200 }
00201
00202
00203
00204
00205
00206 namespace internal {
00207 template<typename Lhs, typename Rhs, int UpLo>
00208 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
00209 : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00210 {
00211 typedef Dense StorageKind;
00212 };
00213 }
00214
00215 template<typename Lhs, typename Rhs, int UpLo>
00216 class SparseSelfAdjointTimeDenseProduct
00217 : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00218 {
00219 public:
00220 EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
00221
00222 SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00223 {}
00224
00225 template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
00226 {
00227
00228 eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
00229 typedef typename internal::remove_all<Lhs>::type _Lhs;
00230 typedef typename internal::remove_all<Rhs>::type _Rhs;
00231 typedef typename _Lhs::InnerIterator LhsInnerIterator;
00232 enum {
00233 LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
00234 ProcessFirstHalf =
00235 ((UpLo&(Upper|Lower))==(Upper|Lower))
00236 || ( (UpLo&Upper) && !LhsIsRowMajor)
00237 || ( (UpLo&Lower) && LhsIsRowMajor),
00238 ProcessSecondHalf = !ProcessFirstHalf
00239 };
00240 for (Index j=0; j<m_lhs.outerSize(); ++j)
00241 {
00242 LhsInnerIterator i(m_lhs,j);
00243 if (ProcessSecondHalf)
00244 {
00245 while (i && i.index()<j) ++i;
00246 if(i && i.index()==j)
00247 {
00248 dest.row(j) += i.value() * m_rhs.row(j);
00249 ++i;
00250 }
00251 }
00252 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
00253 {
00254 Index a = LhsIsRowMajor ? j : i.index();
00255 Index b = LhsIsRowMajor ? i.index() : j;
00256 typename Lhs::Scalar v = i.value();
00257 dest.row(a) += (v) * m_rhs.row(b);
00258 dest.row(b) += internal::conj(v) * m_rhs.row(a);
00259 }
00260 if (ProcessFirstHalf && i && (i.index()==j))
00261 dest.row(j) += i.value() * m_rhs.row(j);
00262 }
00263 }
00264
00265 private:
00266 SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
00267 };
00268
00269 namespace internal {
00270 template<typename Lhs, typename Rhs, int UpLo>
00271 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
00272 : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00273 {};
00274 }
00275
00276 template<typename Lhs, typename Rhs, int UpLo>
00277 class DenseTimeSparseSelfAdjointProduct
00278 : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00279 {
00280 public:
00281 EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
00282
00283 DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00284 {}
00285
00286 template<typename Dest> void scaleAndAddTo(Dest& , Scalar ) const
00287 {
00288
00289 }
00290
00291 private:
00292 DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
00293 };
00294
00295
00296
00297
00298 namespace internal {
00299
00300 template<typename MatrixType, int UpLo>
00301 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
00302 };
00303
00304 template<int UpLo,typename MatrixType,int DestOrder>
00305 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00306 {
00307 typedef typename MatrixType::Index Index;
00308 typedef typename MatrixType::Scalar Scalar;
00309 typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
00310 typedef Matrix<Index,Dynamic,1> VectorI;
00311
00312 Dest& dest(_dest.derived());
00313 enum {
00314 StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
00315 };
00316
00317 Index size = mat.rows();
00318 VectorI count;
00319 count.resize(size);
00320 count.setZero();
00321 dest.resize(size,size);
00322 for(Index j = 0; j<size; ++j)
00323 {
00324 Index jp = perm ? perm[j] : j;
00325 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00326 {
00327 Index i = it.index();
00328 Index r = it.row();
00329 Index c = it.col();
00330 Index ip = perm ? perm[i] : i;
00331 if(UpLo==(Upper|Lower))
00332 count[StorageOrderMatch ? jp : ip]++;
00333 else if(r==c)
00334 count[ip]++;
00335 else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c))
00336 {
00337 count[ip]++;
00338 count[jp]++;
00339 }
00340 }
00341 }
00342 Index nnz = count.sum();
00343
00344
00345 dest.resizeNonZeros(nnz);
00346 dest.outerIndexPtr()[0] = 0;
00347 for(Index j=0; j<size; ++j)
00348 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
00349 for(Index j=0; j<size; ++j)
00350 count[j] = dest.outerIndexPtr()[j];
00351
00352
00353 for(Index j = 0; j<size; ++j)
00354 {
00355 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00356 {
00357 Index i = it.index();
00358 Index r = it.row();
00359 Index c = it.col();
00360
00361 Index jp = perm ? perm[j] : j;
00362 Index ip = perm ? perm[i] : i;
00363
00364 if(UpLo==(Upper|Lower))
00365 {
00366 Index k = count[StorageOrderMatch ? jp : ip]++;
00367 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
00368 dest.valuePtr()[k] = it.value();
00369 }
00370 else if(r==c)
00371 {
00372 Index k = count[ip]++;
00373 dest.innerIndexPtr()[k] = ip;
00374 dest.valuePtr()[k] = it.value();
00375 }
00376 else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c))
00377 {
00378 if(!StorageOrderMatch)
00379 std::swap(ip,jp);
00380 Index k = count[jp]++;
00381 dest.innerIndexPtr()[k] = ip;
00382 dest.valuePtr()[k] = it.value();
00383 k = count[ip]++;
00384 dest.innerIndexPtr()[k] = jp;
00385 dest.valuePtr()[k] = internal::conj(it.value());
00386 }
00387 }
00388 }
00389 }
00390
00391 template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
00392 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00393 {
00394 typedef typename MatrixType::Index Index;
00395 typedef typename MatrixType::Scalar Scalar;
00396 SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived());
00397 typedef Matrix<Index,Dynamic,1> VectorI;
00398 enum {
00399 SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor,
00400 StorageOrderMatch = int(SrcOrder) == int(DstOrder),
00401 DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo,
00402 SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
00403 };
00404
00405 Index size = mat.rows();
00406 VectorI count(size);
00407 count.setZero();
00408 dest.resize(size,size);
00409 for(Index j = 0; j<size; ++j)
00410 {
00411 Index jp = perm ? perm[j] : j;
00412 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00413 {
00414 Index i = it.index();
00415 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
00416 continue;
00417
00418 Index ip = perm ? perm[i] : i;
00419 count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00420 }
00421 }
00422 dest.outerIndexPtr()[0] = 0;
00423 for(Index j=0; j<size; ++j)
00424 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
00425 dest.resizeNonZeros(dest.outerIndexPtr()[size]);
00426 for(Index j=0; j<size; ++j)
00427 count[j] = dest.outerIndexPtr()[j];
00428
00429 for(Index j = 0; j<size; ++j)
00430 {
00431
00432 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00433 {
00434 Index i = it.index();
00435 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
00436 continue;
00437
00438 Index jp = perm ? perm[j] : j;
00439 Index ip = perm? perm[i] : i;
00440
00441 Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00442 dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
00443
00444 if(!StorageOrderMatch) std::swap(ip,jp);
00445 if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp)))
00446 dest.valuePtr()[k] = conj(it.value());
00447 else
00448 dest.valuePtr()[k] = it.value();
00449 }
00450 }
00451 }
00452
00453 }
00454
00455 template<typename MatrixType,int UpLo>
00456 class SparseSymmetricPermutationProduct
00457 : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
00458 {
00459 public:
00460 typedef typename MatrixType::Scalar Scalar;
00461 typedef typename MatrixType::Index Index;
00462 protected:
00463 typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm;
00464 public:
00465 typedef Matrix<Index,Dynamic,1> VectorI;
00466 typedef typename MatrixType::Nested MatrixTypeNested;
00467 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00468
00469 SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
00470 : m_matrix(mat), m_perm(perm)
00471 {}
00472
00473 inline Index rows() const { return m_matrix.rows(); }
00474 inline Index cols() const { return m_matrix.cols(); }
00475
00476 template<typename DestScalar, int Options, typename DstIndex>
00477 void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const
00478 {
00479 internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
00480 }
00481
00482 template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
00483 {
00484 internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
00485 }
00486
00487 protected:
00488 MatrixTypeNested m_matrix;
00489 const Perm& m_perm;
00490
00491 };
00492
00493 }
00494
00495 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H