从本征序列化分解矩阵（SparseLU对象） [英] Serializing Decomposed Matrix from eigen (SparseLU object)

问题描述

我想解决 Ax = b ，其中矩阵A可以大到接近 1M x 1M 大小，是稀疏和对称的，但可能没有正确定义。

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.

问题是，可能需要很长时间才能使用 sparseLU对象，并且将会想到一个存储sparseLU矩阵而不是原始矩阵，使得每当我们使用相同的矩阵A执行类似的操作时，我们可以不需要重新计算

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

快速搜索stackoverflow和google返回这，这和这稀疏矩阵用于特征矩阵的序列化。但是，我不确定是否相同的代码可以应用于sparseLU对象。

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:

如何将分解的矩阵存储到文件？

当前方法都集中在存储原始矩阵，但我想存储分解矩阵。有什么办法吗？

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.

推荐答案

以下示例将帮助您实现自己的序列化。

The following example should help you implement your own serialization.

EDIT 已更改示例以回答重新编写的问题。

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*)&inSz , sizeof(IND));
        readFile.read((char*)&outSz, 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;
}

此外， nonZeros 是 Matrix 的成员，而不是 SparseLU 的成员。

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

这篇关于从本征序列化分解矩阵（SparseLU对象）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！