数值分析LU分解法C++

### 使用C++实现LU分解

在数值分析领域，LU分解是一种重要的工具，用于高效求解线性方程组。对于C++而言，可以利用标准模板库(STL)以及Eigen等第三方库来简化矩阵操作。

#### 不借助外部库的简单实现

下面展示了一个不依赖任何额外库的基础版本：

#include <iostream>
#include <vector>

using namespace std;

void luDecompose(vector<vector<double>>& A, vector<int>& index) {
    int n = A.size();
    double maxA;
    int imax;
    
    for (int i = 0; i < n; ++i) {
        maxA = 0.0;
        imax = i;
        
        for (int k = i; k < n; ++k)
            if (abs(A[k][i]) > maxA) {
                maxA = abs(A[k][i]);
                imax = k;
            }

        if (!maxA) { cout << "Error: Singular matrix"; exit(1); }
        
        if (imax != i) {
            swap(index[i], index[imax]);

            for (int j = 0; j < n; ++j)
                swap(A[i][j], A[imax][j]);
            
            // Keep track of row swaps.
            for (int k = 0; k < n; ++k)
                if (index[k] == i) index[k] = imax;
                else if (index[k] == imax) index[k] = i;
        }
        
        index[i] = imax;
        
        for (int j = i + 1; j < n; ++j) {
            A[j][i] /= A[i][i];
            for (int k = i + 1; k < n; ++k)
                A[j][k] -= A[j][i]*A[i][k];
        }
    }
}

// Solve Lz=b and Ux=z using forward substitution followed by backward substitution.

void solve(const vector<vector<double>>& LU, const vector<int>& idx, const vector<double>& b, vector<double>& x) {
    size_t n = LU.size();

    for (size_t i = 0; i < n; ++i){
        x[i]=b[idx[i]];
        for(size_t j=0;j<i;++j)x[i]-=LU[i][j]*x[j];
    } 

    for(int i=n-1;i>=0;--i){
        for(int j=i+1;j<n;++j)x[i]-=LU[i][j]*x[j];
        x[i]/=LU[i][i];
    }
}

此代码实现了部分选主元策略下的Doolittle形式的LU分解算法[^1]。通过上述函数`luDecompose()`完成输入矩阵\( \textbf{A} \)，将其转换成下三角矩阵\(\textbf{L}\)(含单位对角线元素)与上三角矩阵\(\textbf{U}\) 的乘积；随后调用`solve()`执行前向替代和回带过程得到最终解向量\[ \vec{x} \].

为了提高效率并减少编程工作量，在实际项目开发过程中推荐采用成熟的开源数学库如[Eigen](https://eigen.tuxfamily.org/dox/), 它们已经内置了优化良好的LU分解功能以及其他丰富的线性代数运算接口。