Search code examples
c++sparse-matrixeigen

Eigen: How to initialize a sparse matrix with some sub sparse matrix


in eigen, we can initialize a matrix or vector with some other matrix or vector like this:

MatrixXf matA(2, 2);
matA << 1, 2, 3, 4;
MatrixXf matB(4, 4);
matB << matA, matA/10, matA/10, matA;
std::cout << matB << std::endl;

what i want to achieve:

SparseMatrix<double> matA(2, 2);
matA.coeffRef(0, 0) = 1;
matA.coeffRef(1, 1) = 1;
SparseMatrix<double> matB(4, 4);
matB << matA, matA/10, matA/10, matA;
std::cout << matB << std::endl;

then i get a matrix like this:

1   0   0.1 0
0   1   0   0.1
0.1 0   1   0
0   0.1 0   0.1

but, it doesn't work for sparse matrix, so does eigen have built-in initializer like this? or i need to write it myself, if so? how?


Solution

  • You cannot have such an initializer because of the storage format. From the manual Sparse matrix manipulations > Block operations:

    However, for performance reasons, writing to a sub-sparse-matrix is much more limited, and currently only contiguous sets of columns (resp. rows) of a column-major (resp. row-major) SparseMatrix are writable. Moreover, this information has to be known at compile-time, leaving out methods such as block(...) and corner*(...).

    The only option you have is to convert everything to dense matrices, use the comma initializer and convert back to sparse.

    #include <iostream>
    #include <Eigen/Sparse>
    
    using namespace Eigen;
    typedef SparseMatrix<double> SparseMatrixXd;
    
    int main()
    {
      SparseMatrixXd matA(2, 2);
      matA.coeffRef(0, 0) = 1;
      matA.coeffRef(1, 1) = 1;
      SparseMatrixXd matB(4, 4);
      MatrixXd matC(4,4);
      matC <<
        MatrixXd(matA),
        MatrixXd(matA)/10,
        MatrixXd(matA)/10,
        MatrixXd(matA);
      matB = matC.sparseView();
      std::cout << matB << std::endl;
    }
    

    Alternatively you can use the unsupported Kronecker product module for this exact example.

    #include <iostream>
    #include <Eigen/Sparse>
    #include <unsupported/Eigen/KroneckerProduct>
    
    using namespace Eigen;
    typedef SparseMatrix<double> SparseMatrixXd;
    
    int main()
    {
      SparseMatrixXd matA(2, 2);
      matA.coeffRef(0, 0) = 1;
      matA.coeffRef(1, 1) = 1;
      SparseMatrixXd matB(4, 4);
      matB =
        kroneckerProduct( (MatrixXd(2,2) << 1,0,0,1).finished(), matA ) +
        kroneckerProduct( (MatrixXd(2,2) << 0,1,1,0).finished(), matA/10);
      std::cout << matB << std::endl;
    }