程序中的 LU 分解 (C++)答案

【问题标题】：LU factorization in a program (C++)程序中的 LU 分解 (C++)
【发布时间】：2013-05-24 10:31:21
【问题描述】：

我需要在我的 C++ 程序中执行 LU 分解（分解，有什么区别？）的函数或类。我正在用作编译器 Dev-c 最新版本（在 Windows 上）。我想出了如何安装犰狳和 lapack，但它似乎真的很复杂并且有一些问题http://icl.cs.utk.edu/lapack-for-windows/lapack/（如页面底部所说）。所以我想要一个运行良好的库（在 C++ 中）并且安装起来并不复杂。我在示例中找到了一些关于 Eigen 的东西，怎么样？还有其他建议吗？

谢谢

附：矩阵在对角线上和附近是密集的，在其余部分是稀疏的，并且靠近角（N-E）的部分是空的。

【问题讨论】：

你的矩阵是密集的还是稀疏的？你用什么格式存储它？
矩阵足够稀疏，我还没有对格式进行推理，但这不是问题。

标签： c++ libraries lapack lapack++ matrix-factorization

【解决方案1】：

我为你写了一些 C 代码。也许它可以提供帮助。这是一个将密集矩阵分解为 L, U 的过程，其中 L*U=A，L - 下三角，U - 上三角，L[i,i]=U[i,i]（对角线元素相等） .这种分解也称为 LU(sq)。

#include <math.h>
// A, L, U each allocates at least N*N doubles
// A contains elements of given matrix, written row by row
void decompose(unsigned N, double *A, double *L, double *U) {
#define _(M,i,j) M[N*i + j]
#define _A(i,j) _(A,i,j)
#define _L(i,j) _(L,i,j)
#define _U(i,j) _(U,i,j)
    _L(0,0) = sqrt(_A(0,0));
    _U(0,0) = sqrt(_A(0,0));
    for ( int i = 0; i < N; ++i ) {
        _L(i,0) = _A(i,0) / _A(0,0);
        _U(0,i) = _A(0,i) / _A(0,0);
        _L(0,i) = _U(i,0) = 0.0;

        double s = 0.0;
        for ( int k = 0; k < i; ++k ) {
            s += _L(i,k) * _U(k,i);
        }
        _L(i,i) = _U(i,i) = sqrt(_A(i,i) - s);

        for ( int j = 1; j < i; ++j ) {
            double s = 0.0;
            for ( int k = 0; k < j; ++k ) {
                s += _L(i,k) * _U(k,j);
            }
            _L(i,j) = (_A(i,j) - s) / _L(i,i);
            _L(j,i) = 0.0;

            double s = 0.0;
            for ( int k = 0; k < i; ++k ) {
                s += _L(i,k) * _U(k,j);
            }
            _U(j,i) = (_A(i,j) - s) / _U(i,i);
            _U(i,j) = 0.0;
        }
    }
}

很遗憾，我现在没有时间检查代码中的错误。我有一些我有信心的公式（几年前使用过）：

我的代码基于这些公式。当然，如果你想使用一些特定的适合你的数据稀疏格式，你需要改变程序，而不是公式。

【讨论】：