Как избежать динамического распределения памяти c на разложении факторизованной матрицы с помощью Eigen

Question

В моем приложении мне нужно избегать динамического выделения памяти c (подобно mallo c), за исключением конструкторов классов. У меня есть разреженная полуопределенная матрица M, элементы которой изменяются во время выполнения программы, но она имеет фиксированный шаблон разреженности.

Чтобы решить многие линейные системы M * x = b как можно быстрее, идея состоит в том, чтобы использовать декомпозиция на месте в конструкторе моего класса, как описано в Декомпозиции матрицы на месте , затем вызывать метод factorize при каждом изменении M:

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_ , поэтому используется динамическое выделение памяти c.

Можно ли каким-то образом избежать этого? Если Eigen API не разрешает эту функцию, одной из идей является создание производного класса SimplicialLDLT, чтобы динамическое выделение памяти c выполнялось только в analyPattern методе, который будет вызываться в конструкторе классов. Предложения приветствуются ...

Answer 1

В конце я нашел обходной путь, используя библиотеку CSparse для получения 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_noallo c - это небольшая рефакторизация cs_symperm функция библиотеки CSparse.

Кажется, это работает, по крайней мере, с моей особой проблемой. В будущем было бы очень полезно, если бы Эйген избежал создания временных (в кучу) для некоторых разреженных операций с матрицами.