ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

c-从本征序列化分解矩阵(SparseLU对象)

2019-10-10 06:08:10  阅读:420  来源: 互联网

标签:c sparse-matrix eigen


我试图解决Ax = b,其中矩阵A可以大到接近1M x 1M的大小,稀疏且对称,但可能没有明确定义.

问题在于,使用本征中的sparseLU object来计算分解可能会花费很长时间,并且可以存储一个sparseLU矩阵而不是原始矩阵,这样,只要我们使用相同的矩阵A执行类似的运算,我们就可以无需重新计算

在stackoverflow上进行的快速搜索和google返回了thisthisthis的稀疏矩阵,用于本征矩阵的序列化.但是,我不确定是否可以将相同的代码应用于sparseLU对象.

也许我应该改一下我的问题:

如何将分解后的矩阵存储到文件中?

当前的方法都集中于存储原始矩阵,但是我想存储分解后的矩阵.有什么办法吗?谢谢.

解决方法:

下面的示例应帮助您实现自己的序列化.

编辑更改示例以回答改写的问题.

#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;
}

另外,nonZeros是Matrix的成员,而不是SparseLU.

标签:c,sparse-matrix,eigen
来源: https://codeday.me/bug/20191010/1884606.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有