稀疏矩阵算法RCM算法三之RCM算法

本章节是介绍 George and Liu 书中 RCM 算法脚本的第三篇，介绍RCM的完整算法。

程序结构

在这一章节中，我们会介绍三个子程序 DEGREE, RCM, GENRCM ，再加上之前提到的子程序，就构成了 RCM 算法的所有脚本。其输入为 ROOT, MASK, XADJ, ADJNCY 和 PERM 。不同子程序之间的关系见下图：

DEGREE

这个子程序会计算一个连通组分中所有节点的度数。其输入为 ROOT, MASK, XADJ, ADJNCY 。

这个子程序会从第一层 (只包含 ROOT ) 开始，同时计算每一层的所有节点的度数，然后产生下一层的所有节点，加入到 LS 数组中。当一个节点加入到 LS 数组中，其在 XADJ 数组中的符号会改变，从而确保这个节点只会被记录一次（这个功能之前在 ROOTLS 中是通过修改 MASK 数组实现的，但是这里 MASK 必须保持不变，以保证所有节点的度数可以正确计算）。

变量 CCSIZE 包含了当前在 LS 中的节点数目。当所有的节点均被发现之后，所有节点的度数也计算完毕，此时 LS 数组中的节点在 XADJ 中的元素会改回原来的符号。

输出的 LS 的结构是一个节点，后面接着这个节点的所有相邻节点，而不是按照上面的那种层次结构。

脚本如下（IABS()是 fortran 内置的对整数求绝对值的函数）

C****************************************************************          1.
C****************************************************************          2.
C**********    DEGREE ..... DEGREE IN MASKED COMPONENT        ***          3.
C****************************************************************          4.
C****************************************************************          5.
C                                                                          6.
C     PURPOSE - THIS ROUTINE COMPUTES THE DEGREES OF THE NODE              7.
C        IN THE CONNECTED COMPONENT SPECIFIED BY MASK AND ROOT.            8.
C        NODES FOR WHICH MASK IS ZERO ARE IGNORED.                         9.
C                                                                         10.
C    INPUT PARAMETERS -                                                   11.
C        ROOT - IS THE INPUT NODE THAT DEFINES THE COMPONENT              12.
C        (XADJ, ADJNCY) - ADJACENCY STRUCTURE PAIR.                       13.
C        MASK - SPECIFIES A SECTION SUBGRAPH.                             14.
C                                                                         15.
C     OUTPUT PARAMETERS -                                                 16.
C        DEG - ARRAY CONTAINING THE DEGREES OF THE NODES IN               17.
C              THE COMPONENT                                              18.
C        CCSIZE - SIZE OF THE COMPONENT SPECIFIED BY MASK AND ROOT        19.
C                                                                         20.
C     WORKING PARAMETER -                                                 21.
C        LS - A TEMPORARY VECTOR USED TO STORE THE NODES OF THE           22.
C               COMPONENT LEVEL BY LEVEL                                  23.
C                                                                         24.
C****************************************************************         25.
C                                                                         26.
      SUBROUTINE  DEGREE ( ROOT, XADJ, ADJNCY, MASK,                      27.
     1                     DEG, CCSIZE, LS )                              28.
C                                                                         29.
C****************************************************************         30.
C                                                                         31.
         INTEGER ADJNCY(1), DEG(1), LS(1), MASK(1)                        32.
         INTEGER XADJ(1), CCSIZE, I, IDEG, J, JSTOP, JSTRT,               33.
     1           LBEGIN, LVLEND, LVSIZE, NBR, NODE, ROOT                  34.
C                                                                         35.
C****************************************************************         36.
C                                                                         37.
C        -----------------------------------------------------            38.
C        INITIALIZATION ...                                               39.
C        THE ARRAY XADJ IS USED AS A TEMPORARY MARKER TO                  40.
C        INDICATE WHICH NODES HAVE BEEN CONSIDERED SO FAR.                41.
C        -----------------------------------------------------            42.
         LS(1) = ROOT                                                     43.
         XADJ(ROOT) = -XADJ(ROOT)                                         44.
         LVLEND = 0                                                       45.
         CCSIZE = 1                                                       46.
C        -----------------------------------------------------            47.
C        LBEGIN IS THE POINTER OT THE BEGINNING OF THE CURRENT            48.
C        LEVEL. AND LVLEND POINTS TO THE END OF THIS LEVEL.               49.
C        -----------------------------------------------------            50.
   100   LBEGIN = LVLEND + 1                                              51.
         LVLEND = CCSIZE                                                  52.
C        -----------------------------------------------------            53.
C        FIND THE DEGREES OF NODES IN THE CURRENT LEVEL,                  54.
C        AND AT THE SAME TIME, GENERATE THE NEXT LEVEL.                   55.
C        -----------------------------------------------------            56.
         DO 400 I = LBEGIN, LVLEND                                        57.
            NODE = LS(I)                                                  58.
            JSTRT = -XADJ(NODE)                                           59.
            JSTOP = IABS(XADJ(NODE+1)) - 1                                60.
            IDEG = 0                                                      61.
            IF ( JSTOP .LT. JSTRT ) GO TO 300                             62.
               DO 200 J = JSTRT, JSTOP                                    63.
                  NBR = ADJNCY(J)                                         64.
                  IF ( MASK(NBR) .EQ. 0 ) GO TO 200                       65.
                     IDEG = IDEG + 1                                      66.
                     IF ( XADJ(NBR) .LT. 0 ) GO TO 200                    67.
                        XADJ(NBR) = -XADJ(NBR)                            68.
                        CCSIZE = CCSIZE + 1                               69.
                        LS(CCSIZE) = NBR                                  70.
  200          CONTINUE                                                   71.
  300       DEG(NODE) = IDEG                                              72.
  400    CONTINUE                                                         73.
C        -----------------------------------------------------            74.
C        COMPUTE THE CURRENT LEVEL WIDTH.                                 75.
C        IF IT IS NONZERO, GENERATE ANOTHER LEVEL.                        76.
C        -----------------------------------------------------            77.
         LVSIZE = CCSIZE - LVLEND                                         78.
         IF ( LVSIZE .GT. 0 ) GO TO 100                                   79.
C        -----------------------------------------------------            80.
C        RESET XADJ TO ITS CORRECT SIGN AND RETURN.                       81.
C        -----------------------------------------------------            82.
         DO 500 I = 1, CCSIZE                                             83.
            NODE = LS(I)                                                  84.
            XADJ(NODE) = -XADJ(NODE)                                      85.
  500    CONTINUE                                                         86.
         RETURN                                                           87.
      END                                                                 88.

下面我们开始逐行分析脚本。

首先是初始化过程，从 ROOT 节点开始，将其加入到 LS 向量中，并将其 XADJ 中的值改为负数。

凡是加入到 LS 向量中的节点，其 XADJ 中的值均要改为相应的负数，防止重复添加到 LS 向量中。

C        -----------------------------------------------------            38.
C        INITIALIZATION ...                                               39.
C        THE ARRAY XADJ IS USED AS A TEMPORARY MARKER TO                  40.
C        INDICATE WHICH NODES HAVE BEEN CONSIDERED SO FAR.                41.
C        -----------------------------------------------------            42.
         LS(1) = ROOT                                                     43.
         XADJ(ROOT) = -XADJ(ROOT)                                         44.
         LVLEND = 0                                                       45.
         CCSIZE = 1                                                       46.

开始进入循环，其中 LBEGIN 是当前层的起点，LVLEND 是当前层的终点

C        -----------------------------------------------------            47.
C        LBEGIN IS THE POINTER OT THE BEGINNING OF THE CURRENT            48.
C        LEVEL. AND LVLEND POINTS TO THE END OF THIS LEVEL.               49.
C        -----------------------------------------------------            50.
   100   LBEGIN = LVLEND + 1                                              51.
         LVLEND = CCSIZE                                                  52.

查找当前层的所有节点的度数，注意这里已经是加入到 LS 向量中的节点，XADJ(NODE) 本身就是负数，而 XADJ(NODE + 1) 不确定是否是负数，因此采用 IABS() 函数取绝对值。

IDEG 记录度数，其中的节点必须满足 MASK 值不为0的条件；如果这些节点进一步满足其 XADJ 值不为负数，则加入到 LS 向量中，LS 向量长度 CCSIZE 加1，其 XADJ 值改为负数，这些节点作为下一层次的节点。

C        -----------------------------------------------------            53.
C        FIND THE DEGREES OF NODES IN THE CURRENT LEVEL,                  54.
C        AND AT THE SAME TIME, GENERATE THE NEXT LEVEL.                   55.
C        -----------------------------------------------------            56.
         DO 400 I = LBEGIN, LVLEND                                        57.
            NODE = LS(I)                                                  58.
            JSTRT = -XADJ(NODE)                                           59.
            JSTOP = IABS(XADJ(NODE+1)) - 1                                60.
            IDEG = 0                                                      61.
            IF ( JSTOP .LT. JSTRT ) GO TO 300                             62.
               DO 200 J = JSTRT, JSTOP                                    63.
                  NBR = ADJNCY(J)                                         64.
                  IF ( MASK(NBR) .EQ. 0 ) GO TO 200                       65.
                     IDEG = IDEG + 1                                      66.
                     IF ( XADJ(NBR) .LT. 0 ) GO TO 200                    67.
                        XADJ(NBR) = -XADJ(NBR)                            68.
                        CCSIZE = CCSIZE + 1                               69.
                        LS(CCSIZE) = NBR                                  70.
  200          CONTINUE                                                   71.
  300       DEG(NODE) = IDEG                                              72.
  400    CONTINUE                                                         73.

我们计算当前层次的节点数目，如果不为0，则继续循环。

C        -----------------------------------------------------            74.
C        COMPUTE THE CURRENT LEVEL WIDTH.                                 75.
C        IF IT IS NONZERO, GENERATE ANOTHER LEVEL.                        76.
C        -----------------------------------------------------            77.
         LVSIZE = CCSIZE - LVLEND                                         78.
         IF ( LVSIZE .GT. 0 ) GO TO 100                                   79.

最后我们将 XADJ 向量“复原”

C        -----------------------------------------------------            80.
C        RESET XADJ TO ITS CORRECT SIGN AND RETURN.                       81.
C        -----------------------------------------------------            82.
         DO 500 I = 1, CCSIZE                                             83.
            NODE = LS(I)                                                  84.
            XADJ(NODE) = -XADJ(NODE)                                      85.
  500    CONTINUE                                                         86.
         RETURN                                                           87.
      END                                                                 88.

RCM

Reverse Cuthill-McKee, 这个子程序就是对一个连通组分执行 RCM 算法，其输入为 ROOT, MASK, XADJ, ADJNCY 。

因为这个算法要求知道连通组分中所有节点的度数，因此其第一步就是通过调用 DEGREE 子程序计算所有节点的度数，所有节点会按照一层接着一层的顺序排列在 LS 数组中（这里是 PERM 数组），但是每一层的节点之间没有按照度数从低到高排序；因此我们通过 DO 600 I = ... 对每一层的节点按照度数排序或者说编码，生成一个新的顺序的层；通过循环 DO 200 I = ... 查找某个节点的所有未编码的相邻节点，按照度数从低到高的顺序排序。最终新的排序会记录在 PERM 数组中，最后再翻转这个顺序，得到我们最终想要的顺序。

注意，这里就像 ROOTLS 子程序，当节点 \(i\) 编码了之后，MASK(i) 会设为0；但是不像 ROOTLS 子程序， RCM 子程序不会将 MASK 数组恢复如初，其编码后的节点始终为 0 。

C****************************************************************          1.
C****************************************************************          2.
C**********    RCM ..... REVERSE CUTHILL-MCKEE ORDERING       ***          3.
C****************************************************************          4.
C****************************************************************          5.
C                                                                          6.
C     PURPOSE - RCM NUMBERS A CONNECTED COMPONENT SPECIFIED BY             7.
C        MASK AND ROOT, USING THE RCM ALGORITHM.                           8.
C        THE NUMBERING IS TO BE STARTED AT THE NODE ROOT.                  9.
C                                                                         10.
C    INPUT PARAMETERS -                                                   11.
C        ROOT - IS THE NODE THAT DEFINES THE CONNECTED                    12.
C               COMPONENT AND IT IS USED AS THE STARTING                  13.
C               NODE FOR THE RCM ORDERING                                 14.
C        (XADJ, ADJNCY) - ADJACENCY STRUCTURE PAIR FOR                    15.
C                         THE GRAPH.                                      16.
C                                                                         17.
C     UPDATED PARAMETERS -                                                18.
C        MASK - ONLY THOSE NODES WITH NONZERO INPUT MASK                  19.
C               VALUES ARE CONSIDERED BY THE ROUTINE. THE                 20.
C               NODES NUMBERED BY RCM WILL HAVE THEIR                     21.
C               MASK VALUES SET TO ZERO.                                  22.
C                                                                         23.
C     OUTPUT PARAMETERS -                                                 24.
C        PERM - WILL CONTAIN THE RCM ORDERING                             25.
C        CCSIZE - IS THE SIZE OF THE CONNECTED COMPONENT                  26.
C                 THAT HAS BEEN NUMBERED BY RCM.                          27.
C                                                                         28.
C     WORKING PARAMETER -                                                 29.
C        DEG - IS A TEMPORARY VECTOR USED TO HOLD THE DEGREE              30.
C              OF THE NODES IN THE SECTION GRAPH SPECIFIED                31.
C              BY MASK AND ROOT                                           32.
C                                                                         33.
C     PROGRAM SUBROUTINES -                                               34.
C        DEGREE.                                                          35.
C                                                                         36.
C****************************************************************         37.
C                                                                         38.
      SUBROUTINE  RCM ( ROOT, XADJ, ADJNCY, MASK,                         39.
     1                     PERM, CCSIZE, DEG )                            40.
C                                                                         41.
C****************************************************************         42.
C                                                                         43.
         INTEGER ADJNCY(1), DEG(1), MASK(1), PERM(1)                      44.
         INTEGER XADJ(1), CCSIZE, FNBR, I, J, JSTOP,                      45.
     1           JSTRT, K, L, LBEGIN, LNBR, LPERM,                        46.
     1           LVLEND, NBR, NODE, ROOT                                  47.
C                                                                         48.
C****************************************************************         49.
C                                                                         50.
C        -----------------------------------------------------            51.
C        FIND THE DEGREES OF THE NODES IN THE                             52.
C        COMPONENT SPECIFIED BY MASK ADN ROOT.                            53.
C        -----------------------------------------------------            54.
         CALL DEGREE( ROOT, XADJ, ADJNCY, MASK, DEG,                      55.
     1                CCSIZE, PERM )                                      56.
         MASK(ROOT) = 0                                                   57.
         IF ( CCSIZE .LE. 1 ) RETURN                                      58.
         LVLEND = 0                                                       59.
         LNBR = 1                                                         60.
C        -----------------------------------------------------            61.
C        LBEGIN AND LVLEND POINT OT THE BEGINNING AND                     62.
C        THE END OF THE CURRENT LEVEL RESPECTIVELY.                       63.
C        -----------------------------------------------------            64.
  100    LBEGIN = LVLEND + 1                                              65.
         LVLEND = LNBR                                                    66.
         DO 600 I = LBEGIN, LVLEND                                        67.
C           --------------------------------------------------            68.
C           FOR EACH NODE IN CURRENT LEVEL ...                            69.
C           --------------------------------------------------            70.
            NODE = PERM(I)                                                71.
            JSTRT = XADJ(NODE)                                            72.
            JSTOP = XADJ(NODE+1) - 1                                      73.
C           --------------------------------------------------            74.
C           FIND THE UNNUMBERED NEIGHBORS OF NODE.                        75.
C           FNBR AND LNBR POINT TO THE FIRST AND LAST                     76.
C           UNNUMBERED NEIGHBORS RESPECTIVELY OF THE CURRENT              77.
C           NODE IN PERM.                                                 78.
C           --------------------------------------------------            79.
            FNBR = LNBR + 1                                               80.
            DO 200 J = JSTRT, JSTOP                                       81.
               NBR = ADJNCY(J)                                            82.
               IF ( MASK(NBR) .EQ. 0 ) GO TO 200                          83.
                  LNBR = LNBR + 1                                         84.
                  MASK(NBR) = 0                                           85.
                  PERM(LNBR) = NBR                                        86.
  200       CONTINUE                                                      87.
            IF ( FNBR .GE. LNBR ) GO TO 600                               88.
C              -----------------------------------------------            89.
C              SORT THE NEIGHBORS OF NODE IN INCREASING                   90.
C              ORDER BY DEGREE. LINEAR INSERTION IS USED.                 91.
C              -----------------------------------------------            92.
               K = FNBR                                                   93.
  300          L = K                                                      94.
                  K = K + 1                                               95.
                  NBR = PERM(K)                                           96.
  400             IF ( L .LT. FNBR ) GO TO 500                            97.
                     LPERM = PERM(L)                                      98.
                     IF ( DEG(LPERM) .LE. DEG(NBR) ) GO TO 500            99.
                        PERM(L+1) = LPERM                                100.
                        L = L -1                                         101.
                        GO TO 400                                        102.
  500             PERM(L+1) = NBR                                        103.
                  IF ( K .LT. LNBR ) GO TO 300                           104.
  600    CONTINUE                                                        105.
         IF (LNBR .GT. LVLEND) GO TO 100                                 106.
C        -----------------------------------------------------           107.
C        WE NOW HAVE THE CUTHILL MCKEE ORDERING.                         108.
C        REVERSE IT BELOW ...                                            109.
C        -----------------------------------------------------           110.
         K = CCSIZE/2                                                    111.
         L = CCSIZE                                                      112.
         DO 700 I = 1, K                                                 113.
            LPERM = PERM(L)                                              114.
            PERM(L) = PERM(I)                                            115.
            PERM(I) = LPERM                                              116.
            L = L - 1                                                    117.
  700    CONTINUE                                                        118.
         RETURN                                                          119.
      END                                                                120.

下面我们开始逐行查看这个脚本

首先我们调用 DEGREE 子程序，找到连通组分中所有节点的度数，注意这里的 LS 向量替换为 PERM 向量

C        -----------------------------------------------------            51.
C        FIND THE DEGREES OF THE NODES IN THE                             52.
C        COMPONENT SPECIFIED BY MASK ADN ROOT.                            53.
C        -----------------------------------------------------            54.
         CALL DEGREE( ROOT, XADJ, ADJNCY, MASK, DEG,                      55.
     1                CCSIZE, PERM )                                      56.

然后我们从 ROOT 出发，将其 MASK 值设为 0

MASK(ROOT) = 0                                                   57.
IF ( CCSIZE .LE. 1 ) RETURN                                      58.
LVLEND = 0                                                       59.
LNBR = 1                                                         60.

LBEGIN 和 LVLEND 指向当前层次的初始和结束位置，进入循环。

FNBR 和 LNBR 指向这个节点的所有未排序的相邻节点的初始和结束位置（由于是连通子图，每个节点至少有一个相邻节点；但是考虑到 MASK 值为0的情况，因此可能有0个相邻节点）。并且将相应的节点重新赋值给 PERM 向量（这些节点之间还没有排序）。

C        -----------------------------------------------------            61.
C        LBEGIN AND LVLEND POINT OT THE BEGINNING AND                     62.
C        THE END OF THE CURRENT LEVEL RESPECTIVELY.                       63.
C        -----------------------------------------------------            64.
  100    LBEGIN = LVLEND + 1                                              65.
         LVLEND = LNBR                                                    66.
         DO 600 I = LBEGIN, LVLEND                                        67.
C           --------------------------------------------------            68.
C           FOR EACH NODE IN CURRENT LEVEL ...                            69.
C           --------------------------------------------------            70.
            NODE = PERM(I)                                                71.
            JSTRT = XADJ(NODE)                                            72.
            JSTOP = XADJ(NODE+1) - 1                                      73.
C           --------------------------------------------------            74.
C           FIND THE UNNUMBERED NEIGHBORS OF NODE.                        75.
C           FNBR AND LNBR POINT TO THE FIRST AND LAST                     76.
C           UNNUMBERED NEIGHBORS RESPECTIVELY OF THE CURRENT              77.
C           NODE IN PERM.                                                 78.
C           --------------------------------------------------            79.
            FNBR = LNBR + 1                                               80.
            DO 200 J = JSTRT, JSTOP                                       81.
               NBR = ADJNCY(J)                                            82.
               IF ( MASK(NBR) .EQ. 0 ) GO TO 200                          83.
                  LNBR = LNBR + 1                                         84.
                  MASK(NBR) = 0                                           85.
                  PERM(LNBR) = NBR                                        86.
  200       CONTINUE                                                      87.

如果没有MASK 不为0的相邻节点，则 FNBR = LNBR + 1；如果只有一个MASK 值不为0的相邻节点，此时 FNBR = LNBR。这两种情况下，相邻节点不需要排序，GO TO 600，查找下一个节点。

IF ( FNBR .GE. LNBR ) GO TO 600                               88.

正常情况下，我们要对刚刚加入到 PERM 向量中的向量节点进行排序（注意这里还在 DO 600 的循环内部）。

下面这段脚本的可读性太差了。看了半天大概看懂了，第 K 个位置和节点 NBR 相当于新增的节点，而之前的 FNBR 到 K-1 默认是已经排好序的节点。首先将第 K 个位置的节点与 K-1 位置的节点比较，如果左边大于等于右边则不用重排；否则将 K-1 位置的节点放到 K 位置，将 NBR 再和 K-2 位置的节点比较，依次比下去，直到找到 NBR 节点合适的插入位置。

其实就是看新增的 NBR 节点适合插入哪个位置，其实就是插入排序算法。

C              -----------------------------------------------            89.
C              SORT THE NEIGHBORS OF NODE IN INCREASING                   90.
C              ORDER BY DEGREE. LINEAR INSERTION IS USED.                 91.
C              -----------------------------------------------            92.
               K = FNBR                                                   93.
  300          L = K                                                      94.
                  K = K + 1                                               95.
                  NBR = PERM(K)                                           96.
  400             IF ( L .LT. FNBR ) GO TO 500                            97.
                     LPERM = PERM(L)                                      98.
                     IF ( DEG(LPERM) .LE. DEG(NBR) ) GO TO 500            99.
                        PERM(L+1) = LPERM                                100.
                        L = L -1                                         101.
                        GO TO 400                                        102.
  500             PERM(L+1) = NBR                                        103.
                  IF ( K .LT. LNBR ) GO TO 300                           104.
  600    CONTINUE                                                        105.

对这一层的所有节点的相邻节点进行排序，加入到 PERM 中后，判断是否新增了节点，即 LNBR 是否大于 LVLEND ，LNBR 是此时 PERM 的长度，LVLEND 是之前 PERM 的长度。如果新增了节点，则继续循环

IF (LNBR .GT. LVLEND) GO TO 100                                 106.

最后我们对 PERM 向量进行翻转，K 是总数除以2取整，然后第一个数和最后一个数互换，第二个数和倒数第二个数互换，以此类推。

举个例子，假设 CCSIZE = 5，那么 K = 5/2=2，那么就是位置1和位置5的数互换，位置2和位置4的数互换，位置3不动。

C        -----------------------------------------------------           107.
C        WE NOW HAVE THE CUTHILL MCKEE ORDERING.                         108.
C        REVERSE IT BELOW ...                                            109.
C        -----------------------------------------------------           110.
         K = CCSIZE/2                                                    111.
         L = CCSIZE                                                      112.
         DO 700 I = 1, K                                                 113.
            LPERM = PERM(L)                                              114.
            PERM(L) = PERM(I)                                            115.
            PERM(I) = LPERM                                              116.
            L = L - 1                                                    117.
  700    CONTINUE                                                        118.
         RETURN                                                          119.
      END                                                                120.

GENRCM

(GENeral RCM) ，这个子程序用于找到一个不连通的图的 RCM 顺序，这个子程序会对每一个连通组分运行子程序 RCM 。这个子程序的输入为 NEQNS (the number of nodes or equations)，(XADJ, ADJNCY)。中间数据为 MASK 和 XLS ，作为中间子程序 FNROOT 和 RCM 的输入。

这个子程序一开始将 MASK 中的所有值均设为 1，之后对 MASK 进行循环知道其发现一个节点 \(i\) ，其 MASK(i) = 1; 节点 \(i\) 和 MASK, (XADJ, ADJNCY) 便可以指定一个连通子图，之后调用子程序 FNROOT 和 RCM 来确定这个子图中所有节点的顺序（注意运行完 RCM 之后，所有排序号的节点其 MASK 值会设为 0）。

注意 NUM 指向数组 PERM 中的第一个自由的位置 (first free position) ，每次运行完 RCM 子程序后就会更新。GENRCM 程序中运行 RCM 中关于 PERM 的实际参数为 PERM(NUM) ，也就是说在 RCM 中的 PERM(NUM) 参数是 GENRCM 中 PERM 数据的最后 NEQNS - NUM + 1 个元素。

最后当一个连通组分排好序后，GENRCM 程序会查找下一个 \(MASK(i) = 1\) 的节点，对下一个连通组分进行排序，直到终止。

C****************************************************************          1.
C****************************************************************          2.
C**********    GENRCM ..... GENERAL REVERSE CUTHILL MCKEE     ***          3.
C****************************************************************          4.
C****************************************************************          5.
C                                                                          6.
C     PURPOSE - GENRCM FINDS THE REVERSE CUTHILL-MCKEE                     7.
C        ORDERING FOR A GENERAL GRAPH. FOR EACH CONNECTED                  8.
C        COMPONENT IN THE GRAPH, GENRCM OBTAINS THE ORDERING               9.
C        BY CALLING THE SUBROUTINE RCM.                                   10.
C                                                                         11.
C     INPUT PARAMETERS -                                                  12.
C        NEQNS - NUMBER OF EQUATIONS                                      13.
C        (XADJ, ADJNCY) - ARRAY PAIR CONTAINING THE ADJACENCY             14.
C               STRUCTURE OF THE GRAPH OF THE MATRIX                      15.
C                                                                         16.
C     OUTPUT PARAMETER -                                                  17.
C        PERM - VECTOR THAT CONTAINS THE RCM ORDERING.                    18.
C                                                                         19.
C     WORKING PARAMETERS -                                                20.
C        MASK - IS USED TO MARK VARIABLES THAT HAVE BEEN                  21.
C               NUMBERED DURING THE ORDERING PROCESS. IT IS               22.
C               INITIALIZED TO 1, AND SET TO ZERO AS EACH NODE            23.
C               IS NUMBERED.                                              24.
C        XLS - THE INDEX VECTOR FOR A LEVEL STRUCTURE. THE                25.
C               LEVEL STRUCTURE IS STORED IN THE CURRENTLY                26.
C               UNUSED SPACES IN THE PERMUTATION VECTOR PERM.             27.
C                                                                         28.
C     PROGRAM SUBROUTINES -                                               29.
C        FNROOT, RCM.                                                     30.
C                                                                         31.
C****************************************************************         32.
C                                                                         33.
      SUBROUTINE  GENRCM ( NEQNS, XADJ, ADJNCY, PERM, MASK, XLS )         34.
C                                                                         35.
C****************************************************************         36.
C                                                                         37.
         INTEGER ADJNCY(1), MASK(1), PERM(1), XLS(1)                      38.
         INTEGER XADJ(1), CCSIZE, I, NEQNS, NLVL                          39.
     1           NUM, ROOT                                                40.
C                                                                         41.
C****************************************************************         42.
C                                                                         43.
         DO 100 I = 1, NEQNS                                              44.
            MASK(I) = 1                                                   45.
  100    CONTINUE                                                         46.
         NUM = 1                                                          47.
         DO 200 I = 1, NEQNS                                              48.
C           -----------------------------------------------------         49.
C           FOR EACH MASKED CONNECTED COMPONENT ...                       50.
C           -----------------------------------------------------         51.
            IF (MASK(I) .EQ. 0) GO TO 200                                 52.
               ROOT = I                                                   53.
C              --------------------------------------------------         54.
C              FIRST FIND A PSEUDO-PERIPHERAL NODE ROOT.                  55.
C              NOTE THAT THE LEVEL STRUCTURE FOUND BY                     56.
C              FNROOT IS STORED STARTING AT PERM(NUM).                    57.
C              THEN RCM IS CALLED TO ORDER THE COMPONENT                  58.
C              USING ROOT AS THE STARTING NODE.                           59.
C              --------------------------------------------------         60.
               CALL FNROOT ( ROOT, XADJ, ADJNCY, MASK,                    61.
     1                       NLVL, XLS, PERM(NUM) )                       62.
               CALL    RCM ( ROOT, XADJ, ADJNCY, MASK,                    63.
     1                       PERM(NUM), CCSIZE, XLS )                     64.
               NUM = NUM + CCSIZE                                         65.
               IF (NUM .GT. NEQNS) RETURN                                 66.
  200    CONTINUE                                                         67.
         RETURN                                                           68.
      END                                                                 69.

下面我们开始逐行解析脚本

首先将 MASK 向量的元素均初始化为 1，NUM 初始设为 1

       DO 100 I = 1, NEQNS                                              44.
          MASK(I) = 1                                                   45.
100    CONTINUE                                                         46.
       NUM = 1                                                          47.

开始循环，找到第一个 MASK 值不为 0 个节点，设为 ROOT

然后调用 FNROOT 查找一个合适的初始节点，然后调用 RCM 算法查找合适的顺序，其中 NUM 表示 PERM 第一额可用的位置。

FNROOT 子程序中的 LS 向量用 PERM(NUM) 替代，RCM 子程序中的 DEG 向量用 XLS 替代，应该是出于节省内存的作用。这两个都只是中间向量，没有必要再创建两个数组。

调用完这两个子程序之后，已经排序的节点的 MASK 值已经设为 0 了。

之后对 NUM 进行更新，如果 NUM 大于 NEQNS ，说明所有节点已经排序，退出程序 (RETURN) ；否则继续对下一个连通组分进行排序，直至所有节点均被排序。

         DO 200 I = 1, NEQNS                                              48.
C           -----------------------------------------------------         49.
C           FOR EACH MASKED CONNECTED COMPONENT ...                       50.
C           -----------------------------------------------------         51.
            IF (MASK(I) .EQ. 0) GO TO 200                                 52.
               ROOT = I                                                   53.
C              --------------------------------------------------         54.
C              FIRST FIND A PSEUDO-PERIPHERAL NODE ROOT.                  55.
C              NOTE THAT THE LEVEL STRUCTURE FOUND BY                     56.
C              FNROOT IS STORED STARTING AT PERM(NUM).                    57.
C              THEN RCM IS CALLED TO ORDER THE COMPONENT                  58.
C              USING ROOT AS THE STARTING NODE.                           59.
C              --------------------------------------------------         60.
               CALL FNROOT ( ROOT, XADJ, ADJNCY, MASK,                    61.
     1                       NLVL, XLS, PERM(NUM) )                       62.
               CALL    RCM ( ROOT, XADJ, ADJNCY, MASK,                    63.
     1                       PERM(NUM), CCSIZE, XLS )                     64.
               NUM = NUM + CCSIZE                                         65.
               IF (NUM .GT. NEQNS) RETURN                                 66.
  200    CONTINUE                                                         67.
         RETURN                                                           68.
      END                                                                 69.

参考文献

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