[eigen] SparseMatrix with non POD-"Scalar" type (e.g. SparseBlockMatrix) yields unnecessary initialisations of Scalar-type?

[<Daniel.Vollmer@xxxxxx>]

Hello all,

(note: my last attempt to send this message did not seem to come through, so this time I've not attached the example but appended it)

I've been experimenting with Eigen (3.2.0-beta1 at the moment, but I don't think I'm using any of the new features)
to see whether I can use its SparseMatrix type to store block-matrices of fixed size (instead of scalars).

I've attached my experiment that essentially uses a dense 5x5 matrix as the "Scalar" type. Basically, this works (which is already pretty nice) but I have noticed two things which seem sub-optimal, and so I was curious whether I'm missing something or doing it wrong.

1) The Scalar-type used in SparseMatrix seems to need an abs() function in order to (I assume) only iterate over non-zero elements (AmbiVector.h:308 & 335). Oddly enough, my abs-version was needed to compile, but never got called.
It might be advantageous to provide an iterator that iterates over all elements without making such a guarantee, but then the point shouldn't matter once it is compressed.
2) Many times the single-argument constructor is called for the Scalar-type of SparseMatrix in instances where I don't think this is necessary. As this becomes more and more expensive with a more heavy-weight Scalar type, I'd like to avoid this where possible, for example
- SparseMatrix.h:1108 (insertUncompressed): Element is initialised to 0 even though it is later overwritten by the inserted value. Worked fine for me if I removed the = 0.
- SparseMatrix.h:382 (insertBackByOuterInner) and SparseMatrix.h:392 (insertBackByOuterInnerUnordered): Element is initialised to 0 even though it is later overwritten by the inserted value. I hacked this by making a variant of CompressedStorage::append() that only takes an Index and not a value-ref and calling that instead (because again, the value is later overwritten via the reference we return)
- AmbiVector.h:31 (constructor): This seems fair enough (needed for m_zero)..
- AmbiVector.h:171 (setZero): Shouldn't / couldn't this use assignment from m_zero instead?
- AmbiVector.h:197/206/239 (coeffRef): I don't quite understand this, but I'm not sure the element needs to be initialised here either.

To verify the above, just comment the relevant constructor in the sample code (lines 30-35).

Is this approach (using a custom scalar type) still the recommended way to create sparse block-matrices? Can we reduce the overhead of these initialisations or make them optional?

I'd welcome any feedback. :)

Thanks!
    Daniel Vollmer

--------------------------
Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR)
German Aerospace Center
Institute of Aerodynamics and Flow Technology | Lilienthalplatz 7 | 38108 Braunschweig | Germany

main.cpp

//
//  main.cpp
//  SparseMatrixTest
//
//  Created by Daniel Vollmer on 15/03/2013.
//  Copyright (c) 2013 DLR. All rights reserved.
//

#include <iostream>

#include "Eigen/Dense"
#include "Eigen/Sparse"
#include "Eigen/Cholesky"
#include "Eigen/IterativeLinearSolvers"

template<typename B>
class BlockMatrixWrapper : public B
{
public:
	enum
	{
		BlockRowsAtCompileTime = B::RowsAtCompileTime,
		BlockColsAtCompileTime = B::ColsAtCompileTime
	};
	
	BlockMatrixWrapper(void):B() {}
	

	// comment to see who tries to create instances via this constructor
	BlockMatrixWrapper(typename B::Scalar V)
	: B(B::Zero())
	{
		std::cout << "BMI " << V << " "; // shows when we get constructed with BlockMatrixWrapper(0) for example
		B::diagonal().setConstant(V);
	}
	
	// This constructor allows you to construct BlockMatrix from Eigen expressions
	template<typename OtherDerived>
	BlockMatrixWrapper(const Eigen::MatrixBase<OtherDerived>& other)
	: B(other)
	{
	}
	
	// This method allows you to assign Eigen expressions to BlockMatrix
	template<typename OtherDerived>
	BlockMatrixWrapper& operator= (const Eigen::MatrixBase <OtherDerived>& other)
	{
		B::operator=(other);
		return *this;
	}
};

namespace Eigen {
	template<typename B> struct NumTraits<BlockMatrixWrapper<B> >
	: NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
	{
		typedef typename B::Scalar Real;
		typedef B NonInteger;
		typedef B Nested;
		enum {
			IsComplex = 0,
			IsInteger = 0,
			IsSigned = 1,
			RequireInitialization = 1,
			ReadCost = B::RowsAtCompileTime * B::ColsAtCompileTime * NumTraits<typename B::Scalar>::ReadCost,
			AddCost = B::RowsAtCompileTime * B::ColsAtCompileTime * NumTraits<typename B::Scalar>::AddCost,
			MulCost = B::RowsAtCompileTime * B::RowsAtCompileTime * B::ColsAtCompileTime * NumTraits<typename B::Scalar>::MulCost
			+ (B::RowsAtCompileTime-1) * B::RowsAtCompileTime * B::ColsAtCompileTime * NumTraits<typename B::Scalar>::AddCost
		};
	};
}

// needed for InnerIterator :(
template<typename B>
inline typename BlockMatrixWrapper<B>::Scalar abs(const BlockMatrixWrapper<B>& x)
{
	std::cout << "ABS " << x << " ";
	return x.norm();
}

typedef Eigen::MatrixXd DenseMatrix;
typedef BlockMatrixWrapper<Eigen::Matrix<double, 5, 5> > BlockMatrix;
typedef Eigen::SparseMatrix<BlockMatrix> SparseBlockMatrix;
typedef Eigen::Matrix<BlockMatrix, Eigen::Dynamic, 1> BlockVector;

DenseMatrix ExpandSparseBlockMatrix(const SparseBlockMatrix &s)
{
	DenseMatrix result(s.rows() * SparseBlockMatrix::Scalar::BlockRowsAtCompileTime, s.cols() * SparseBlockMatrix::Scalar::BlockColsAtCompileTime);
	
	result.setZero();
	for (SparseBlockMatrix::Index k = 0; k < s.outerSize(); k++)
	{
		for (SparseBlockMatrix::InnerIterator it(s, k); it; ++it)
		{
			result.block<SparseBlockMatrix::Scalar::BlockRowsAtCompileTime, SparseBlockMatrix::Scalar::BlockColsAtCompileTime>(it.row() * SparseBlockMatrix::Scalar::BlockRowsAtCompileTime, it.col() * SparseBlockMatrix::Scalar::BlockColsAtCompileTime) = it.value();
		}
	}
	return result;
}

int main(int argc, const char * argv[])
{
	const int size = 4;
	SparseBlockMatrix sm(size, size);
	BlockMatrix block(BlockMatrix::Zero());
//	BlockVector rhs(1024), soln(1024);
	
	block.diagonal().setOnes();
	for (int i = 0; i < size; i++)
	{
		sm.insert(i, i) = block;
		if (i > 0)
		{
			sm.insert(i, i - 1) = block;
		}
		if (i < size - 1)
		{
			sm.insert(i, i + 1) = block;
		}
//		rhs[i] = block;
	}
//	sm.makeCompressed();
	DenseMatrix dense = ExpandSparseBlockMatrix(sm);
	
	sm = sm * sm; // SparseMatrix doesn't seem to have operator*=
	dense *= dense;
	std::cout << "ErrorNorm: " << (dense - ExpandSparseBlockMatrix(sm)).norm() << std::endl << dense << std::endl;
	
//	Eigen::SimplicialLDLT<SparseBlockMatrix> solver;
//	solver.compute(sm);
//	soln = solver.solve(rhs);
	
	const BlockMatrix &a = sm.sum();
	std::cout << "nonZeros: " << sm.nonZeros() << std::endl << a;
	
    return 0;
}