Improve Performance of repmat for sparse matrices

[paul0noah]

I recently wanted to use the repmat function for sparse matrices and encountered that the function is very slow. Additionally, the function does not work for row-major format.
I came up with a faster implementation which also supports row-major type, so maybe you want to integrate that one into the project :)

template <typename T, int majorType>
void repmat(const Eigen::SparseMatrix<T, majorType> &A, const int r, const int c, Eigen::SparseMatrix<T, majorType> B) {
    assert(r>0);
    assert(c>0);
    B.resize(r*A.rows(), c*A.cols());
    std::vector<Eigen::Triplet<T>> b;
    b.reserve(r*c*A.nonZeros());
    
    for(int i = 0; i < r; i++) {
        for(int j = 0; j < c; j++) {
            // loop outer level
            for (int k = 0; k < A.outerSize(); ++k) {
                // loop inner level
                for (typename Eigen::SparseMatrix<T, majorType>::InnerIterator it(A,k); it; ++it) {
                    Eigen::Triplet<T> triplet(i * A.rows() + it.row(), j * A.cols() + it.col(), it.value());
                    b.push_back(triplet);
                }
            }
        }
    }
    
    B.setFromTriplets(b.begin(), b.end());
}

[alecjacobson]

It's probably a good idea to switch to this. Want to submit a PR?

I remember (perhaps many Eigens ago) that setFromTriplets was much slower than getting into the guts of the output format. That may no longer be true and is always a pain because the rowmajor/colmajor combinations need to be handled more carefully.

[jdumas]

setFromTriplet does two pass on the data, using a bucket sort to handle duplicate/sort indices/transpose in place. If you know the input indices are ordered in advance, you will save on the bucket sort by filing the sparse entries manually (provided you reserved the correct number of nonzeros per rows/cols). Overall I'd say the perf cost of setFromTriplet is about 2x due to the additional pass on the data/memalloc (but it nicely handles all the rowmajor/colmajor headaches).

[paul0noah]

So you'd say it is preferably to allocate the rows and cols before (which is quite hard for general non-zero patterns) and then fill the matrix somehow in order?
i guess from a code-understandability my solution is maybe nicer and less error prone but i totally understand if you don't like it because of the additional sorting.

[alecjacobson]

If input and output are the same w.r.t. RowMajor/ColMajor, then it should be straight-forward to write as repeated loops over the input inserting directly along rows/columns. When they're different it's a bit annoying and maybe at that point it's worth bailing into setfromtriplets.

[alecjacobson]

Actually, I take it back. Even in the case they're different it should be fine. It's just a matter of sprinkling `if RowMajor` in all the right places. In any case, this is a performance optimization that we could do after accepting a PR that at least worked _correctly_ in all cases.

…

[jdumas]

Agreed. Ideally we should have some unit tests to make sure we're not breaking anything when updating the implementation...

[paul0noah]

I just wanted to push a unit test and the improved code but i unfortunately have no access rights.

[jdumas]

Please use a fork in order to create a PR.