在我的应用程序中,除了在类构造函数中,我需要避免动态内存分配(类似malloc).我有一个稀疏的半定矩阵M,它的元素在程序执行过程中会发生变化,但它保持着固定的稀疏性模式。
为了尽可能快地解决许多线性系统M * x = b,我的想法是在我的类构造函数中使用inplace分解,如在 置换矩阵分解,然后叫 因式化 每当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 );
}
}
但是 因式化 方法仍然会创建一个临时的、具有相同大小的 As_,所以使用动态内存分配。
有没有可能用某种方式避免呢?如果Eigen API不允许这个功能,一个想法是创建一个SimplicialLDLT的派生类,使动态内存分配只在 分析模式 方法,将在类构造函数中调用。欢迎提出建议...
最后,我找到了一个使用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_noalloc 的一个小重构。cs_symperm CSparse库的函数。
至少在我的特殊问题上,它似乎可以工作。将来,如果Eigen能避免为一些稀疏矩阵操作创建临时变量(到堆中),那将是非常有用的。