Eigen:: Cast sparse matrix as, specifically, row order or col order

Question

I have in a function internal,

Eigen::SparseMatrix<double> & M;

if (M.IsRowMajor)
     return my_func_template<Eigen::SparseMatrix<double,Eigen::RowMajor>&>(M,M.rows());

However, this does not compile, as the compiler does not believe M is an Eigen::SparseMatrix<double,Eigen::RowMajor>. How do I cast my reference as, specifically, Eigen::SparseMatrix<double,Eigen::RowMajor>, in the type-safe environment of C++11?

---

For example:

typedef Eigen::SparseMatrix<double> Smat;
typedef Eigen::SparseMatrix<double,Eigen::RowMajor> RMSmat;
typedef Eigen::SparseMatrix<double,Eigen::ColMajor> CMSmat;    

enum direction { row, col};

template<class Mat>
vector<double> sum_along_inner(Mat &M){
    vector<double> sums(M.innerSize(),0);
    for(auto i = 0; i < M.outerSize(); i++){
        for(typename M::InnerIterator it(M,i); it;++it){
            sums[i] += it.value();
        }
    }
}
vector<double> sum_along_axis(Smat &M, direction dir){

    // If I could solve this problem, 
    // 
    // I could also function off these if components, 
    // and re-use them for other order-dependent functions I write
    // so that my top level functions are only about 2-4 lines long

    if(dir == direction::row){
        if(M.IsRowMajor)
            return sum_along_inner<RMSmat>((my question) M);

        //else
        RMsmat Mrowmajor = M;
        return sum_along_inner<RMSmat>(Mrowmajor);
    }
    else { 
       if(!M.IsRowMajor)
            return sum_along_inner<CMSmat>(M);

       // else
       CMSmat Mcolmajor = M;
       return sum_along_inner<CMSmat>((my_question) Mcolmajor);
    }
}

And if I do more than just sum_along_axis, then the code complexity in terms of number of lines, readability, etc. is double what it needs to be if only I could solve this problem that I am asking about.

Otherwise, I can't abstract the loop, and I have to repeat it for column major and row major...because I can't just assume I wont call sum_along_axis from a function that hasn't already swapped the major-order from the default Eigen::ColMajor to Eigen::RowMajor...

Further, if I am operating at the order of mb-sized sparse matrices with dimensions too unwieldy to represent in dense matrix form, I am going to notice a major slowdown (which defeats the purpose of using a sparse matrix to begin with) if I don't write composable functions which are order agnostic, and transition the major-order only when needed.

So, unless I solve for this, my line count and/or function count, more or less, starts to go combinatorial.

Answer 1

As I wrote in my first comment M.IsRowMajor will always be false. This is because Eigen::SparseMatrix has always two template arguments, where the second defaults to Eigen::ColMajor

If you want to write a function which accepts both row- and column-major matrices, you need to write something like

template<int mode>
vector<double> sum_along_axis(Eigen::SparseMatrix<double,mode> const &M, direction dir)
    if(dir == direction::row){
        return sum_along_inner<RMSmat>(M); // implicit conversion if necessary
    }
    else { 
        return sum_along_inner<CMSmat>(M); // implicit conversion if necessary
    }
}

You need to rewrite sum_along_inner to accept a const reference to make the implicit conversion work:

template<class Mat>
vector<double> sum_along_inner(Mat const &M){
    vector<double> sums(M.outerSize(),0); // sums needs to have size M.outerSize()
    for(auto i = 0; i < M.outerSize(); i++){
        for(typename M::InnerIterator it(M,i); it;++it){
            sums[i] += it.value();
        }
    }
}

If you want to avoid the conversion from row- to column-major (and vice versa) you should write a function which sums along the outer dimension and decide in your main function which function to call.