避免在使用Eigen分解稀疏矩阵上进行动态内存分配 [英] Avoiding dynamic memory allocation on factorizing sparse matrix with Eigen

问题描述

在我的应用程序中，除类构造函数外，我需要避免动态内存分配(类似malloc).我有一个稀疏的半定矩阵M，其元素在程序执行过程中会发生变化，但它具有固定的稀疏性模式.

In my application I need to avoid dynamic memory allocation (malloc like) except in the class constructors. I have a sparse semidefinite matrix M whose elements change during the program execution but it mantains a fixed sparsity pattern.

为了尽可能快地求解许多线性系统M * x = b，其想法是在类构造函数中使用就地分解，如

In order to solve many linear systems M * x = b as fast as possible, the idea is to use inplace decomposition in my class constructor as described in Inplace matrix decompositions, then call factorize method whenever M changes:

struct MyClass {
private:
    SparseMatrix<double> As_;
    SimplicialLDLT<Ref<SparseMatrix<double>>> solver_;
public:
    /** Constructor */
    MyClass( const SparseMatrix<double> &As ) 
        : As_( As )
        , solver_( As_ ) // Inplace decomposition
    {}

    void assign( const SparseMatrix<double> &As_new ) {
        // Here As_new has the same sparsity pattern of As_
        solver_.factorize( As_new );
    }

    void solve( const VectorXd &b, VectorXd &x )
    {
        x = solver_.solve( b );
    }
}

但是 factorize 方法仍会创建一个大小与 As _ 相同的临时文件，因此请使用动态内存分配.

However factorize method still creates one temporary with the same size of As_, so using dynamic memory allocation.

有可能以某种方式避免它吗?如果Eigen API不允许此功能，则一个想法是创建SimplicialLDLT的派生类，以便仅在将在类构造函数中调用的 analyzePattern 方法中执行动态内存分配.欢迎提出建议...

Is it possible to avoid it in some way? If the Eigen API does not allow this feature, one idea is to create a derived class of SimplicialLDLT so that dynamic memory allocation is performed only in analyzePattern method that will be called in class constructor. Suggestions are welcome...

推荐答案

最后，我找到了使用CSparse库获取H = P * A * P'的解决方法:

At the end I found a workaround using CSparse library to get H = P * A * P':

class SparseLDLTLinearSolver {
private:
    /** Ordering algorithm */
    AMDOrdering<int> ordering_;
    /** Ordering P matrix */
    PermutationMatrix<Dynamic, Dynamic, int> P_;
    /** Inverse of P matrix */
    PermutationMatrix<Dynamic, Dynamic, int> P_inv_;
    /** Permuted matrix H = P * A * P' */
    SparseMatrix<double> H_;
    /** H matrix CSparse structure */
    cs H_cs_;
    /** Support vector for solve */
    VectorXd y_;
    /** Support permutation vector */
    VectorXi w_;
    /** LDLT sparse linear solver without ordering */
    SimplicialLDLT<SparseMatrix<double>, Upper, NaturalOrdering<int>> solver_;
public:
    int SparseLDLTLinearSolver( const SparseMatrix<double> &A )
        : P_( A.rows() )
        , P_inv_( A.rows() )
        , H_( A.rows(), A.rows() )
        , y_( A.rows() )
        , w_( A.rows() )
    {
        assert( ( A.rows() == A.cols() ) && "Invalid matrix" );
        ordering_( A.selfadjointView<Upper>(), P_inv_ );
        P_ = P_inv_.inverse();
        H_ = A.triangularView<Upper>();
        H_.makeCompressed();
        // Fill CSparse structure
        H_cs_.nzmax = H_.nonZeros();
        H_cs_.m = H_.rows();
        H_cs_.n = H_.cols();
        H_cs_.p = H_.outerIndexPtr();
        H_cs_.i = H_.innerIndexPtr();
        H_cs_.x = H_.valuePtr();
        H_cs_.nz = -1;
        const cs_sparse A_cs{
            A.nonZeros(), A.rows(), A.cols(),
            const_cast<int*>( A.outerIndexPtr() ),
            const_cast<int*>( A.innerIndexPtr() ),
            const_cast<double*>( A.valuePtr() ),
            -1 };
        cs_symperm_noalloc( &A_cs, P_.indices().data(), &H_cs_, w_.data() );
        solver_.analyzePattern( H_ );
        // Factorize in order to allocate internal data and avoid it on next factorization
        solver_.factorize( H_ );
        /*.*/
        return -solver_.info();
    }

    int factorize( const Eigen::SparseMatrix<double> &A )
    {
        assert( ( A.rows() == P_.size() ) && ( A.cols() == P_.size() ) &&
            "Invalid matrix size" );
        // Fill CSparse structure
        const cs_sparse A_cs{ 
            A.nonZeros(), A.rows(), A.cols(),
            const_cast<int*>( A.outerIndexPtr() ), 
            const_cast<int*>( A.innerIndexPtr() ), 
            const_cast<double*>( A.valuePtr() ), 
            -1 };
        cs_symperm_noalloc( &A_cs, P_.indices().data(), &H_cs_, w_.data() );
        solver_.factorize( H_ );
        /*.*/
        return -solver_.info();
    }

    void solve( const VectorXd &rhs, VectorXd &x )
    {
        assert( ( rhs.size() == P_.size() ) && ( x.size() == P_.size() ) &&
            "Invalid vector size" );
        // Solve (P * A * P') * y = P * b, then return x = P' * y
        y_ = solver_.solve( P_ * rhs );
        x.noalias() = P_inv_ * y_;
    }
};

cs_symperm_noalloc 是CSparse库的 cs_symperm 函数的次要重构.

cs_symperm_noalloc is a minor refactorization of cs_symperm function of CSparse library.

这似乎行得通，至少与我的特殊问题有关.将来，如果Eigen避免为某些稀疏矩阵操作创建临时变量(进入堆)，将非常有用.

It seems to work, at least with my special problem. In the future it would be very useful if Eigen avoided creating temporaries (into the heap) for some sparse matrix operations.

这篇关于避免在使用Eigen分解稀疏矩阵上进行动态内存分配的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！