稀疏矩阵算法Envelope方法二之三角方程组求解

本章节是介绍 George and Liu 书中 Envelope Method 的应用

在这一章节中，我们利用上面章节提到的 envelope storage scheme，描述执行 Cholesky 分解和求解的子程序。我们会先介绍三角方程组求解子程序 ELSLV (Envelope-Lower-Solve) 和 EUSLV (Envelope-Upper-Solve) ，再介绍矩阵分解的子程序 ESFCT (Envelope-Symmetric-Factorization) 。

ELSLV

这些子程序用于求解下面的下三角方程组和上三角方程组

其中 是一个下三角矩阵，其存储方式就是之前提到的 envelope storage scheme 。

ELSLV 脚本如下，该流程首先查找右手项 RHS 中的第一个非零元素 (IFIRST) ，之前的未知数就全部是 0。然后从 IFIRST 行循环至最后一行，采用内积原理，逐个计算 ，覆盖 RHS 相应元素的值。

变量 LAST 记录上一次计算的非零解的索引，这样如果本方程计算时左手项第一个非零元素大于 LAST ，那么就可以直接采用 计算解。

最终输出 RHS 和 OPS，这里 RHS 就是最终的解，OPS 是执行的算术操作的次数（乘法和除法）。我感觉OPS 应该采用长整型，而不是用双精度浮点数。

注意，下面脚本中的 ，其中的 1.0D0 就是 ， 是双精度的意思。

下面我们开始逐行解析脚本，定义变量时将 OPS 定义成了 SPKOPS 组的第一个全局变量，用于记录乘除法的操作次数。

COMMON /SPKOPS/ OPS

第一步，我们查找右手项中第一个非零元素的位置

如果右手项全部元素均为0，则程序终止，全为0的 RHS 就是方程组的解

	  IFIRST = 0
100	  IFIRST = IFIRST + 1
	    IF( RHS(IFIRST) .NE. 0.0E0 ) GO TO 200
	    IF ( IFIRST .LT. NEQNS ) GO TO 100
	    RETURN
200	  LAST = 0

从 IFIRST 行开始逐行求解 (之前的行解也全部为0)，采用内积原理进行求解。

IBAND 是这一行的 bandwidth ， IF ( IBAND .GE. I) IBAND = I - 1 这句目前看是废话，之后会有进一步的解释。

L = I - IBAND 是这一行第一个非零元素的位置；RHS(I) = 0.0E0 将对应的解首先设置为 0 ，之后再覆盖。

LAST 记录的是上一次计算的非零解的位置。

IF (IBAND .EQ. 0 .OR. LAST .LT. L) GO TO 400 ，这句话第一个条件很容易理解，如果 IBAND 等于 0 ，那么就只需要考虑对角线元素；第二个条件是 LAST 小于 L ，LAST 是上一次计算的非零解的位置，那就说明从 LAST + 1 到 I - 1，这些位置的解均为 0；而这一行第一个非零元素的位置 L 大于 LAST ，说明这一行非对角线元素对应的解均为0，因此也只需要考虑对角线元素。

如果不满足上面两个条件，接着往下走，KSTRT 和 KSTOP 是 IBAND 的第一个和最后一个元素，DO 300 循环根据内积原理从 S 中减去非对角线元素的乘积 (S = S - ENV(K)*RHS(L)) 。

最后，COUNT 是本次的乘除法计算次数，加入到 OPS 中

	  DO 500 I = IFIRST, NEQNS
	    IBAND = XENV(I+1) - XENV(I)
	    IF ( IBAND .GE. I) IBAND = I - 1
	    S = RHS(I)
	    L = I - IBAND
	    RHS(I) = 0.0E0
	    
	    IF (IBAND .EQ. 0 .OR. LAST .LT. L) GO TO 400
	      KSTRT = XENV(I+1) - IBAND
	      KSTOP = XENV(I+1) - 1
	      DO 300 K = KSTRT, KSTOP
	        S = S - ENV(K)*RHS(L)
	        L = L + 1
300	      CONTINUE
	      COUNT = IBAND
	      OPS = OPS + COUNT

最后计算 RHS(I) ，如果此时 S 已经等于 0 ，那么 RHS(I) = 0 , 这就是我们预设的值，不用做任何事情，直接求解下一个未知数；否则的话，RHS(I) = S/DIAG(I) ，OPS 加1，LAST 更新为 I 。

400     IF( S .EQ. 0.0E0) GO TO 500
		  RHS(I) = S/DIAG(I)
		  OPS = OPS + 1.0D0
		  LAST = I
500   CONTINUE
      RETURN
    END

下图描述了本子程序的计算逻辑，其中圆圈圈起来的元素是ELSLV 实际使用的元素，可以按照上面的子程序手动模拟一下计算流程，看看是不是只用了圆圈圈起来的元素。

注意，在这个子流程中，我们在 DO 300 循环中可能执行了 一些乘以 0 的乘法，但是在大部分机器中通过做一个检查来避免这种乘法可能比不管它更加耗时。

子程序中 54 行对 IBAND 的处理（ IF ( IBAND .GE. I) IBAND = I - 1 ）解释如下，在某些环境中，ELSLV 子程序只用于求解导入的系数矩阵**(XADJ, ADJNCY)的子矩阵**的解，如下图所示， 实际是一个 的矩阵，而我们想要求解的方程组是框起来的 的子矩阵，伴随着其相应的右手项 RHS ，我们求解的脚本为

在这个子程序中，XENV(5) 被解释为 XENV(1)， XENV(6) 被解释为 XENV(2)，依此类推。

EUSLV

原理

我们现在查看子程序 EUSLV ，其求解 ，其中 就是 ELSLV 中的系数矩阵，这就意味着我们可以获取 的列，此时我们采用 outer product form，该方法描述如下

假设 的方程组 ，其中 是一个非奇异的下三角矩阵，并且按列存储，其计算公式如下

这种算法适用于 向量很稀疏的情况，如果在第 步一开始 就等于 0 ，那么 等于 0，我们可以整个跳过这一步。

同样类似于 ELSLV ，我们可以对提供的系数矩阵的子矩阵进行计算求解。

其求解顺序是反的，从 到 。

脚本

现在对该脚本逐行解析如下，首先这个是按照从ｎ到１的顺序，设 I 的初值设为 NEQNS + 1，进入循环。

如果 RHS(I) 等于0，直接跳过，进入下一个循环

	  I = NEQNS + 1
100   I = I - 1
	  IF (I .EQ. 0) RETURN
	    IF (RHS(I) .EQ. 0.0E0) GO TO 100

不然的话，按照上面的公式进行处理，下面这两行将 RHS(I) 赋值为 RHS(I)/DIAG(I) ，OPS 加1

S = RHS(I)/DIAG(I)
RHS(I) = S
OPS = OPS + 1.0D0

接下来处理第二步， 的第 列，就是 矩阵的第 行，而且只需要针对第 行的非零元素。

IBAND 是第 行的 bandwidth ，KSTRT 和 KSTOP 是第 行相对于 IBAND 的第一个位置和最后一个位置。

L = XENV(I+1) - IBAND 其实就是 XENV(I)，就是 ENV 向量中第 行的第一个元素位置 。

DO 200 循环中，RHS(K) = RHS(K) - S*ENV(L) ，就是右手项相应元素减去 RHS(I) 乘以 第 行的相应元素。

	  IBAND = XENV(I+1) - XENV(I)
	  IF (IBAND .GE. I) IBAND = I-1
	  IF (IBAND .EQ. 0) GO TO 100
	    KSTRT = I - IBAND
	    KSTOP = I - 1
	    L = XENV(I+1) - IBAND
	    DO 200 K = KSTRT, KSTOP
	      RHS(K) = RHS(K) - S*ENV(L)
	      L = L + 1
	    CONTINUE
	    COUNT = IBAND
	    OPS = OPS + COUNT
	  GO TO 100
END

一直遍历到 为止，便得到了所有未知数的解。

参考文献

1. George A, Liu J, Ng E. Computer solution of sparse linear systems[J]. Oak Ridge National Laboratory, 1994.