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

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

ESFCT

原理

ESFCT，用于计算 Cholesky 分解 \(\mathbf{LL^{T}}\) ，将 \(\mathbf{L}\) 矩阵按照 envelope storage scheme 格式存储。

Cholesky 分解的方法是 bordering method ，假设矩阵 \(\mathbf{A}\) 可以分块为 \[ \boldsymbol{A}=\left(\begin{array}{ll} \boldsymbol{M} & \boldsymbol{u} \\ \boldsymbol{u}^{T} & s \end{array}\right) \] 其中 \(\mathbf{M}\) 可以分解为 \(\mathbf{L_{M}L_{M}^{T}}\) ，因此 \(\mathbf{A}\) 的分解因子为 \[ \boldsymbol{L}=\left(\begin{array}{ll} \boldsymbol{L_{M}} & \boldsymbol{0} \\ \boldsymbol{w}^{T} & t \end{array}\right) \] 其中 \(\mathbf{L_{M}w = u}\) ，\(\mathbf{t = (s - w^{T}w)^{1/2}}\) ，因此 \(\mathbf{A}\) 的分解因子 \(\mathbf{L}\) 可以逐行计算。以下图为例，假设我们已经完成了前 \(i-1\) 步的分解，因此 \(\mathbf{A}\) 的左上角的 \((i-1) \times (i-1)\) 的子矩阵已经完成了分解。为了计算 \(\mathbf{L}\) 矩阵的第 \(i\) 行，我们需要计算 \(\mathbf{L_{M}w = u}\) 。但是，实际上 \(\mathbf{u}\) 只有一部分 envelope 区域需要计算，因此实际使用和计算的 \(\mathbf{\bar{L}}\) 和 \(\mathbf{\bar{u}^{T}}\) 如右图，因此此时 ELSLV 可以对这个 IBAND 大小的三角方程组进行求解

脚本

具体脚本如下：

现在开始逐行解析。

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

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

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

IFIRST 是第 \(i\) 行第一个元素所在列。

之后调用 ELSLV 求解得到 \(\mathbf{w}\) 。 这里就用到了子矩阵求解的方法，方程数目为 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 ，首先我们从对角线元素中减去 \(\mathbf{w}\) 的内积。

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

	      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.