在第二节中，我们学习了行阶梯矩阵与行化简矩阵，并且使用程序实现了相应的算法，但是也挖了一个坑：要如何使用行阶梯型矩阵来判断线性方程组解的情况，在本节将给出答案并使用C++实现相应的判断解情况的函数。

线性代数部分

1.主元

  在第二节中，我们了解到一个矩阵可以通过行初等变化变化出无数个行阶梯矩阵(零矩阵只能变化出一个)，不过矩阵能且只能通过行初等变化变化出唯一一个行化简矩阵，但是注意到无论是行阶梯矩阵还是行化简矩阵，每一行先导元素在矩阵中所处的位置都是固定的，这一点可由行化简矩阵的唯一性推出，首先每一个矩阵通过初等行变化变化出的行化简矩阵都是唯一的，那么在行化简矩阵中每一行先导元素的位置就必然是固定的，而行化简矩阵比行阶梯型矩阵多出来的两条性质：1.先导元素必须唯一；2.先导元素所在列的除先导元素以外的元素必须为零，这两条性质都不改变先导元素本身的位置，由此可以推知，在行阶梯型矩阵中先导元素的位置也必然是固定的，我们称先导元素所在的固定的位置为主元位置，有主元位置的列称为主元列比如下方的行阶梯型矩阵：
$$\left[\begin{array}{cccc}2& 3& 1& 4\\ 0& 1& 0& 4\\ 0& 0& 2& 8\\ 0& 0& 0& 1\end{array}\right]$$
  在该矩阵中斜对角线上的所有位置都是主元位置，由于每一列都存在主元位置，因此每一列都是主元列，又比如下方的矩阵：
$$\left[\begin{array}{ccccc}\ast & x& x& x& x\\ 0& 0& \ast & x& x\\ 0& 0& 0& \ast & x\\ 0& 0& 0& 0& 0\end{array}\right]$$
  在该矩阵中，最后一行没有先导元素，因此不存在主元位置，在第一，第二，第三行中的$\ast$位置就是主元位置，注意到第二列与第五列不存在主元位置，因此它们就不是主元列，上方是笔者的理解，在教材的标准定义中：主元位置是一个矩阵对应其行化简矩阵中元素1所在的位置，主元列就是含有主元位置的列，这里有一个原矩阵与行化简矩阵的对应关系，也就是说在标准定义中，主元位置更加强调的是在原矩阵中的位置，比如下方的矩阵与其行化简矩阵：
$$\left[\begin{array}{ccc}x& y& y\\ y& x& y\\ y& y& x\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$
  在上方的原矩阵中，$x$的位置就是主元位置，在左侧原矩阵中的每一列都含有主元位置，因此每一列都是主元列，要注意的是，只要最后化简出来的行化简矩阵是右侧那样的，那么$x$的位置就一定是主元位置，与原矩阵中$x$的值没有任何关系，并且$x$存在的列一定是主元列。但是我们可以推知，在行阶梯矩阵与行化简矩阵中，主元位置上的元素一定是非零的，在此时，我们称主元位置上的非零元素为主元，可以认为主元其实就是行阶梯型矩阵与行化简矩阵中每一行的先导元素，要说明的是，主元是一个非常重要的基础概念，后面许多复杂的矩阵性质都与主元有关。

2.线性方程组的解

  在第二节中，我们总结出在第一节中解线性方程组的步骤本质上就是先把增广矩阵变化为行阶梯型矩阵，然后再把行阶梯型矩阵化简为行化简矩阵，我们称这种解线性方程组的方法为行化简算法，称把增广矩阵变化为行阶梯型矩阵的步骤为解线性方程组的向前步骤，把行阶梯型矩阵化简为行化简矩阵的步骤称为解线性方程组的向后步骤，在实际设计算法时，我们发现向前步骤要比向后步骤更加复杂。使用行化简算法，我们可以解出任意有解的线性方程组，在本小节讨论的重点就是解的情况以及如何使用行阶梯矩阵判断解的情况，我们先来讨论解的情况。
  总体上，一个线性方程组解的情况可以分为有解与无解，也就是线性方程组相容与不相容，我们先来讨论线性方程组相容的情况，下面先来看一个具体的例子：
$$\{\begin{array}{l}{x}_1+x_2+x_3=3\\ 2x_1+2x_2+2x_3=6\\ 2x_1+2x_2=2\end{array}$$
增广矩阵对应的行化简矩阵如下：
$$\left[\begin{array}{cccc}1& 1& 0& 1\\ 0& 0& 1& 2\\ 0& 0& 0& 0\end{array}\right]$$
回代到线性方程组中得：
$$\{\begin{array}{l}{x}_1+x_2=1\\ x_3=2\end{array}$$
  现在我们得知，该矩阵有两个主元，分别是第一行第一个元素和第二行第三个元素，那么对应的就存在两个主元列，我们称线性方程组中对应于主元列的未知数为基本变量，也就是位于主元位置的未知数就是基本变量，其余位置的未知数都称为自由变量，比如在上方的例子中，$x_1，x_3$就是基本变量，$x_2$就是自由变量，在解线性方程组时，我们会使用自由变量表示基本变量，由此得出的就是线性方程组的解，也就是说上方例子实际上的解是：
$$\{\begin{array}{l}{x}_1=1-x_2\\ x_3=2\end{array}$$
  在线性方程组有解时，我们得出的最终解就是上方这种基本变量位于等号左侧，自由变量位于等号右侧的形式，当然实际上我们也可以使用基本变量表示自由变量，不过这不符合线性代数表示的标准，注意到$x_3$的值是固定的，不能用自由变量表示，对于这种情况，不用自由变量就可以了，也就是说在只存在唯一解时，我们得出的解实际上也是上方的这种表示形式，只不过此时只存在基本变量，不存在自由变量，由此可得，解线性方程组的最终目的就是把线性方程组变化为每一个基本变量都独立位于等号左侧的形式。
  对于自由变量，还可以从几何的视角来看，我们先从最简单的几何直觉中查找规律：如果线性方程组有解且没有自由变量，那么就是唯一解，在几何上就是一个固定的点，对应的维度是零维，如果线性方程组的有解且存在一个自由变量，那么解在几何上就是一条直线，对应的维度就是一维，从几何直觉上来看，自由变量的个数肯定与解的维度有关，并且解的维度应该是等于自由变量个数的：零个自由变量的解对应零维点，一个自由变量的解对应一维线，二个自由变量的解对应二维面…并且在第之前我们推测过，线性方程组本身的维度是由线性方程组中维度最高的那个线性方程决定的，也就是未知数最多的那个，显然，解的维度必然是小于原线性方程组的维度的，但是解本身所在的维度必然是原线性方程组的维度，也就是说解集的维度是原线性方程组维度的子集，比如在平面上的两条直线相交，此时这两条直线位于平面上，那么其对应的线性方程组就是二维的，相交后得出一个点，这个点本身是零维的，但是其所处的维度是二维是，也就是二维中的一个零维点，二维中显然包含无数个点，因此就称这一个二维中的点是二维面的一个子集。
  在教材的标准定义中，我们将使用自由变量表示基本变量的线性方程组的形式称为线性方程组的通解，也就是在等号左侧只存在基本变量，等号右侧存在自由变量与已知量的形式，该形式显示表示出了解的所有情况，解线性方程组就是把线性方程组化为该形式或确定线性方程组无解。

2.1解的存在性与唯一性问题

  在本小节正式说明要如何使用行阶梯型矩阵判断线性方程组是否相容，如果相容的话，解是否唯一，先来看下面的三个增广矩阵对应的行化简矩阵：
$$\left[\begin{array}{cccc}1& 0& 0& 5\\ 0& 1& 0& 3\\ 0& 0& 1& 8\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 3\\ 0& 1& 0& 4\\ 0& 0& 0& 5\end{array}\right]\left[\begin{array}{cccc}1& 0& 1& 8\\ 0& 1& 2& 5\\ 0& 0& 0& 0\end{array}\right]$$
对应的线性方程组为：
$$\{\begin{array}{l}{x}_1=5\\ x_2=3\\ x_3=8\end{array}\{\begin{array}{l}{x}_1=3\\ x_2=4\\ 0=8\end{array}\{\begin{array}{l}{x}_1=8-x_3\\ x_2=5-2x_3\\ 0=0\end{array}$$
  上方的情况分别对应唯一解，无解与无数解，重点观察无解情况下的增广矩阵，发现增广矩阵的最右侧一列有主元位置，是主元列，但是在线性方程组中，最右侧一列显然没有未知数，因此就得出了$0=8$，由此推测存在性定理与唯一性定理：

定理
1.如果增广矩阵的最右侧一列不是主元列，也就是没有这种形式的行：$(0,\mathrm{0..},0,b)$那么该增广矩阵对应的线性方程组相容，否则不相容
2.在线性方程组相容的基础上，如果不存在自由变量，那么线性方程组有唯一解

  由于当矩阵变化为行阶梯型矩阵时，就已经能够判断出主元位置和主元列了，只需要把增广矩阵变化为行阶梯型矩阵，然后再判断最右侧一列是否是主元列就能够判断出线性方程组是否相容了，在相容的基础上，如果线性方程组其余的列都是主元列，那么就不存在自由变量，此时线性方程组存在唯一解。

笔者的奇思妙想(不保证正确性):
1：未知数个数与线性方程组维度的关系
  在前方的章节中，我们推断一个线性方程组本身的维度是由线性方程组中维度最高的线性方程所决定的，但是具体是怎么决定的？该推断是正确的吗？考虑下面两个线性方程组：
$$（1）\{\begin{array}{l}{x}_1+x_2=5\\ 3x_1+x_2=8\end{array}（2）\{\begin{array}{l}{x}_1+x_2=5\end{array}$$
这两个线性方程组中，未知数最多的线性方程都有两个未知数，那么这两个线性方程组所处的维度就是相同的吗？先来看看第二个线性方程组，该线性方程组中只有一个线性方程，该线性方程存在两个未知数，很明显其几何形式是一条线，线性方程组所处的维度就是一维的，再来看看第一个线性方程组，其中有两个线性方程，都含有两个未知数，那么几何形式就是两条线，在一维中显然无法表示出两条线，因此该线性方程组所处的维度是二维的，那么为什么不能是三维的，在三维体中显然是能够表示出二条线的，说到底，线性方程组的维度是什么，在教材中并没有该术语的定义，这是笔者为了便于理解而造出来的一个名词，在最初我们根据直觉推测，线性方程组的维度就是线性方程组中维度最高的线性方程所处的维度，但是此时显然出现问题了，即线性方程的个数似乎也会影响到线性方程组的维度，并且高维能够表示低维的信息，在二维中能够表示出两条线，在三维中同样也可以，那么线性方程组的维度存在明确的上限吗？
  首先要搞明白第二个问题，否则线性方程组的维度就是无意义的，从代数的角度进行分析，考虑一个实际问题：使用两个未知数能够表示出三维空间中的一条线吗？我们可以认为，任意维度的基本空间单位都是点，所谓表示出任意维度的一个点就是定位到该点，比如在二维平面中，配合二维数轴我们可以使用$(x,y)$唯一的定位到二维空间中的一个点，而无数个存在关系的点就可以表示出表示出一条线，这里的关系就是指未知数$x$与$y$的关系，也就是所谓的函数，从几何的角度含，函数本质上是一种几何约束，比如$(x,y)$本身是能够表示出二维平面上任意一个点的，那么就相当于表示出了这个平面，这意味着$(x,y)$就本身表示了一个平面，但是加上函数约束后变为$y=kx+b$，此时表示的显然就是一条线了，可以推出一个结论：函数关系会让未知数表示的维度降低一维，又比如$(x,y,z)$表示的显然就是三维空间中的一个点，在没有函数约束关系前，$(x,y,z)$能够表示出三维空间中任意一个点，就相当于表示出了整个三维空间，但是在加入了函数关系后变为：$ax+by+c=z$，或者换个形式：$a{x}_1+b{x}_2+c{x}_3=0$显然表示的确实是一个平面，那么再提出一个问题：$(x,y)$能够表示出三维的一个点吗？显然是不能的，三维的点显然要使用三个未知数表示，也就是说零维不需要未知数表示，一维要使用一个未知数表示，二维要使用两个未知数表示…n维要使用n个未知数表示，既然如此，那么线性方程组的维度必然就是有上限的，并且一定和未知数的个数有关。在得出结论前，我们先将线性方程组的维度这个概念修正为：线性方程组所处的维度，此时概念就明确多了，先给这个名词一个具体的定义：线性方程组所处的维度就是能够表示出线性方程组中每一个线性方程几何表示的维度，然后再来看看上方的线性方程组：
$$（1）\{\begin{array}{l}{x}_1+x_2=5\\ 3x_1+x_2=8\end{array}（2）\{\begin{array}{l}{x}_1+x_2=5\end{array}$$
  首先第二个线性方程组真的位于一维吗，还是说它是二维中的一条线，这是两个完全不同的概念，从未知数的个数来看，存在两个未知数就意味着其所处的维度必然是二维是，只不过受到了函数关系的约束，这两个未知数从能够表示完整的二维退化成了表示二维的一条线，因此第二个线性方程组所处的维度是二维的，从几何的角度来看，在一维是不可能完整的表示出一条线的，一个未知数$x$能够表示出一维中任意一个点，但是在任意时刻该未知数只能表示出点，单个未知数绝对无法在同一时刻表示出完整的整个一维信息，要在任意时刻完整的表示出一个完整的一维信息，必须使用两个未知数，一但使用了二个未知数，那么维度就上升到了二维，因为$(x,y)$表示的一定是二维的点，因此该线性方程组是位于二维的，其中的线性方程表示的是二维中的一条线，至此两个问题都解决了，我们可以得出线性方程组所处的维度是由线性方程组中维度最高的线性方程所决定的，并且与该线性方程的维度相同，具体的维度值就等于该线性方程未知数的个数，比如线性方程组：
$$\{\begin{array}{l}{x}_1+x_2+x_3=5\\ x_3=8\end{array}$$
该线性方程组所处的维度就是三维，整个线性方程组共同表示出了三维中的一条直线。
  2.所处维度与表示维度，一个点可以处于零维，一维，二维，三维…，但是该点表示的维度只有零维，由此笔者将维度划分为所处维度与表示维度，显然在线性方程与线性方程组中表示维度绝对小于所处维度，这是由函数的约束关系导致的，比如线性方程$x_1+x_2+x_3=3$，其表示一条直线，表示维度是二维，但是其处于三维，所处维度是三维，完整的表述就是这是一条处于三维空间中的直线，下面来考虑线性方程组的情况：
$$\{\begin{array}{l}{x}_1+x_2+x_3=5\\ x_3=8\end{array}$$
在该线性方程组中，线性方程组所处维度是三维，其表示维度就是最终解的维度，显然解是一条直线，因此表示维度就是一维，对于该线性方程组的完整描述就是：三维中的一条直线，由此得到启发，将线性方程组的维度的概念进行一次扩充，得出：线性方程组的维度分为所处维度与表示维度，其所处维度等于方程组中未知数最多的线性方程拥有的未知数的个数，其表示维度等于解的维度，而解的维度等于自由变量的个数。
3.维度单位与维度信息单位
  对于任意维度，点都是其基本单位，但是在各个维度中，表示一个点所需的信息显然是不同的，在零维中，表示点不需要任何信息，零维本身就是点，在一维中表示点需要一个未知数，在二维中表示点需要两个未知数…在n维中，表示点需要n个未知数，并且可以得出，在各个维度表示点所需的未知数的个数与维度是严格一一对应的，比如两个未知数无论如何也无法表示三维中的一个点，并且在一维中也绝对无法容纳两个未知数的信息，也就是说对于高维度，由于信息不足而无法表示，对于第维度，由于信息量太大而无法表示，这或许就是世界的一种保护机制吧，严格的隔离开了各个维度之间的信息，就像虽然三维是由无数个二维组成的，但是我们无论如何也无法真正的在三维中抽出一个完整的二维，也就是说本维度与本维度的信息是一一对应的，三个未知数就能且只能表示三维的点，无论如何也表示表示其它维度的信息，由此我们可以得出维度的信息单位 = 表示本维度的一个点所需要的未知数的个数。一维的维度信息单位就是一，二维的就是二…n维的就是n。
4.升维与降维
  线性代数中具体的升维与降维是非常复杂的，在这里我们不谈复杂的算法与感念，仅仅提出一个疑问：升维与降维可能吗？从维度信息与维度一一对应的关系来看，这个世界似乎不存在升维与降维的可能性，但是不谈这么死板的底层原理，在思维的世界中有无限可能，从人类认识到升维与降维的那一刻起，就一定有可能实现升维与降维，现在让我们放飞思想，来探讨升维与降维的方案吧，首先我们得假定我们所处在二维平面世界，因为我们无法想象四维，但是我们可以想象出零，一，二，三维，这是我们已知的，遵循从已知推出未知的原则，想象此时我们处在二维平面，往下是一维，往上是三维，我们都可以想象，在二维平面中没有高度信息，高度对于二维生物就像是我们无法感知到四维的那一个独有信息一样，此时身在二维世界的我们发现，就算我们能够表示出整个二维所有的信息，依据无法达到三维，因为单个二维平面绝对不可能出现高度信息，因此我们考虑复制整个二维世界，以已经表示出来的二维世界信息为模板，不断的复制新的平面出来，从第一个世界复制成功开始，我们就能够成功感知到高了，然后不断的复制，最终实现升维，由此得到启发，不考虑能量，只考虑信息，通过Copy或许就能够实现升维，显然最终的答案又回到了计算机，最为人类目前最强的工具，计算机是唯一有可能实现在信息层面Copy世界的，笔者将其称为Copy升维方案，并且在物理学中我们得知，宇宙是在以光速不断扩充的，谈点"浪漫"的，这是否是宇宙在Copy自己呢，这个宇宙本身是不是也想实现升维？我们不得而知，下面来看看降维，我们是能够想象出二维的，并且也能够知道我们所处的三维世界就是由无数个二维世界构成的，因此可以考虑两个方案，方案一是从无数个二维中抽出来一个，方案二是将三维进行压缩，把信息压缩到二维，看到方案二，笔者联想到的就是引力与黑洞，黑洞不就向是一个宇宙的压缩包吗，连光信息都无法逃出黑洞的引力，但是黑洞的压缩是存在上限的，当"内存"满了，就会触发黑洞喷流，被压缩的信息就会被释放出来，由此得到启发，如果能够实现一个人工黑洞，该黑洞可以由人类投入指定量的信息与能量，此后就开始不断的压缩，由于没有不断的投入能量与信息，由此黑洞喷流就不会发生，或许到达某一个压缩极限时，二维世界就诞生了。然后是抽离方案，我们首先得捕获到三维中的一个二维平面，然后是存储，我们得先构造出一个"空维“空间，在该空间中不存在然后信息与能量，然后把捕获到的二维平面投入该空维，由于能量会逸散，因此还得加固我们捕获到的二维平面，这个方案显然要复杂不少，因为没有实际的物理现象供我们参考，笔者这两个降维方案称为压缩降维与抽离降维，当然知名的降维还有质子展开，这些东西听起来像科幻，但是要知道《海底两万里》在几十年前也是科幻，登月那更是无法想象的事情，但是现在的人类早已远远超过了两万里，登的是火星了，笔者一直相信，数理公式无法约束人类的上限，思维所想即是人类所达。

编程部分

  在上两节中，我们使用C++实现了表示矩阵，解线性方程组，判断行阶梯型矩阵与行化简矩阵等关键函数，现在让我们以新学习到的线性代数知识为指导，实现一批新的函数，具体函数如下：

函数
1.给定一个矩阵并指定行，判断该行的主元位置
2.给定一个矩阵并指定列，判断该列是否是主元列
3.给定一个矩阵，判断矩阵解的情况

  首先的第一个函数：给定一个矩阵并指定行，判断该行的主元位置，显然我们得先把矩阵变化为行阶梯型矩阵，然后就能够找到主元位置，接口设计如下：

ssize_t FindPivotPosition(matrix& m, size_t row)

如果指定行存在主元位置，那么就返回主元位置在行中的下标，如果不存在就返回$-1$，函数实现如下：

ssize_t FindPivotPosition(matrix &m, size_t row)
{
    int row_sz = m.size();
    int col_sz = m[0].size();
    if (row > row_sz)
    {
        return -1;
    }

    // 如果是行阶梯型矩阵，那么直接在指定行中找主元
    if (IsREF(m))
    {
        for (int i = 0; i < col_sz; i++)
        {
            if (m[row][i])
            {
                return i;
            }
        }
        // 是全零矩阵，不存在主元位置，返回-1
        return -1;
    }
    // 拷贝出一份新矩阵，不改变原矩阵
    matrix tmp_m = m;

    // 将矩阵化为行阶梯型矩阵
    PivotDiagonal(m);

    for (int i = 0; i < col_sz; i++)
    {
        if (m[row][i])
        {
            return i;
        }
    }
    // 是全零行，不存在主元位置，返回-1
    return -1;
}

显然中间有一部分可以提取出来，便于理解的完整版代码如上，简化版代码如下：

ssize_t FindOneRowPP(matrix& m, size_t row, size_t col_sz)
{
    for (int i = 0; i < col_sz; i++)
    {
        if (m[row][i])
        {
            return i;
        }
    }
    return -1;
}

ssize_t FindPivotPosition(matrix &m, size_t row)
{
    int row_sz = m.size();
    int col_sz = m[0].size();
    if (row > row_sz)
    {
        return -1;
    }
    if (IsREF(m))
    {
        return FindOneRowPP(m, row, col_sz);
    }
    matrix tmp_m = m;
    PivotDiagonal(m);
    return FindOneRowPP(m, row, col_sz);
}

测试用例如下(测试简化版)：
$$(1)\left[\begin{array}{cccc}1& 0& 0& 5\\ 0& 1& 0& 3\\ 0& 0& 1& 8\end{array}\right](2)\left[\begin{array}{ccc}1& 3& 5\\ 2& 1& 3\\ 3& 4& 5\\ 8& 1& 9\end{array}\right](3)\left[\begin{array}{cccc}1& 0& 1& 8\\ 0& 2& 2& 5\\ 0& 0& 0& 0\end{array}\right]$$
先来人眼判断一下，对于矩阵一，传入第零行会返回0，第一行返回1，第二行返回2，对于矩阵二，我们得先把它变化为行阶梯型矩阵，结果如下：
$$\left[\begin{array}{ccc}1& 3& 5\\ 0& 5& 7\\ 0& 0& 3\\ 0& 0& 0\end{array}\right]$$
因此传入第零行会返回0，第一行返回1，第二行返回2，第三行返回-1，最后是矩阵三，显然传入第零行返回0，第一行返回1，第二行返回-1，由于第二个矩阵最具有普遍性，先使用矩阵二进行测试，测试代码如下：

void TestFindPivotPosition()
{
    matrix m = {{1, 3, 5}, {2, 1, 3}, {3, 4, 5}, {8, 1, 9}};
    for(int i = 0; i < m.size(); i++)
    {
        int tmp = FindPivotPosition(m, i);
        if(tmp != -1)
        {
            std::cout << "第" << i << "行" << "主元位置为: " << tmp << std::endl;
        }
        else{
             std::cout << "第" << i << "行" << "不存在主元位置" << std::endl;
        }
    }
}

显然结果错误，原因出在行阶梯型矩阵变化函数上，该函数代码如下：

void PivotDiagonal(matrix &m)
        {
            int row_size = m.size();
            int column_size = m[0].size();

            for (int i = 0; i < row_size; i++)
            {
                // 1.找到对角线元素为0的行
                if (m[i][i] == 0)
                {
                    for (int j = 0; j < column_size; j++)
                    {
                        // 2.找到该行的其它非零元素所在的列数，这些数字就是可以与该行进行对换的行号
                        if (m[i][j])
                        {
                            // 3.在可对换行中找到能够让对角线元素变为非0的行，最后进行对换
                            if (m[j][i])
                            {
                                RowSwapping(m, i, j);
                                break;
                            }
                        }
                    }
                }
            }
        }

  在第一节中，我们假设线性方程组有解并且使用的是不具备普遍性的测试用例，仅仅只实现最基本的功能，现在我们得知，行阶梯型矩阵变化不是使用对换变化就可以完成的，我们必须实现出一个通用的函数，下面对函数进行修改。我们第一版的行阶梯型变化矩阵函数由于没有根据行阶梯矩阵的性质来设计算法，因此能够变化的是只需要通过对换变化就能够得出行阶梯型的矩阵，并且只能处理对角线上的元素，太过于与绝对了，改版都很难改，因此根据行阶梯矩阵的性质直接设计新的具有普适性的算法(从这里也可以看出数学性质的指定对于算法设计有多么重要了)，行阶梯矩阵的性质如下：

性质
1.非零行必须位于零行的上方
2.先导元素所在列下方的元素必须为全零
3.每一行先导元素所在列都必须位于上一行先导元素的后方

并且还要注意：全零矩阵也是行阶梯型矩阵，现在我们得根据上方的性质实现一个具有普适性的行阶梯矩阵的变化算法，我们一步步实现，首先是：非零行必须位于零行的上方，考虑使用双指针，上方指针从第一行开始遍历，下方指针从最后一行开始遍历，上方指针变量到零行时停止，下方指针遍历到非零行时停止，然后交换，交换后继续遍历，直到上下指针遍历到同一行时，此时该行上方的必然全是非零行，下方的必然的零行，那么该行是什么行就不重要了，因此遍历的条件是：上方指针 < 下方指针，等于时不需要遍历，函数实现如下：

void SwapRow(matrix& m, size_t row1, size_t row2)
{
    if(row1 > m.size() || row2 > m.size())
    {
        return;
    }
    std::vector<double> tmp = std::move(m[row1]);
    m[row1] = std::move(m[row2]);
    m[row2] = std::move(tmp);
}

//将零行全部移动到非零行的下方，返回移动后最下方非零行的行号
int SwapZeroRowAndNotZeroRow(matrix& m)
{
    int row_sz = m.size();
    int col_sz = m[0].size();
    int up = 0;
    int down = row_sz - 1;

    //如果矩阵本身就是行阶梯型矩阵，那么找出最下方非零行的行号，然后返回，如果是全零矩阵，那么返回行号0
    if(IsREF(m))
    {
        for(int i = row_sz - 1; i >= 0; i--)
        {
            if(!IsZeroRow(m, i))
            {
                return i;
            }
        }
        return 0;
    }

    int up_flag = 0;
    int down_flag = 0;

    while(up < down)
    {
        //上方指针遍历查找零行，如果找到了就停止并保持值不变，如果没有找到那么就加一遍历下一行
        if(IsZeroRow(m, up))
        {
            up_flag = 1;
        }
        else
        {
            up++;
        }

        //下方指针遍历查找非零行，如果找到了就停止并保持值不变，如果没有找到那么就减一遍历上一行
        if(IsZeroRow(m, down))
        {
            down--;
        }
        else
        {
            down_flag = 1;
        }

        //如果两个都找到了，那么交换两行的位置，交换后让up++,down--,继续往后遍历
        if(up_flag == 1 && down_flag == 1)
        {
            SwapRow(m, up, down);
            up++;
            down--;
            up_flag = 0;
            down_flag = 0;
        }
    }
    return down;
}

其中的IsZeroRow在第二节有详细的实现介绍，测试用例如下：

void TestSwapZeroRowAndNotZeroRow()
{
    matrix m = {{1, 1, 2, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 1}};
    std::cout << "交换前: " << std::endl;
    ShowMatrix(m);

    int row = SwapZeroRowAndNotZeroRow(m);

    std::cout << "交换后: " << std::endl;
    ShowMatrix(m);

    std::cout << "最下方非零行行号: " << row << std::endl;
}

测试结果如下：

结果正确，下面来实现第二个性质：先导元素所在列下方的元素必须为全零，在查找先导元素前，必须做一些准备工作，以下方矩阵为例：
$$\left[\begin{array}{cccc}1& 0& 3& 4\\ 2& 0& 6& 8\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$
  显然上方的矩阵是符合性质一的，但是存在一些非常恶心的情况，比如第二行的元素是第一行的倍数，也就是说使用倍加变化是可以把第二行或第一行变为零行的，并且每有一个只存在一个非零元素的列，就会有一行被定死，情况非常多且复杂，我们必须考虑周全，首先我们不能让性质一被破坏，因此得先找出经过倍加变化会变为零行的行，把它们变为零行后再与最下方的非零行交换位置，这也是上方函数设计那个返回值的原因，下面来设计相应的算法：
  先看第一列，我们可以从上到下遍历，找出第一个非零元素并记录其值与行号，然后继续往下遍历，再次遍历到非零元素时对该行使用倍加变化，将其首元素变化成与记录元素相同的值后再与记录行号的行相减，相减后判断该行是不是零行，如果是，那么就与最下方的非零行交换位置，然后记录最下方非零的非零行号减一，要注意的是在交换后不能直接继续往下遍历，因为交换过来的非零行也得检测，直到检测到进行了倍加变化后还是非零行的行，那么才能够继续往下遍历，直到遍历到最下方的非零行。对于上方的算法，考虑下面两个矩阵：
$$\left[\begin{array}{ccc}0& 2& 1\\ 0& 4& 2\\ 0& 2& 3\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}0& 2& 1\\ 0& 4& 2\\ 1& 2& 3\\ 0& 0& 0\end{array}\right]$$
  显然无法处理，左侧的矩阵的第一列元素全是零，而右侧的矩阵上面两行经过倍加变化会得出零行，但是上方的算法只能够遍历找到第三行，因此考虑第二版的算法，我们得同时利用行与列，首先遍历第一行，找到第一个非零元素所在的列，记录该元素并遍历该列，如果找到了其它行在该列拥有非零元素，那么就让第一行与该行进行倍加变化，如果变化后第一行变为零行，那么就与最下方的非零行交换，然后记录最下方非零行号的变量减一并使用相同的操作重新遍历交换上来的行，如果第一行没有变为非零行，那么继续往后遍历该列下方的元素，如果又找到了非零元素，那么就让第一行与该行再次进行倍加变化，直到遍历到零行，然后以相同的操作继续遍历第二行，总结成步骤就是：

步骤
1.从第一行开始遍历每一行，找到第一个非零元素所处的
2.遍历该列，对于遍历到的元素分为两种情况
2.1遍历到零元素，那么继续往下遍历
2.2遍历到非零元素，那么让上方行与该行进行倍加变化，倍数加的结果分为两种
2.2.1如果倍加后，上方行变为零行，那么让上方行与最下方的非零行交换位置，对于交换上来的非零行，采用相同的操作重新遍历一次
2.2.2如果倍加变化后上方行还是非零行，那么让遍历列的变量继续往下遍历，如果再次遇到非零元素，那么采用相同的操作进行处理
3.直到该列遍历到第一个非零行，停止遍历，然后继续对下一行执行相同的操作

  可以感觉到，时间复杂度高的离谱，但是有实现的可能性，我们先实现，然后再来考虑优化，函数的接口设计如下：

int ProcessMonZeroRows(matrix& m)

实现如下：

//处理所有非零行，将倍加变化后变为零行的行移动到下面，返回变化后矩阵最下方非零行的行号
int ProcessMonZeroRows(matrix& m)
{
    //先将零行全部移动到非零行的下方并得出最下方非零行的行号
    int down = SwapZeroRowAndNotZeroRow(m);
    //如果交换后第零行是最下方的非零行，那么直接返回
    if(down == 0)
    {
        return 0;
    }

    int row_sz = m.size();
    int col_sz = m[0].size();

    //从第一行开始遍历每一行，找到第一个非零元素所处的列
    for(int rowup = 0; rowup < row_sz; rowup++)
    {
        for(int col = 0; col < col_sz; col++)
        {
            if(m[rowup][col])  
            {
                double upnum = m[rowup][col];
                for(int rowdown = rowup + 1; rowdown <= down; rowdown++)
                {
                    if(m[rowdown][col])   
                    {
                        double downnum = m[rowdown][col];
                        //遍历到非零元素，那么让上方行与该行进行倍加变化，倍数加的结果分为两种
                        RowAddition(m, rowdown, -upnum/downnum, rowup);
                        if(IsZeroRow(m, rowup))
                        {
                            //如果倍加后，上方行变为零行，那么让上方行与最下方的非零行交换位置，对于交换上来的非零行，采用相同的操作重新遍历一次
                            SwapRow(m, rowup, down);
                            down--;
                            //上方行成为最下方的非零行，退出函数
                            if(down == rowup)
                            {
                                return down;
                            }
                            else
                            {
                                rowup--;
                                goto out;
                            }
                        }
                    }
                }
            }
        }
out:
    }
    return down;
}

函数非常复杂，得多设计几组测试用例，测试用例如下：
$$\left[\begin{array}{ccc}0& 2& 1\\ 0& 4& 2\\ 0& 2& 3\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}0& 2& 1\\ 0& 4& 2\\ 1& 2& 3\\ 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}1& 0& 3& 4\\ 2& 0& 6& 8\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$
先测试函数，然后再通过手算确定结果是否正确，测试代码如下：

void TestProcessMonZeroRows()
{
    matrix m1 = {{0, 2, 1}, {0, 4, 2}, {0, 2, 3}, {0, 0, 0}};
    std::cout << "m1变化前: " << std::endl;
    ShowMatrix(m1);

    int d1 = ProcessMonZeroRows(m1);

    std::cout << "m1变化后: " << std::endl;
    ShowMatrix(m1);
    std::cout << "m1最下方非零行号: " << d1 << std::endl;

    matrix m2 = {{0, 2, 1}, {0, 4, 2}, {1, 2, 3}, {0, 0, 0}};
    std::cout << "m2变化前: " << std::endl;
    ShowMatrix(m2);

    int d2 = ProcessMonZeroRows(m2);

    std::cout << "m2变化后: " << std::endl;
    ShowMatrix(m2);
    std::cout << "m2最下方非零行号: " << d2 << std::endl;

    matrix m3 = {{1, 0, 3, 4}, {2, 0, 6, 8}, {0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}};
    std::cout << "m3变化前: " << std::endl;
    ShowMatrix(m3);

    int d3 = ProcessMonZeroRows(m3);

    std::cout << "m3变化后: " << std::endl;
    ShowMatrix(m3);
    std::cout << "m3最下方非零行号: " << d3 << std::endl;
}

结果如下：

  首先可以发现，能够通过倍加变化变为零行的行确实变为零行了，并且也没有破坏非零行必须位于零行上方的性质，验证正确性的方法是利用行初等变化的可逆性，如果变化后的矩阵能够通过行初等变化变为变化前的矩阵，那么本次变化就是正确的，具体的验证过程就不写出来水字数了，读者可以在草稿纸上自行验证，经过笔者的验证，上方的三个变化都是正确的，下面再来看一个测试矩阵：
$$\left[\begin{array}{ccc}4& 2& 1\\ 1& 3& 9\\ 5& 1& 3\\ 8& 4& 2\end{array}\right]$$
测试结果如下：

  答案错误了，是函数出bug了吗，应该说是也不是，我们遇到了在数据处理中一个非常恶心的问题，那就是小数精度的丢失，众所周知，计算机是无法精确的表示出无限不循环小数与无限循环小数的，应该说这是不可能被精确表示的，几次运算还不明显，但是运算次数越多，精度丢失就越严重，自然就无法得出正确的结果了，在此处采取一个简单的处理方案，那就是在进行了倍加变化后，如果没有得出零行，那么就使用方向的倍加变化把原行给变回去，修改后的代码如下：

//处理所有非零行，将倍加变化后变为零行的行移动到下面，返回变化后矩阵最下方非零行的行号
int ProcessMonZeroRows(matrix& m)
{
    //先将零行全部移动到非零行的下方并得出最下方非零行的行号
    int down = SwapZeroRowAndNotZeroRow(m);
    //如果交换后第零行是最下方的非零行，那么直接返回
    if(down == 0)
    {
        return 0;
    }

    int row_sz = m.size();
    int col_sz = m[0].size();

    //从第一行开始遍历每一行，找到第一个非零元素所处的列
    for(int rowup = 0; rowup < row_sz; rowup++)
    {
        for(int col = 0; col < col_sz; col++)
        {
            if(m[rowup][col])  
            {
                double upnum = m[rowup][col];
                for(int rowdown = rowup + 1; rowdown <= down; rowdown++)
                {
                    if(m[rowdown][col])   
                    {
                        double downnum = m[rowdown][col];
                        //遍历到非零元素，那么让上方行与该行进行倍加变化，倍数加的结果分为两种
                        RowAddition(m, rowdown, -upnum/downnum, rowup);
                        if(IsZeroRow(m, rowup))
                        {
                            //如果倍加后，上方行变为零行，那么让上方行与最下方的非零行交换位置，对于交换上来的非零行，采用相同的操作重新遍历一次
                            SwapRow(m, rowup, down);
                            down--;
                            //上方行成为最下方的非零行，退出函数
                            if(down == rowup)
                            {
                                return down;
                            }
                            else
                            {
                                rowup--;
                                goto out;
                            }
                        }
                        else
                        {
                            //防止精度丢失!!!
                            RowAddition(m, rowdown, upnum/downnum, rowup);
                        }
                    }
                }
            }
        }
out:
    }
    return down;
}

再来看看结果：

  现在就正确了，但是这个方法是治标不治本的，当数据量大的离谱的时候，每次矩阵变化都会出现一部分的精度丢失，不过现在先不管，等到学习了后面的线性代数知识后再来解决该问题。来回顾一下，我们的目的是实现一个函数，该函数能够将传入的任意矩阵变为行阶梯型矩阵，我们已经完成的步骤是：1.将所有零行移动到非零行的下方，2.处理了经过倍加变化会变为零行的行。那么下一步得思考一个问题：对于剩下的非零行，是否每一行都存在先导元素？在线性代数中我们得知：先导元素就是行阶梯型矩阵中非零行中的第一个非零元素，目前对于剩下的所有非零行，都无法使用倍加变化使其变为零行，并且显然对换变化和倍乘变化也无法使一个非零行变为零行，这意味着对于剩下的所有非零行使用任意的行初等变化，都无法使其变为零行，而一个矩阵通过行初等变化又一定能够变为行阶梯型矩阵，因此对于剩下的非零行，使用行初等变化一定能够使其变为阶梯型并且保持非零性质，所以可得剩下的行就是最终行阶梯矩阵中的所有非零行，每一行都存在先导元素，而我们的下一步就是把使用行初等变化把剩余行变化为阶梯型，然后就可以得出结果了，考虑下方矩阵：
$$\left[\begin{array}{cccc}9& 0& 2& 0\\ 1& 0& 3& 4\\ 1& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$
  显然上方的三行就是最终行阶梯矩阵的非零行，每一行都存在先导元素，有一列只存在一个元素，那么会出现某一行被定死的情况吗？答案是不会，因为使用倍加变化能够让该列其它行的元素也变为非零的，比如上方的矩阵，让第二行减去第三行，那么第二行的第二个元素就会变为非零的，因此我们可以使用按列遍历，从第一列开始，找出其中的一个非零元素，如果该列所有元素都为零，那么就继续遍历下一列，对于找出非零元素的那行，我们将其移动到第一行，让该元素变为第一行的先导元素，然后再使用倍加变化把该元素所在列下方的所有元素全部变为零，在操作结束后，第一行固定，因此在遍历第二列时得从第二行开始，如果找到了非零元素，那么也是移动到第二行，后面的所有行同理。
  在查找某一列的非零元素时，不一定就得找第一个遍历到的非零元素，比如上方的矩阵，第一列存在三个非零元素，那么这些非零元素所在的行都是可以移动到第一行的，我们可以选择任意一个元素，但是我们应该尽可能的减小精度损失，考虑下方矩阵：
$$\left[\begin{array}{cc}9& 0\\ 0.3& -9\\ 1& 1\\ 5& 0.01\\ 3& 2\\ 0& -0.001\end{array}\right]$$
  首先对于精度损失得认识到一点：使用一个很小的精度损失乘一个很大的数，那么精度损失就会被放大，比如我们在计算机中存储$1/3$，假设计算机只能够存储小数点后3位，那么$1/3$在计算机中存储的值就是$0.333$，失去的$\mathrm{0.0003333....}$就是精度损失，那么如果使用$100000$乘$1/3$，原本得出的精确结果应该是$100000/3$，整数部分是$33333$，但是由于在计算机中只存储了小数点后三位的数，因此的实际计算的结果是$33300$，发现了吗，精度的损失从$\mathrm{0.000333...}$变为了$33$，如果后面继续乘呢？那么损失就会继续被放大，最终导致程序出现严重偏差，以此为指导思想，我们选择的元素必须尽可能避免这种情况。
  对于上方的矩阵，在第一列中如果我们选择最小值$0.3$作为交换到第一行的元素，那么就得把另外四个元素通过倍加变化变为零，比如对于元素1就是：$1-0.3\ast (1/0.3)$，显然精度损失被放大了，因为$1/0.3=(1/3)\ast 10$，精度损失放大了十倍，再看看选择最大元素9的情况，此时对于元素1就是：$1-9\ast (1/9)$，我们发现，选取成为主元的元素总是会在倍加变化中变为除数，而除数越小，在除法中的精度损失就会被放的越大，下面来看看第二列的情况，此时出现了负数，如果我们选择最小元素$-9$作为主元，那么对于元素1就是：$1+9\ast (-1/9)$显然没有出现精度损失放大的情况，由此可得，我们在选取列中的一个元素作为主元时，要选择绝对值最大的那一个元素，该方法也叫做列主元法。
  由此实现我们的算法，具体步骤如下：

步骤
1.从第一列开始遍历每一列，遍历的情况分为两种
1.1.情况1：遍历列中没有非零元素，那么从同一行开始继续遍历下一列
1.2.情况2：遍历列中存在非零元素，那么先找出该列绝对值最大的非零元素，然后将该元素所在行与遍历列的起始行进行交换
1.2.1.交换后，使用倍加变化将该主元所在列下方的所有元素全部变为零，然后从下一行开始继续使用相同的操作遍历下一列

函数实现如下：

void RowChange::PivotDiagonal(matrix &m)
{
    int down = ProcessMonZeroRows(m);
    // 第一行就是最下方的非零行，直接就是行阶梯矩阵，因此直接返回
    if (down == 0)
    {
        return;
    }

    int row_sz = m.size();
    int col_sz = m[0].size();

    // 遍历每一列，最初从第一行开始
    int start_row = 0;

    for (int col = 0; col < col_sz; col++)
    {
        int flag = 0; // 标记在本列是否遍历到了非零元素
        double abs_max = 0;
        double val_max = 0;
        int abs_max_row = -1;
        for (int row = start_row; row <= down; row++)
        {
            // 遍历到非零元素
            if (m[row][col])
            {
                flag = 1;
                double tmp = abs_max;
                abs_max = std::max(abs_max, std::abs(m[row][col]));
                // abs_max值改变，说明新元素的绝对值更大，因此将abs_max_row修改成新元素所在行
                if (abs_max != tmp)
                {
                    val_max = m[row][col];
                    abs_max_row = row;
                }
            }
        }
        // 不是全零列，进行交换操作
        if (flag)
        {
            MatrixUtil::SwapRow(m, start_row, abs_max_row);
            abs_max_row = start_row; // 交换后，原本的start_row成为拥有绝对值最大的行
            start_row++;             // 本行被确定了，下一次遍历列从下一行开始
            // 使用倍加变化将下方元素全部变为零
            for (int row = start_row; row <= down; row++)
            {
                RowAddition(m, abs_max_row, -m[row][col] / val_max, row);
            }
        }
    }
}

测试用例如下：
$$\left[\begin{array}{ccc}0& 2& 1\\ 0& 4& 2\\ 0& 0& 0\\ 0& 8& 4\end{array}\right]\left[\begin{array}{ccc}4& 2& 1\\ 2& 4& 2\\ 1& 2& 3\\ 3& 5& 7\end{array}\right]\left[\begin{array}{cccc}1& 0& 3& 4\\ 2& 0& 6& 8\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$
测试代码如下：

static void TestPivotDiagonal()
{
    std::cout << "-------------------------------------------------" << std::endl;
    std::cout << "TestPivotDiagonal: " << std::endl;
    matrix m1 = {{0, 2, 1}, {0, 4, 2}, {0, 0, 0}, {0, 8, 4}};
    std::cout << "m1变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m1);

    _rc.PivotDiagonal(m1);

    std::cout << "m1变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m1);

    matrix m2 = {{4, 2, 1}, {2, 4, 2}, {1, 2, 3}, {3, 5, 7}};
    std::cout << "m2变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m2);

    _rc.PivotDiagonal(m2);

    std::cout << "m2变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m2);

    matrix m3 = {{1, 0, 3, 4}, {2, 0, 6, 8}, {0, 1, 0, 0}, {1, 0, 0, 0}};
    std::cout << "m3变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m3);

    _rc.PivotDiagonal(m3);

    std::cout << "m3变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m3);
}

结果如下：

  答案不知道对不对，但是形状是对了，同样的可以利用初等行变化的可逆性进行验证，看到m2矩阵，最后一行变为了零行，但是显然对于m2原矩阵的任意两行进行倍加变化都无法得出零行，也就是说使用ProcessMonZeroRows函数是无法变化出零行的，测试如下：

static void TestProcessMonZeroRows()
{
    std::cout << "-------------------------------------------------" << std::endl;
    std::cout << "TestProcessMonZeroRows: " << std::endl;
    matrix m = {{4, 2, 1}, {2, 4, 2}, {1, 2, 3}, {3, 5, 7}};
    std::cout << "变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m);

    _rc.ProcessMonZeroRows(m);

    std::cout << "变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m);
}

结果如下：

  显然没有出现零行，所以是我们的函数出bug了吗，显然不是的，上方的矩阵总共就只有四行，那么要变为行阶梯型矩阵就必须有一行变为零行，因此在上方我们的推理是有些问题的，在经过了ProcessMonZeroRows函数变化后的非零行使用行初等变化也是有可能变化出非零行的，在上方的推理中只考虑了两行一一对应进行倍加变化的情况，但是同一行是可以与多行进行倍加变化的，因此是可能会出现零行的，不过我们的行阶梯型矩阵变化函数误打误撞的处理这个问题，因为在算法设计中是让每一行同时与多行进行倍加变化的，因此得出了正确的答案，这只能说是运气比较好，一般而言，在出现这种情况时都是要重新设计算法的，我们得汲取这次经验，在之后设计算法时必须尽可能的考虑到所有情况，目前我们就完成了解线性方程组中的前向步骤，也就是把矩阵变化为行阶梯型矩阵，其实笔者的第一章中设计的行化简矩阵算法也是有问题的，只能够处理单一情况，因此重新设计新的后向步骤的算法，实现更加普适的函数。下面来验证一下之前的行化简函数是真的有bug的，第一节的行化简函数如下：

void ClearCol(matrix& m)
{
    int row_size = m.size();
    int column_size = m[0].size();

    for(int i = 0; i < row_size; i++)
    {
        // 1.从第一列开始，确定对角线元素是否为1，如果不是，那么就使用倍乘变化将其变为1
        if(m[i][i] != 1)
        {
            RowScaling(m, i, 1.0/m[i][i]);
        }
        // 2.将该列其它行的元素变为0
        for(int j = 0; j < row_size; j++)
        {
            if(m[j][i] && j != i)
            {
                RowAddition(m, i, -m[j][i], j);
            }
        }
    }
}

测试相同的矩阵，得出的结果如下：

  显然出bug了，nan表示数学运算无意义，比如0.0 除以 0.0，对负数进行开平方，对负数取对数等操作，之前设计的函数只考虑了对角线上的元素，不出bug才怪，下面重新设计新的算法，考虑下方矩阵：
$$\left[\begin{array}{cccc}4& 2& 1& 3\\ 0& 3& 9& 3\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right]$$
  我们要把这个行阶梯型矩阵化为行化简矩阵，算法基础就是那两条行化简矩阵多出来的性质：

性质
1.先导元素必须为1
2.先导元素所在列除了先导元素之外的元素必须全部为零

关键就是先导元素所处的位置，也就是主元位置，在成功将矩阵化为行阶梯型矩阵后，我们可以使用本节实现的接口来查找，也就是这个函数：

ssize_t FindOneRowPP(matrix& m, size_t row, size_t col_sz)
{
    for (int i = 0; i < col_sz; i++)
    {
        if (m[row][i])
        {
            return i;
        }
    }
    return -1;
}

ssize_t FindPivotPosition(matrix &m, size_t row)
{
    int row_sz = m.size();
    int col_sz = m[0].size();
    if (row > row_sz)
    {
        return -1;
    }
    if (IsREF(m))
    {
        return FindOneRowPP(m, row, col_sz);
    }
    matrix tmp_m = m;
    PivotDiagonal(m);
    return FindOneRowPP(m, row, col_sz);
}

在修改了行阶梯型矩阵变化函数后，让我们重新来测试一下该函数，测试代码如下：

void TestFindPivotPosition()
{
    matrix m = {{1, 3, 5}, {2, 1, 3}, {3, 4, 5}, {8, 1, 9}};
    for(int i = 0; i < m.size(); i++)
    {
        int tmp = FindPivotPosition(m, i);
        if(tmp != -1)
        {
            std::cout << "第" << i << "行" << "主元位置为: " << tmp << std::endl;
        }
        else{
             std::cout << "第" << i << "行" << "不存在主元位置" << std::endl;
        }
    }
}

结果如下：

结果正确，我们可以利用该函数先找出每一行的主元，然后再把使用倍乘变化把主元变为1，最后再使用倍加变化把其它行中该列的元素全部变为零，具体的步骤就是：

步骤
1.从第一行开始找出每一行主元所在列
2.使用倍乘变化将主元变为一
3.使用倍加变化将同列的其它元素全部变为零
函数实现如下：

void RowChange::ChangeMatrixToClearCol(matrix &m)
{
    // 如果矩阵不是行阶梯型矩阵，那么就把其变为行阶梯型矩阵
    if (!MatrixUtil::IsREF(m))
    {
        PivotDiagonal(m);
    }

    int row_sz = m.size();

    // 从第一行开始找出每一行的主元
    for (int row = 0; row < row_sz; row++)
    {
        int pcol = FindPivotPosition(m, row);
        // 遍历到零行，退出遍历
        if (pcol == -1)
        {
            break;
        }

        // 使用倍加变化将同列的其它元素全部变为零，由行阶梯型矩阵的性质可知：只有上方的元素可能是非零的，下方元素必定为零
        double tmp = m[row][pcol];
        if (std::abs(m[row][pcol]) > 1)
        {
            // 考虑到精度损失，在倍加变化时要尽可能的让除数的绝对值较大
            for (int i = 0; i < row; i++)
            {
                RowAddition(m, row, -m[i][pcol] / m[row][pcol], i);
            }
            // 使用倍乘变化将主元变为一
            RowScaling(m, row, 1 / m[row][pcol]);
        }
        else
        {
            RowScaling(m, row, 1 / m[row][pcol]);
            for (int i = 0; i < row; i++)
            {
                RowAddition(m, row, -m[i][pcol] / m[row][pcol], i);
            }
        }
    }
}

测试代码如下：

static void TestChangeMatrixToClearCol()
{
    std::cout << "-------------------------------------------------" << std::endl;
    std::cout << "TestPivotDiagonal: " << std::endl;
    matrix m1 = {{0, 2, 1}, {0, 4, 2}, {0, 0, 0}, {0, 8, 4}};
    std::cout << "m1变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m1);

    _rc.ChangeMatrixToClearCol(m1);

    std::cout << "m1变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m1);

    matrix m2 = {{4, 2, 1}, {2, 4, 2}, {1, 2, 3}, {3, 5, 7}};
    std::cout << "m2变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m2);

    _rc.ChangeMatrixToClearCol(m2);

    std::cout << "m2变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m2);

    matrix m3 = {{1, 0, 3, 4}, {2, 0, 6, 8}, {0, 1, 0, 0}, {1, 0, 0, 0}};
    std::cout << "m3变化前: " << std::endl;
    MatrixUtil::ShowMatrix(m3);

    _rc.ChangeMatrixToClearCol(m3);

    std::cout << "m3变化后: " << std::endl;
    MatrixUtil::ShowMatrix(m3);
}

结果如下：

  可以使用手算验证一下，经过笔者的验证，结果是正确的，对于小数部分，只保留了一位，我们已经实现了解线性方程组的向前步骤和向后步骤，由于没有经过大型矩阵的大量测试，笔者不敢保证上方的函数就是完全正确的，或者应该说大概率是有bug的，不过本次是算法设计是严格的以行阶梯矩阵与行化简矩阵的性质为依据的，单说普适性肯定是提升了不少的。在实现本节剩下的两个函数前，先来看看行阶梯型矩阵变化函数与行化简矩阵变化函数的时空复杂度，首先分析行阶梯型矩阵变化函数：
  在书中说前向步骤要在计算机中要经过更多的计算，在设计了具体的算法并实现后，这是非常明显的，前向步骤的算法要比后向步骤复杂的多。我们先来分析空间复杂度，首先是前向步骤的第一步，把零行移动到非零行下方，也就是函数SwapZeroRowAndNotZeroRow，在该函数中并没有任何额外的空间开销，空间复杂度显然是$O(1)$的，下面来看看时间复杂度，首先进入函数就是一个IsREF判断，在上一节分析过了，IsREF的时间复杂度是$O(n^2)$的，然后是下方的循环，在IsZeroRow中是存在一个循环的，加上外层的循环，因此循环的时间复杂度是$O(n^2)$的，那么函数整体的时间复杂度就是$O(n^2)$的，然后是函数ProcessMonZeroRows，其实可以发现，这个函数是没有必要的，完全可以在PivotDiagonal函数倍加变化时加上一步对换变化来达成相同的效果，改版后的PivotDiagonal函数如下：

void RowChange::PivotDiagonal(matrix &m)
{
    int down = SwapZeroRowAndNotZeroRow(m);
    // 第一行就是最下方的非零行，直接就是行阶梯矩阵，因此直接返回
    if (down == 0)
    {
        return;
    }

    int row_sz = m.size();
    int col_sz = m[0].size();

    // 遍历每一列，最初从第一行开始
    int start_row = 0;

    for (int col = 0; col < col_sz; col++)
    {
        int flag = 0; // 标记在本列是否遍历到了非零元素
        double abs_max = 0;
        double val_max = 0;
        int abs_max_row = -1;
        for (int row = start_row; row <= down; row++)
        {
            // 遍历到非零元素
            if (m[row][col])
            {
                flag = 1;
                double tmp = abs_max;
                abs_max = std::max(abs_max, std::abs(m[row][col]));
                // abs_max值改变，说明新元素的绝对值更大，因此将abs_max_row修改成新元素所在行
                if (abs_max != tmp)
                {
                    val_max = m[row][col];
                    abs_max_row = row;
                }
            }
        }
        // 不是全零列，进行交换操作
        if (flag)
        {
            MatrixUtil::SwapRow(m, start_row, abs_max_row);
            abs_max_row = start_row; // 交换后，原本的start_row成为拥有绝对值最大的行
            start_row++;             // 本行被确定了，下一次遍历列从下一行开始
            // 使用倍加变化将下方元素全部变为零
            for (int row = start_row; row <= down; row++)
            {
                RowAddition(m, abs_max_row, -m[row][col] / val_max, row);
                // 如果变化出零行，那么与最下方非零行交换
                if (MatrixUtil::IsZeroRow(m, row))
                {
                    RowSwapping(m, row, down);
                    down--;
                    // ！重新处理交换上来的非零行
                    row--;
                }
            }
        }
    }
}

  笔者测试过了，效果是一模一样的，这是一次质变的优化，直接把最复杂，时间复杂度最高的函数给去掉了，变成了一个简单的if判断，所以才说算法是迭代出来的，当你设计到后面，自然会发现前面的设计缺陷，在使用时自然也会不断的进行优化调整，所以对于一个算法，在初期最重要的永远是先实现，让代码跑起来，就算复杂度高的离谱也无所谓，随着设计的不断深入，不断的测试，自然就会去调整优化了，在脑子里面胡思乱想，不真正的上手实现一个可以跑的函数，是想不出好算法的。
  下面来分析一下PivotDiagonal的时空复杂度，首先是空间复杂度，函数中调用了RowAddition，而RowAddition的空间复杂度是$O(n)$的，在其它地方并没有额外的空间开销，并且在RowAddition也是临时使用，马上释放，因此函数整体的空间复杂度是$O(n)$的，然后是时间复杂度，首先进入函数就是一个SwapZeroRowAndNotZeroRow函数调用，该函数的时间复杂度是$O(n^2)$的，然后进入循环，最深嵌套了三层循环，每层在最坏情况下都是遍历完整的一行或一列，因此循环中的时间复杂度是$(O^3)$的，函数整体的时间复杂度就是$O(n^3)$，有点丑陋，后面再优化吧，现在先用着，下面来看看最后ChangeMatrixToClearCol函数的情况：
  首先在函数中调用了PivotDiagonal，时间复杂度最低也是$O(n^3)$的，然后开始循环，最深嵌套同样是三层循环，并且在最坏情况下同样是每层遍历完一行或一列，因此时间复杂度是$O(n^3)$的，函数整体的时间复杂度就是$O(n^3)$，然后是空间复杂度，同样只有RowAddition消耗了空间，因此空间复杂度是$O(n)$的。
  终于，我们填平了第一节挖的坑，变化行阶梯型矩阵函数与变化行化简矩阵函数可以说是整个线性地址中最基础，最重要的步骤，必须要好好的设计，下面来实现本节的函数，还有两个没有实现：

函数
1.给定一个矩阵并指定列，判断该列是否是主元列
2.给定一个矩阵，判断矩阵解的情况

  先来看看第二个，显然只需要先把矩阵变化为行阶梯型矩阵，然后再判断该列是否存在主元位置就可以了，具体步骤如下：

步骤
1.将每一行投入主元列查找函数，找出主元所在列
2.比较指定列与主元所在列，得出结果
函数实现如下：

bool RowChange::IsPivotCol(matrix m, size_t col)
{
    // 获取最下方非零行行号
    int down = MatrixUtil::DownNotZrorRow(m);
    int col_sz = m[0].size();

    if (col > col_sz || down == -1)
    {
        return false;
    }
	
	//将矩阵变为行阶梯型矩阵
    PivotDiagonal(m);

    // 将每一行投入主元列查找函数，找出主元所在列
    for (int row = 0; row <= down; row++)
    {
        // 比较指定列与主元所在列，得出结果
        if (col == FindPivotPosition(m, row))
        {
            return true;
        }
    }
    return false;
}

测试代码如下：

static void TestIsPivotCol()
{
    matrix m1 = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 0}};
    std::cout << "矩阵m1: " << std::endl;
    MatrixUtil::ShowMatrix(m1);
    for (int i = 0; i < m1[0].size(); i++)
    {
        if (_rc.IsPivotCol(m1, i))
        {
            std::cout << "m1第" << i << "列是主元列" << std::endl;
        }
    }

    matrix m2 = {{1, 4, 3, 6, 9}, {3, 1, 2, 4, 5}, {2, 3, 4, 5, 9}, {1, 2, 2, 3, 4}};
    std::cout << "矩阵m2: " << std::endl;
    MatrixUtil::ShowMatrix(m2);
    for (int i = 0; i < m2[0].size(); i++)
    {
        if (_rc.IsPivotCol(m2, i))
        {
            std::cout << "m2第" << i << "列是主元列" << std::endl;
        }
    }
}

结果如下：

经笔者验算，上方的结果是正确的，该函数不改变原矩阵，因此使用拷贝传参，空间复杂度就是$O(n^2)$的，在函数中调用了PivotDiagonal，时间复杂度是$O(n^3)$的，其下方的循环显然时间复杂度是$O(n^2)$的，因此总体的时间复杂度就是$O(n^3)$，最后来实现函数：给定一个矩阵，判断矩阵解的情况，对于该函数，显然要利用上方在解的存在性与唯一性中学到的线性代数知识，首先要把矩阵化为行阶梯型矩阵，然后判断最后一列是否是主元列，如果是，那么矩阵对应的线性方程组无解，如果不是，那么线性方程组有解，在有解的情况下，将除最后一列的每一列都进行判断，如果每一列都是主元列，那么有唯一解，否则就是无数解，具体步骤入下：

步骤
1.将矩阵化为行阶梯型矩阵
2.判断最后一列是否是主元列
2.1如果是，那么得出无解
2.2如果不是，那么再次分情况讨论
2.2.1将其余每一列都进行判断，如果全部都是主元列，那么存在唯一解，否则存在无数解
函数接口如下：

enum class SOLUTION
{
     ONE_SOLUTION,     //唯一解
     MYRIAD_SOLUTION,  //无数解
     NONE_SOLUTION     //无解
};
PUBLIC_MOD::SOLUTION LinearEquations(matrix m);        

函数实现如下：

PUBLIC_MOD::SOLUTION RowChange::LinearEquations(matrix m)
{
    // 将矩阵化为行阶梯型矩阵
    PivotDiagonal(m);

    int row_sz = m.size();
    int col_sz = m[0].size();
    // 判断最后一列是否是主元列
    if (IsPivotCol(m, col_sz - 1))
    {
        return PUBLIC_MOD::SOLUTION::NONE_SOLUTION;
    }
    else
    {
        // 将其余每一列都进行判断，如果全部都是主元列，那么存在唯一解，否则存在无数解
        for (int col = 0; col < col_sz - 1; col++)
        {
            if (!IsPivotCol(m, col))
            {
                return PUBLIC_MOD::SOLUTION::MYRIAD_SOLUTION;
            }
        }
        return PUBLIC_MOD::SOLUTION::ONE_SOLUTION;
    }
}

测试矩阵如下：
$$(1)\left[\begin{array}{ccc}1& 0& 1\\ 0& 4& 2\\ 0& 0& 0\\ 0& 8& 4\end{array}\right](2)\left[\begin{array}{ccc}4& 2& 1\\ 2& 4& 2\\ 1& 2& 3\\ 3& 5& 7\end{array}\right](3)\left[\begin{array}{cccc}1& 0& 3& 4\\ 0& 1& 7& 0\end{array}\right]$$
先来手算一下，第一个矩阵显然只有唯一解，第三个矩阵显然有无数解，第二个矩阵变化行阶梯矩阵后情况如下：
$$(2)\left[\begin{array}{ccc}1& 2& 1\\ 0& -1& 4\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]$$
显然是无解的，测试代码如下:

static void TestLinearEquations()
{
    matrix m1 = {{1, 0, 1}, {0, 4, 2}, {0, 0, 0}, {0, 8, 4}};
    if (_rc.LinearEquations(m1) == PUBLIC_MOD::SOLUTION::MYRIAD_SOLUTION)
    {
        std::cout << "m1有无数解" << std::endl;
    }
    else if (_rc.LinearEquations(m1) == PUBLIC_MOD::SOLUTION::NONE_SOLUTION)
    {
        std::cout << "m1无解" << std::endl;
    }
    else
        std::cout << "m1有唯一解" << std::endl;

    matrix m2 = {{4, 2, 1}, {2, 4, 2}, {1, 2, 3}, {3, 5, 7}};
    if (_rc.LinearEquations(m2) == PUBLIC_MOD::SOLUTION::MYRIAD_SOLUTION)
    {
        std::cout << "m2有无数解" << std::endl;
    }
    else if (_rc.LinearEquations(m2) == PUBLIC_MOD::SOLUTION::NONE_SOLUTION)
    {
        std::cout << "m2无解" << std::endl;
    }
    else
        std::cout << "m2有唯一解" << std::endl;

    matrix m3 = {{1, 0, 3, 4}, {0, 1, 7, 0}};
    if (_rc.LinearEquations(m3) == PUBLIC_MOD::SOLUTION::MYRIAD_SOLUTION)
    {
        std::cout << "m3有无数解" << std::endl;
    }
    else if (_rc.LinearEquations(m3) == PUBLIC_MOD::SOLUTION::NONE_SOLUTION)
    {
        std::cout << "m3无解" << std::endl;
    }
    else
        std::cout << "m3有唯一解" << std::endl;
}

结果如下：

结果正确，最后再来额外实现一个函数，该函数接收一个增广矩阵，返回自由变量的个数，如果不存在自由变量，那么返回-1，步骤如下：

步骤
1.判断矩阵的解是否是无数解的情况，如果不是，直接返回-1
2.将矩阵除去最后一列的其余列减去主元列，得出的就是自由变量的个数

函数实现如下：

ssize_t RowChange::FreeXSize(matrix m)
{
    // 将矩阵化为行阶梯型矩阵
    PivotDiagonal(m);

    int row_sz = m.size();
    int col_sz = m[0].size();
    // 判断最后一列是否是主元列
    if (IsPivotCol(m, col_sz - 1))
    {
        return -1;
    }
    int sum = 0;
    for (int col = 0; col < col_sz - 1; col++)
    {
        if (IsPivotCol(m, col))
        {
            sum++;
        }
    }

    return sum == col_sz - 1 ? -1 : col_sz - 1 - sum;
}

没有复用代码是因为会多一次拷贝传参，测试矩阵如下：
$$(1)\left[\begin{array}{ccccc}1& 0& 1& 4& 8\\ 0& 4& 2& 2& 3\end{array}\right](2)\left[\begin{array}{ccc}4& 2& 1\\ 2& 4& 2\end{array}\right](3)\left[\begin{array}{cccc}1& 0& 3& 4\\ 0& 1& 7& 0\end{array}\right]$$
显然第一个线性方程组的解中存在两个自由变量，解是一个平面，第二个线性方程组的解中不存在自由变量，解是一个点，第三个线性方程组的解中存在一个自由变量，解是一条直线，测试代码如下：

static void TestFreeXSize()
{
    matrix m1 = {{1, 0, 1, 4, 8}, {0, 4, 2, 2, 3}};
    int j1 = _rc.FreeXSize(m1);
    if (j1 != -1)
    {
        std::cout << "m1自由变量个数: " << j1 << std::endl;
    }

    matrix m2 = {{4, 2, 1}, {2, 4, 2}};
    int j2 = _rc.FreeXSize(m2);
    if (j2 != -1)
    {
        std::cout << "m2自由变量个数: " << j2 << std::endl;
    }

    matrix m3 = {{1, 0, 3, 4}, {0, 1, 7, 0}};
    int j3 = _rc.FreeXSize(m3);
    if (j3 != -1)
    {
        std::cout << "m3自由变量个数: " << j3 << std::endl;
    }
}

结果如下：

结果正确。
  至此本小节结束，其实明显能够感觉到代码是有些bug的，因此在下一节中将整体性的优化目前实现的所有函数，同样的笔者在写完后感觉收获巨大，愈发的体会到数学定理对算法重要的指导意义，也希望这篇文章能够给读者一些收获。