25 #ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H
26 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
44 template<
typename Lhs,
typename Rhs,
int UpLo>
45 class SparseSelfAdjointTimeDenseProduct;
47 template<
typename Lhs,
typename Rhs,
int UpLo>
48 class DenseTimeSparseSelfAdjointProduct;
52 template<
typename MatrixType,
unsigned int UpLo>
53 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > :
traits<MatrixType> {
56 template<
int SrcUpLo,
int DstUpLo,
typename MatrixType,
int DestOrder>
57 void permute_symm_to_symm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm = 0);
59 template<
int UpLo,
typename MatrixType,
int DestOrder>
60 void permute_symm_to_fullsymm(
const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest,
const typename MatrixType::Index* perm = 0);
65 :
public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
69 typedef typename MatrixType::Scalar
Scalar;
70 typedef typename MatrixType::Index
Index;
88 template<
typename OtherDerived>
96 template<
typename OtherDerived>
friend
111 template<
typename DerivedU>
117 internal::permute_symm_to_fullsymm<UpLo>(
m_matrix, _dest);
124 internal::permute_symm_to_fullsymm<UpLo>(
m_matrix, tmp);
134 template<
typename SrcMatrixType,
int SrcUpLo>
137 permutedMatrix.
evalTo(*
this);
148 template<
typename SrcMatrixType,
unsigned int SrcUpLo>
170 template<
typename Derived>
171 template<
unsigned int UpLo>
177 template<
typename Derived>
178 template<
unsigned int UpLo>
188 template<
typename MatrixType,
unsigned int UpLo>
189 template<
typename DerivedU>
195 m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
207 template<
typename Lhs,
typename Rhs,
int UpLo>
209 :
traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
211 typedef Dense StorageKind;
215 template<
typename Lhs,
typename Rhs,
int UpLo>
217 :
public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
229 typedef typename internal::remove_all<Lhs>::type _Lhs;
230 typedef typename internal::remove_all<Rhs>::type _Rhs;
231 typedef typename _Lhs::InnerIterator LhsInnerIterator;
236 || ( (UpLo&
Upper) && !LhsIsRowMajor)
237 || ( (UpLo&
Lower) && LhsIsRowMajor),
238 ProcessSecondHalf = !ProcessFirstHalf
242 LhsInnerIterator i(
m_lhs,j);
243 if (ProcessSecondHalf)
245 while (i && i.index()<j) ++i;
246 if(i && i.index()==j)
248 dest.row(j) += i.value() *
m_rhs.row(j);
252 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
254 Index a = LhsIsRowMajor ? j : i.index();
255 Index b = LhsIsRowMajor ? i.index() : j;
256 typename Lhs::Scalar v = i.value();
257 dest.row(a) += (v) *
m_rhs.row(b);
260 if (ProcessFirstHalf && i && (i.index()==j))
261 dest.row(j) += i.value() *
m_rhs.row(j);
270 template<
typename Lhs,
typename Rhs,
int UpLo>
271 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
272 :
traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
276 template<
typename Lhs,
typename Rhs,
int UpLo>
278 :
public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
300 template<
typename MatrixType,
int UpLo>
301 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > :
traits<MatrixType> {
304 template<
int UpLo,
typename MatrixType,
int DestOrder>
307 typedef typename MatrixType::Index Index;
308 typedef typename MatrixType::Scalar Scalar;
314 StorageOrderMatch =
int(Dest::IsRowMajor) ==
int(MatrixType::IsRowMajor)
317 Index size = mat.rows();
321 dest.resize(size,size);
322 for(Index j = 0; j<size; ++j)
324 Index jp = perm ? perm[j] : j;
325 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
327 Index i = it.index();
330 Index ip = perm ? perm[i] : i;
332 count[StorageOrderMatch ? jp : ip]++;
335 else if(( UpLo==
Lower && r>c) || ( UpLo==
Upper && r<c))
342 Index nnz = count.sum();
345 dest.resizeNonZeros(nnz);
346 dest.outerIndexPtr()[0] = 0;
347 for(Index j=0; j<size; ++j)
348 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
349 for(Index j=0; j<size; ++j)
350 count[j] = dest.outerIndexPtr()[j];
353 for(Index j = 0; j<size; ++j)
355 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
357 Index i = it.index();
361 Index jp = perm ? perm[j] : j;
362 Index ip = perm ? perm[i] : i;
366 Index k = count[StorageOrderMatch ? jp : ip]++;
367 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
368 dest.valuePtr()[k] = it.value();
372 Index k = count[ip]++;
373 dest.innerIndexPtr()[k] = ip;
374 dest.valuePtr()[k] = it.value();
376 else if(( (UpLo&
Lower)==Lower && r>c) || ( (UpLo&
Upper)==Upper && r<c))
378 if(!StorageOrderMatch)
380 Index k = count[jp]++;
381 dest.innerIndexPtr()[k] = ip;
382 dest.valuePtr()[k] = it.value();
384 dest.innerIndexPtr()[k] = jp;
391 template<
int _SrcUpLo,
int _DstUpLo,
typename MatrixType,
int DstOrder>
394 typedef typename MatrixType::Index Index;
395 typedef typename MatrixType::Scalar Scalar;
400 StorageOrderMatch =
int(SrcOrder) ==
int(DstOrder),
405 Index size = mat.rows();
408 dest.resize(size,size);
409 for(Index j = 0; j<size; ++j)
411 Index jp = perm ? perm[j] : j;
412 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
414 Index i = it.index();
415 if((
int(SrcUpLo)==
int(
Lower) && i<j) || (
int(SrcUpLo)==
int(
Upper) && i>j))
418 Index ip = perm ? perm[i] : i;
419 count[
int(DstUpLo)==
int(
Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
422 dest.outerIndexPtr()[0] = 0;
423 for(Index j=0; j<size; ++j)
424 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
425 dest.resizeNonZeros(dest.outerIndexPtr()[size]);
426 for(Index j=0; j<size; ++j)
427 count[j] = dest.outerIndexPtr()[j];
429 for(Index j = 0; j<size; ++j)
432 for(
typename MatrixType::InnerIterator it(mat,j); it; ++it)
434 Index i = it.index();
435 if((
int(SrcUpLo)==
int(
Lower) && i<j) || (
int(SrcUpLo)==
int(
Upper) && i>j))
438 Index jp = perm ? perm[j] : j;
439 Index ip = perm? perm[i] : i;
441 Index k = count[
int(DstUpLo)==
int(
Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
442 dest.innerIndexPtr()[k] =
int(DstUpLo)==
int(
Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
444 if(!StorageOrderMatch) std::swap(ip,jp);
445 if( ((
int(DstUpLo)==
int(
Lower) && ip<jp) || (
int(DstUpLo)==
int(
Upper) && ip>jp)))
446 dest.valuePtr()[k] =
conj(it.value());
448 dest.valuePtr()[k] = it.value();
455 template<
typename MatrixType,
int UpLo>
457 :
public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
460 typedef typename MatrixType::Scalar
Scalar;
461 typedef typename MatrixType::Index
Index;
476 template<
typename DestScalar,
int Options,
typename DstIndex>
495 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H