稀疏矩阵算法Envelope方法一之存储格式

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

Envelope 存储格式

这个存储格式中，数组 ENV 包含矩阵中每一行的 envelope entries ，一个长度为 的索引向量 XENV ，其用于指示每一行的开始位置，而 ，如同之前的索引向量。

中的 元素对应关系如下：

换句话说，矩阵 的 envelope 区域的某个元素 可以表示为 。下图是一个例子，对于其中的 ，我们有 。

更一般的操作是提取 envelope 区域的一行，这个可以按照下面的方式进行

The Storage Allocation Subroutine FNENV

在这一章节中我们讨论子程序 FNENV (FiNd ENVelope) ，这个子程序输入为矩阵 的图，存储在数组对**(XADJ, ADJNCY)** 中，伴随着其置换矩阵 PERM 以及其逆向量 INVP 。假如 ，那么就说明原始编号 在新顺序中其编号是 。我们通常会用一个相应的长度为 的置换向量 INVP (the inverse permutation) ，其满足 。也就是说，INVP(k) 表示原始编号 在新顺序中的位置。

这个子程序的目的是计算上面讨论的数组 XENV ，这可以用于存储 的 Cholesky 分解因子 ；同时返回 ENVSZE，这是 的 envelope size ，等于 ，其中 NEQNS 为节点或者方程的数目。

这里是说， 和 的 envelope 的结构是一样的（证明略）。

这个子程序的内容很简单，循环 DO 200 I= ... 对每一行进行处理，其中循环 DO 100 I= ... 会找到 第 行的第一个非零元素的索引 (IFIRST) 。

下面开始逐行解析脚本，首先初始化 BANDW=0, ENVSZE=1，注意这里的 ENVSZE 其实就是每行的初始位置。

对于每一行进行循环，XENV(I) = ENVSZE 设置每一行的 XENV 的元素，这一行对应的置换矩阵的原始节点编号为 IPERM = PERM(I)，设这一行第一个非零元素位置为 IFIRST ，初始值设为 i

通过 DO 100 遍历其相邻的原始节点，NABOR = ADJNCY(J) 得到其原始编号，NABOR = INVP(NABOR) 进一步得到其在置换矩阵中的新编号，如果新编号低于 IFIRST ，则用这个编号覆盖 IFIRST ，最终得到 IFIRST 的真实值。

通过 I - IFIRST 得到这一行的 bandwidth ，记为 IBAND ，下一行的初始位置更新为 ENVSZE = ENVSZE + IBAND ，如果 BANDW 小于 IBAND ，则 BANDW = IBAND ，最终得到矩阵的 bandwidth 。

	  BANDW = 0
	  ENVSZE = 1
	  DO 200 I=1, NEQNS
	    XENV(I) = ENVSZE
	    IPERM = PERM(I)
	    JSTRT = XADJ(IPERM)
	    JSTOP = XADJ(IPERM + 1) - 1
	    IF (JSTOP .LT. JSTRT) GO TO 200
          IFIRST = I
          DO 100 J = JSTRT, JSTOP
            NABOR = ADJNCY(J)
            NABOR = INVP(NABOR)
            IF ( NABOR .LT. IFIRST ) IFIRST = NABOR
100       CONTINUE
          IBAND = I - IFIRST
          ENVSZE = ENVSZE + IBAND
          IF (BANDW .LT. IBAND) BANDW = IBAND
200   CONTINUE

最后，XENV 需要补上最后一个元素 NEQNS+1

实际的 ENVSZE 需要减1，程序终止

  XENV(NEQNS+1) = ENVSZE
  ENVSZE = ENVSZE - 1
  RETURN
END

参考文献

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