Serializing Decomposed Matrix from eigen (SparseLU object)

Question

I am trying to solve Ax=b where the matrix A can be big close to 1M x 1M in size, is sparse and symmetric but might not be positively defined.

Problem is that it might take a long time to compute the decomposition using the sparseLU object in eigen and it will be idea for one to store the sparseLU matrix instead of the original matrix such that whenever we perform similar operation using the same matrix A, we can speed things up by not needing to re-compute the

A quick search on stackoverflow and google returned this, this and this for sparse matrix for the serialization of eigen matrix. However, I am not sure if the same code can be apply for the sparseLU object.

Maybe I should rephrase my question:

How can I store the decomposed matrix into a file?

The current methods all focus on storing the original matrix but I want to store the decomposed matrix. Is there any way to do that? Thank you.

Answer 1

The following example should help you implement your own serialization.

EDIT Changed the example to answer the reworded question.

#include <Eigen/Dense>
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/SparseLU>
#include <iostream>
#include <fstream>

using namespace Eigen;
using namespace std;

typedef Triplet<int> Trip;

template <typename T, int whatever, typename IND>
void Serialize(SparseMatrix<T, whatever, IND>& m) {
    std::vector<Trip> res;
    int sz = m.nonZeros();
    m.makeCompressed();

    fstream writeFile;
    writeFile.open("matrix", ios::binary | ios::out);

    if(writeFile.is_open())
    {
        IND rows, cols, nnzs, outS, innS;
        rows = m.rows()     ;
        cols = m.cols()     ;
        nnzs = m.nonZeros() ;
        outS = m.outerSize();
        innS = m.innerSize();

        writeFile.write((const char *)&(rows), sizeof(IND));
        writeFile.write((const char *)&(cols), sizeof(IND));
        writeFile.write((const char *)&(nnzs), sizeof(IND));
        writeFile.write((const char *)&(outS), sizeof(IND));
        writeFile.write((const char *)&(innS), sizeof(IND));

        writeFile.write((const char *)(m.valuePtr()),       sizeof(T  ) * m.nonZeros());
        writeFile.write((const char *)(m.outerIndexPtr()),  sizeof(IND) * m.outerSize());
        writeFile.write((const char *)(m.innerIndexPtr()),  sizeof(IND) * m.nonZeros());

        writeFile.close();
    }
}

template <typename T, int whatever, typename IND>
void Deserialize(SparseMatrix<T, whatever, IND>& m) {
    fstream readFile;
    readFile.open("matrix", ios::binary | ios::in);
    if(readFile.is_open())
    {
        IND rows, cols, nnz, inSz, outSz;
        readFile.read((char*)&rows , sizeof(IND));
        readFile.read((char*)&cols , sizeof(IND));
        readFile.read((char*)&nnz  , sizeof(IND));
        readFile.read((char*)&outSz, sizeof(IND));
        readFile.read((char*)&inSz , sizeof(IND));

        m.resize(rows, cols);
        m.makeCompressed();
        m.resizeNonZeros(nnz);

        readFile.read((char*)(m.valuePtr())     , sizeof(T  ) * nnz  );
        readFile.read((char*)(m.outerIndexPtr()), sizeof(IND) * outSz);
        readFile.read((char*)(m.innerIndexPtr()), sizeof(IND) * nnz );

        m.finalize();
        readFile.close();

    } // file is open
}


int main(int argc, char *argv[]){
    int rows, cols;
    rows = cols = 6;
    SparseMatrix<double> A(rows,cols), B;

    std::vector<Trip> trp, tmp;

    trp.push_back(Trip(0, 0, rand()));
    trp.push_back(Trip(1, 1, rand()));
    trp.push_back(Trip(2, 2, rand()));
    trp.push_back(Trip(3, 3, rand()));
    trp.push_back(Trip(4, 4, rand()));
    trp.push_back(Trip(5, 5, rand()));
    trp.push_back(Trip(2, 4, rand()));
    trp.push_back(Trip(3, 1, rand()));

    A.setFromTriplets(trp.begin(), trp.end());
    cout << A.nonZeros() << endl;   // Prints 8
    cout << A.size() << endl;       // Prints 36
    cout << A << endl;              // Prints the matrix along with the sparse matrix stuff

    Serialize(A);

    Deserialize(B);

    cout << B.nonZeros() << endl;   // Prints 8
    cout << B.size() << endl;       // Prints 36
    cout << B << endl;              // Prints the reconstructed matrix along with the sparse matrix stuff


    SparseLU<SparseMatrix<double>, COLAMDOrdering<int> > solver;
    solver.isSymmetric(true);
    solver.compute(A);  // Works...
    /*
    ...
    */

    return 0;
}

Also, nonZeros is a member of Matrix, not SparseLU.