GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Sparse module: add support for sparse selfadjoint * dense

Browse Source
This commit is contained in:
Gael Guennebaud 2009-01-15 18:52:14 +00:00
parent 9a4b7998cf
commit ccdcebcf03
4 changed files with 87 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Sparse/SparseMatrixBase.h
View File

@ -88,7 +88,7 @@ template<typename Derived> class SparseMatrixBase
/** \internal the return type of MatrixBase::imag() */
typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
/** \internal the return type of MatrixBase::adjoint() */
typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> >
typedef SparseTranspose</*NestByValue<*/typename ei_cleantype<ConjugateReturnType>::type> /*>*/
AdjointReturnType;
#ifndef EIGEN_PARSED_BY_DOXYGEN
@ -322,7 +322,7 @@ template<typename Derived> class SparseMatrixBase
SparseTranspose<Derived> transpose() { return derived(); }
const SparseTranspose<Derived> transpose() const { return derived(); }
// void transposeInPlace();
// const AdjointReturnType adjoint() const;
const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
SparseInnerVector<Derived> innerVector(int outer);
const SparseInnerVector<Derived> innerVector(int outer) const;

47
Eigen/src/Sparse/SparseProduct.h
View File

@ -294,17 +294,60 @@ inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs
}
// dense = sparse * dense
// template<typename Derived>
// template<typename Lhs, typename Rhs>
// Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
// {
// typedef typename ei_cleantype<Lhs>::type _Lhs;
// typedef typename _Lhs::InnerIterator LhsInnerIterator;
// enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
// derived().setZero();
// for (int j=0; j<product.lhs().outerSize(); ++j)
// for (LhsInnerIterator i(product.lhs(),j); i; ++i)
// derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
// return derived();
// }
template<typename Derived>
template<typename Lhs, typename Rhs>
Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
{
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename _Lhs::InnerIterator LhsInnerIterator;
enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
enum {
LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
LhsIsSelfAdjoint = (_Lhs::Flags&SelfAdjointBit)==SelfAdjointBit,
ProcessFirstHalf = LhsIsSelfAdjoint
&& ( ((_Lhs::Flags&(UpperTriangularBit|LowerTriangularBit))==0)
|| ( (_Lhs::Flags&UpperTriangularBit) && !LhsIsRowMajor)
|| ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ),
ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
};
derived().setZero();
for (int j=0; j<product.lhs().outerSize(); ++j)
for (LhsInnerIterator i(product.lhs(),j); i; ++i)
{
LhsInnerIterator i(product.lhs(),j);
if (ProcessSecondHalf && i && (i.index()==j))
{
derived().row(j) += i.value() * product.rhs().row(j);
++i;
}
for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
{
if (LhsIsSelfAdjoint)
{
int a = LhsIsRowMajor ? j : i.index();
int b = LhsIsRowMajor ? i.index() : j;
Scalar v = i.value();
derived().row(a) += (v) * product.rhs().row(b);
derived().row(b) += ei_conj(v) * product.rhs().row(a);
}
else
derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
}
if (ProcessFirstHalf && i && (i.index()==j))
derived().row(j) += i.value() * product.rhs().row(j);
}
return derived();
}

7
test/sparse.h
View File

@ -40,7 +40,8 @@
enum {
ForceNonZeroDiag = 1,
MakeLowerTriangular = 2,
MakeUpperTriangular = 4
MakeUpperTriangular = 4,
ForceRealDiag = 8
};
/* Initializes both a sparse and dense matrix with same random values,
@ -73,6 +74,10 @@ initSparse(double density,
v = Scalar(0);
else if ((flags & MakeUpperTriangular) && j<i)
v = Scalar(0);
if ((flags&ForceRealDiag) && (i==j))
v = ei_real(v);
if (v!=Scalar(0))
{
sparseMat.fill(i,j) = v;

33
test/sparse_basic.cpp
View File

@ -269,6 +269,39 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
}
// test self adjoint products
{
DenseMatrix b = DenseMatrix::Random(rows, rows);
DenseMatrix x = DenseMatrix::Random(rows, rows);
DenseMatrix refX = DenseMatrix::Random(rows, rows);
DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
DenseMatrix refS = DenseMatrix::Zero(rows, rows);
SparseMatrix<Scalar> mUp(rows, rows);
SparseMatrix<Scalar> mLo(rows, rows);
SparseMatrix<Scalar> mS(rows, rows);
do {
initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
} while (refUp.isZero());
refLo = refUp.transpose().conjugate();
mLo = mUp.transpose().conjugate();
refS = refUp + refLo;
refS.diagonal() *= 0.5;
mS = mUp + mLo;
for (int k=0; k<mS.outerSize(); ++k)
for (typename SparseMatrix<Scalar>::InnerIterator it(mS,k); it; ++it)
if (it.index() == k)
it.valueRef() *= 0.5;
VERIFY_IS_APPROX(refS.adjoint(), refS);
VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
VERIFY_IS_APPROX(mS, refS);
VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
}
}
void test_sparse_basic()