线性方程组（A是上三角矩阵时）的C++求解

将矩阵$A$通过一系列的线性变换转换为三角矩阵（上三角：upper triangular，下三角：lower triangular），然后便可轻松的求解联立线性方程组。这也是诸如LU分解存在的理由。

---

0、 思路

这里，我们将会看到所谓思路（人脑的思维方式）和计算机的执行方式，并不完全一致。

我们对$A$为上三角矩阵时，线性方程组的解法是这样的，从最后一行开始计算，因为最后一行非零元个数为1，便可轻易地计算得到解$x$的最后一个值（$x[n]$）。利用$x[n]$的值，然后从此向前，逐次计算向量$x$的每一位元素值。

我们以如下的方程组为例：

$$\[ \begin{bmatrix} 1 & 2 & 3 & 4\\ 0 & 5 & 6 & 7 \\ 0 & 0 & 8 & 9 \\ 0 & 0 & 0 & 10\\ \end{bmatrix} \begin{bmatrix} 1\\ 1\\ 1\\ 1\\ \end{bmatrix} = \begin{bmatrix} 10\\ 18\\ 17\\ 10\\ \end{bmatrix} \]$$

其对应的计算流程应当是：
- x[4] = b[4]/A[4][4]
- x[3] = (b[3] - A[3][4]*x[4])/A[3][3]
- x[2] = (b[2] - A[2][4]*x[4] - A[2][3]*x[3])/A[2][2]
- x[1] = (b[1] - A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]) / A[1][1]

一、 matlab解法

clear all; close all; clc
A = [1, 2, 3, 4; 0, 5, 6, 7; 0, 0, 8, 9; 0, 0, 0, 10];
b = [10, 18, 17, 10]';
x = zeros(size(A, 2), 1);
for i = size(A, 1):-1:1,
    x(i) =b(i)/A(i, i);
    j = 1:(i-1);
    b(j) = b(j) - A(j, i)*x(i);
end
x

二、 C++实现

2.1 数据结构的定义：

typedef vector<int> Vec;
typedef vector<Vec > Mat;

这里需要注意Mat矩阵的初始化不同于Vec向量的初始化，Vec可以直接动态的增加元素，Mat却需要首先添加Vec作为其每一行，然后再在Vec中增加元素。

对向量Vec和矩阵Mat的初始化分别为：

int arr[] = { 10, 18, 17, 10 };
Vec b(arr, arr + 4);

Mat A;
for (int i = 0; i < 4; ++i)
    A.push_back(Vec(4));
int count = 1;
for (int i = 0; i < 4; ++i)
    for (int j = i; j < 4; ++j)
        A[i][j] = count++;

2.2 函数接口设计

数学形式：

$$Ax=b$$

为了模拟数学的感觉，变量名我们沿用数学的记号，在不引起混淆的前提下；

void solve(const Mat& A, const Vec& b, Vec& x);

2.3 implementation

void solve(const Mat& A, const Vec& b, Vec& x)
{
    size_t row = A.size();
    Vec b2 = b;                 // 为了保持输入的原始数据不被修改
    for (int i = row - 1; i >= 0; --i)
    {
        x[i] = b2[i] / A[i][i];
        for (int j = 0; j < i; ++j)
            b2[j] -= A[j][i] * x[i];
    }
}

这里需要注意两个问题：

unsigned类型变量为负时：

for (size_t i = row-1; i >= 0; --i)     
// 或者unsigned i;
// typedef _W64 unsigned int   size_t;

当i减少到0，再执行减一的动作时，不会变成-1，而是4294967295，也即($2^32-1$)。此时循环不会退出，但如果试图索引x[4294967295]，则会报错。

返回局部类对象的引用：

切勿返回局部类对象的引用（如下代码所示），因为函数调用完毕，编译器执行清栈操作，会调用局部对象的析构函数。对象如果不存在的话，何谈其引用呢。

Vec& solve(const Mat& A, const Vec& b)
{
    size_t row = A.size();
    Vec x(row);
    ...
    return x;
}

Vec x = solve(A, b);    // x.size() == 0

对问题2的改造：

如2.3节所示：

void solve(const Mat& A, const Vec& b, Vec& x);     
// 在客户端完成对x的初始化

或者：

Vec solve(const Mat& A, const Vec& b);              
// 返回的是类对象，而不是类对象的引用，这样在函数退出时，会调用Vec类的拷贝构造函数，创建并赋值一个临时（未名）对象

三、 实现方法之二

我们再来观察当 A <script type="math/tex" id="MathJax-Element-119">A</script>是四阶矩阵的计算流程：

• x[4] = b[4]/A[4][4]
• x[3] = (b[3] - A[3][4]*x[4])/A[3][3]
• x[2] = (b[2] - A[2][4]*x[4] - A[2][3]*x[3])/A[2][2]
• x[1] = (b[1] - A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]) / A[1][1]

A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]这种A的行和x的乘积再加和的形式对应着线性代数的一个经典的概念就叫内积。我们便由此可以利用A的行和x向量已解部分的内积来进行相关的计算。

3.1 matlab

clear all; close all; clc
A = [1, 2, 3, 4; 0, 5, 6, 7; 0, 0, 8, 9; 0, 0, 0, 10];
b = [10, 18, 17, 10]';
x = zeros(size(A, 2), 1);
r = size(A, 1);
for i = r:-1:1,
    j = i+1:r;
    x(i) = (b(i) - A(i, j)*x(j)) / A(i, i);
end
x

因为未解部分已被初始化为0，可进一步将程序简化为：

for i = r:-1:1,
    x(i) = (b(i) - A(i, :)*x) / A(i, i);
end

3.2 C++

void solve2(const Mat& A, const Vec& b, Vec& x)
{
    size_t r = A.size();
    for (int i = r - 1; i >= 0; --i)
    {
        int tmp = b[i];
        for (int j = 0; j < r; ++j)
            tmp -= A[i][j] * x[j];
        x[i] = tmp / A[i][i];
    }
}