define and filling a sparse matrix using Eigen Library in C++

Question

I am trying to build a spars Matrix using a Eigen or Armadillo library in C++ to solve a system of linear equations Ax=b. A is the coefficient matrix with a dimension of n*n, and B is a vector of right hand side with a dimension of n the Spars Matrix A is like this, see the figure

I had a look though the Eigen document but I have a problem with defining and filling the Spars Matrix in C++.

could you please give me an example code to define the spars matrix and how to fill the values into the matrix using Eigen library in c++?

consider for example a simple spars matrix A:

1 2 0 0

0 3 0 0

0 0 4 5

0 0 6 7

int main()
{
     SparseMatrix<double> A;

     // fill the A matrix ????

     VectorXd b, x;
     SparseCholesky<SparseMatrix<double> > solver;
     solver.compute(A);

     x = solver.solve(b);
     return 0;
} 

Answer 1

#include <vector>
#include <iostream>                                  
#include <Eigen/Dense>                               
#include <Eigen/Sparse>                              
#include <Eigen/Core>
#include <cstdlib>

using namespace Eigen;
using namespace std;

int main()
{
    double L = 5; // Length
    const int N = 120; // No of cells
    double L_cell = L / N;
    double k = 100; // Thermal Conductivity
    double T_A = 100.;
    double T_B = 200.;
    double S = 1000.;

    Vector<double, N> d, D, A, aL, aR, aP, S_u, S_p;
    vector<double> xp;


    xp.push_back((0 + L_cell) / 2.0);
    double xm = xp[0];

    for (int i = 0; i < N - 1; i++)
    {
        xm = xm + L_cell;
        xp.push_back(xm);
    }

    for (int i = 0; i < N; i++)
    {

        A(i) = .1;
        d(i) = L_cell;
        D(i) = k / d(i);
    }

    aL(0) = 0;
    aR(0) = D(0) * A(0);
    S_p(0) = -2 * D(0) * A(0);
    aP(0) = aL(0) + aR(0) - S_p(0);
    S_u(0) = 2 * D(0) * A(0) * T_A + S * L_cell * A(0);

    for (int i = 1; i < N - 1; i++)
    {
        aL(i) = D(i) * A(i);
        aR(i) = D(i) * A(i);
        S_p(i) = 0;
        aP(i) = aL(i) + aR(i) - S_p(i);
        S_u(i) = S * A(i) * L_cell;
    }

    aL(N - 1) = D(N - 1) * A(N - 1);
    aR(N - 1) = 0;
    S_p(N - 1) = -2 * D(N - 1) * A(N - 1);
    aP(N - 1) = aL(N - 1) + aR(N - 1) - S_p(N - 1);
    S_u(N - 1) = 2 * D(N - 1) * A(N - 1) * T_B + S * L_cell * A(N - 1);

    typedef Eigen::Triplet<double> T;
    std::vector<T> tripletList;
    tripletList.reserve(N * 3);

    Matrix<double, N, 3> v; // v is declared here

    v << (-1) * aL, aP, (-1) * aR;

    for (int i = 0, j = 0; i < N && j < N; i++, j++)
    {
        tripletList.push_back(T(i, j, v(i, 1)));

        if (i + 1 < N && j + 1 < N)
        {
            tripletList.push_back(T(i + 1, j, v(i + 1, 0)));
            tripletList.push_back(T(i, j + 1, v(i, 2)));
        }
    }
    SparseMatrix<double> coeff(N, N);
 
    coeff.setFromTriplets(tripletList.begin(), tripletList.end());

    
    SimplicialLDLT<SparseMatrix<double> > solver;
    solver.compute(coeff);

    if (solver.info() != Success) {
        cout << "decomposition failed" << endl;
        return;
    }

    Vector<double, N> temparature;

    temparature = solver.solve(S_u);

    if (solver.info() != Success)
    {
        cout << "Solving failed" << endl;
        return;
    }

    vector<double> Te = {}, x = {};

    Te.push_back(T_A);
    x.push_back(0);

    for (int i = 0; i < N; i++)
    {
        Te.push_back(temparature(i));
        x.push_back(xp[i]);
    }

    Te.push_back(T_B);
    x.push_back(L);
    for (int i = 0; i < N + 2; i++)
    {
        cout << x[i] << " " << Te[i] << endl;
    }
    
    return 0;
}

Here is a full code of a solution to numerical problem which uses SparseMatrix. Look at the matrix v. It has the values of all the nonzero elements of coeff matrix yet to be defined. In the next loop I made a series of tripletList.push_back(...) adding a triplet consisting of row and column index and corresponding value taken from v for each non-zero element of coeff. Now declare a Sparse Matrix coeff with appropriate size and use the method setFromTriplets (documentation) to set its non-zero elements from tripletList triplets.

Answer 2

Since this question also asks about Armadillo, here is the corresponding Armadillo-based code. Best to use Armadillo version 9.100+ or later, and link with SuperLU.

#include <armadillo>

using namespace arma; 

int main()
{
  sp_mat A(4,4);   // don't need to explicitly reserve the number of non-zeros

  // fill with direct element access
  A(0,0) = 1.0;
  A(0,1) = 2.0;
  A(1,1) = 3.0;
  A(2,2) = 4.0;
  A(2,3) = 5.0;
  A(3,2) = 6.0;
  A(3,3) = 7.0;  // etc

  // or load the sparse matrix from a text file with the data stored in coord format
  sp_mat AA;
  AA.load("my_sparse_matrix.txt", coord_ascii)

  vec b; // ... fill b here ...

  vec x = spsolve(A,b);  // solve sparse system

  return 0;
}

See also the documentation for SpMat, element access, .load(), spsolve().

The coord file format is simple. It stores non-zeros values. Each line contains:

row col value

The row and column counts start at zero. Example:

  0 0 1.0
  0 1 2.0
  1 1 3.0
  2 2 4.0
  2 3 5.0
  3 2 6.0
  3 3 7.0
  1000 2000 9.0

Values not explicitly listed are assumed to be zero.

Answer 3

The sparse matrix could be filled with the values mentioned in the post by using the .coeffRef() member function, as shown in this routine:

SparseMatrix<double> fillMatrix() {
  int N = 4;
  int M = 4;
  SparseMatrix<double> m1(N,M);
  m1.reserve(VectorXi::Constant(M, 4)); // 4: estimated number of non-zero enties per column
  m1.coeffRef(0,0) = 1;
  m1.coeffRef(0,1) = 2.;
  m1.coeffRef(1,1) = 3.;
  m1.coeffRef(2,2) = 4.;
  m1.coeffRef(2,3) = 5.;
  m1.coeffRef(3,2) = 6.;
  m1.coeffRef(3,3) = 7.;
  m1.makeCompressed();
  return m1;
}

However, the SparseCholesky module (SimplicialCholesky<SparseMatrix<double> >) won't work in this case because the matrix is not Hermitian. The system could be solved with a LU or BiCGStab solver. Also note that sizes ofx and b need to be defined: VectorXd b(A.rows()), x(A.cols());

In case of larger sparse matrices you may also want to look at the .reserve() function in order to allocate memory before filling the elements. The .reserve() function can be used to provide an estimate of the number of non-zero entries per column (or row, depending on the storage order. The default is comumn-major). In the example above that estimate is 4, but it does not make sense in such a small matrix. The documentation states that it is preferable to overestimate the number of non-zeros per column.