稀疏矩阵算法Envelope方法三之Cholesky分解

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

ESFCT

原理

ESFCT，用于计算 Cholesky 分解 ，将 矩阵按照 envelope storage scheme 格式存储。

Cholesky 分解的方法是 bordering method ，假设矩阵 可以分块为

其中 可以分解为 ，因此 的分解因子为

其中 ， ，因此 的分解因子 可以逐行计算。以下图为例，假设我们已经完成了前 步的分解，因此 的左上角的 的子矩阵已经完成了分解。为了计算 矩阵的第 行，我们需要计算 。但是，实际上 只有一部分 envelope 区域需要计算，因此实际使用和计算的 和 如右图，因此此时 ELSLV 可以对这个 IBAND 大小的三角方程组进行求解

脚本

具体脚本如下：

现在开始逐行解析。

首先对第一行进行分解，如果对角线元素小于等于0则程序中止

IF ( DIAG(1) .LE. 0.0E0) GO TO 400
DIAG(1) = SQRT(DIAG(1))
IF( NEQNS .EQ. 1) RETURN

对第二行及之后的行进行分解，IXENV 是 ENV 向量中第 行第一个元素的位置，TEMP 是矩阵第 行的对角线元素，如果 IBAND 等于 0，则跳过三角方程组求解这一步，GO TO 200 。

IFIRST 是第 行第一个元素所在列。

之后调用 ELSLV 求解得到 。 这里就用到了子矩阵求解的方法，方程数目为 IBAND，XENV(IFIRST) 表示其第一个未知数是 IFIRST所在行/列的未知数，DIAG(IFIRST)，其右手项是本行的非零元素，第一个元素为 ENV(IXENV) 。

因此调用结束后会更新右手项的值，即更新 ENV 向量

DO 300 I = 2, NEQNS
  IXENV = XENV(I)
  IBAND = XENV(I+1) - IXENV
  TEMP = DIAG(I)
  IF (IBAND .EQ. 0) GO TO 200
    IFIRST = I - IBAND
    CALL ELSLV( IBAND, XENV(IFIRST), ENV, DIAG(IFIRST), ENV(IXENV))

下一步，我们计算 t ，首先我们从对角线元素中减去 的内积。

JSTOP 是 ENV 向量中第 行最后一个元素的位置，DO 100 循环遍历了ENV 向量中第 行的全部位置，TEMP = TEMP - S * S 便是对角线元素减去了 的内积

	      JSTOP = XENV(I+1) - 1
	      DO 100 J = IXENV, JSTOP
	      	S = ENV(J)
	      	TEMP = TEMP - S * S
100	      CONTINUE

然后求开方，如果 TEMP 小于等于0，程序中止。

200    	IF(TEMP .LE.  0.0E0) GO TO 400
		  DIAG(I) = SQRT(TEMP)
		  COUNT = IBAND 
		  OPS = OPS + COUNT 
300	  CONTINUE
	  RETURN 

如果存在对角线元素小于等于0的情况，程序中止，IFLAG 设为 1

400   IFLAG = 1
	  RETURN
	  END

参考文献

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