1

我是Eigen的新手，我正在处理稀疏 LU 问题。我发现如果我创建一个向量 b(n)，Eigen可以计算 Ax=b 方程的 x(n)。

问题：

  1. 原矩阵A的因式分解结果L&U如何显示？

  2. 如何在Eigen中插入非零值？现在我只是用一些小的稀疏矩阵进行测试，所以我一个一个地插入非零，但是如果我有一个大规模的矩阵，我怎么能在我的程序中输入矩阵呢？

4

2 回答 2

1

我意识到这个问题是很久以前提出的。显然，参考Eigen 文档：

矩阵 L 的表达式，内部存储为超节点 此表达式可用的唯一操作是三角求解

所以没有办法将它实际转换为实际的稀疏矩阵来显示它。Eigen::FullPivLU执行密集分解，在这里对我们没有用处。在大型稀疏矩阵上使用它，我们会在尝试将其转换为密集矩阵时很快耗尽内存，并且计算分解所需的时间会增加几个数量级。

另一种解决方案是使用Suite Sparse中的 CSparse 库：

extern "C" { // we are in C++ now, since you are using Eigen
#include <csparse/cs.h>
}

const cs *p_matrix = ...; // perhaps possible to use Eigen::internal::viewAsCholmod()

css *p_symbolic_decomposition;
csn *p_factor;

p_symbolic_decomposition = cs_sqr(2, p_matrix, 0); // 1 = ordering A + AT, 2 = ATA
p_factor = cs_lu(p_matrix, m_p_symbolic_decomposition, 1.0); // tol = 1.0 for ATA ordering, or use A + AT with a small tol if the matrix has amostly symmetric nonzero pattern and large enough entries on its diagonal
// calculate ordering, symbolic decomposition and numerical decomposition

cs *L = p_factor->L, *U = p_factor->U;
// there they are (perhaps can use Eigen::internal::viewAsEigen())

cs_sfree(p_symbolic_decomposition); cs_nfree(p_factor);
// clean up (deletes the L and U matrices)

请注意，尽管它不像某些 Eigen 函数那样使用显式矢量化，但它仍然相当快。CSparse 也非常紧凑，它只有一个标头和大约 30 个.c文件，没有外部依赖项。易于合并到任何 C++ 项目中。没有必要实际包含所有 Suite Sparse。

于 2015-09-09T03:57:47.650 回答
-1

如果您使用Eigen::FullPivLU::matrixLU()原始矩阵，您将收到 LU 分解矩阵。要分别显示 L 和 U，您可以使用方法triangularView<模式>。在 Eigen wiki 你可以找到它的好例子。将非零值插入矩阵取决于您不想放入的数字。Eigen 具有方便的语法，因此您可以轻松地在循环中插入值：

for(int i=0;i<size;i++)
  {
  for(int j=size;j>someNumber;j--)
    {
    matrix(i,j)=yourClass.getNextNumber();
    }
  }
于 2014-08-22T20:08:35.900 回答