稀疏矩阵算法RCM算法二之查找初始位点程序

理论学习 / 稀疏矩阵
字数统计: 3k   |   阅读时长≈ 14 分钟

本章节是介绍 George and Liu 书中 RCM 算法脚本的第二篇，介绍查找初始节点的算法。

ROOTLS

(ROOTed Level Structure) 这个子程序的目的是产生一个连通组分的 level structure 。

其输入参数为 ROOT, MASK, XADJ, ADJNCY。(XADJ, ADJNCY) 为整个图的相邻结构， MASK 和 ROOT 指定了连通子图，MASK(i)=0 的节点会直接忽略，ROOT 为起始节点。

其输出的以ROOT为根的 level structure ，保存在数组对 (XLS, LS) 中，其第 层的节点为 。总层数存储在 NLVL 。

注意由于 Fortran 是 1-index ，因此算法中 0 水平在代码中是 1，因此这里的 其实是上面算法里的 层 ，而 。

这个子程序会逐层查找节点，所有找到的节点会存储在数组 LS 中，然后其相应的 MASK 的值设置为 0 ，因此这些节点不会再出现在 LS 中。当所有的水平都找完，MASK 中所有 level structure 中的节点的值会重新设置为 1 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
C****************************************************************          1.
C****************************************************************          2.
C**********    ROOTLS ..... ROOTED LEVEL STRUCTURE            ***          3.
C****************************************************************          4.
C****************************************************************          5.
C                                                                          6.
C     PURPOSE - ROOTLS GENERATES THE LEVEL STRUCTURE ROOTED                7.
C        AT THE INPUT NODE CALLED ROOT. ONLY THOSE NODES FOR               8.
C        WHICH MASK IS NONZERO WILL BE CONSIDERED.                         9.
C                                                                         10.
C    INPUT PARAMETERS -                                                   11.
C        ROOT - THE NODE AT WHICH THE LEVEL STRUCTURE IS TO               12.
C               BE ROOTED.                                                13.
C        (XADJ, ADJNCY) - ADJACENCY STRUCTURE PAIR FOR THE                14.
C                         GIVEN GRAPH.                                    15.
C        MASK - IS USED TO SPECIFY A SECTION SUBGRAPH. NODES              16.
C               WITH MASK(I)=0 ARE IGNORED                                17.
C                                                                         18.
C     OUTPUT PARAMETERS -                                                 19.
C        NLVL - IS THE NUMBER OF LEVELS IN THE LEVEL STRUCTURE.           20.
C        (XLS, LS) - ARRAY PAIR FOR THE ROOTED LEVEL STRUCTURE.           21.
C                                                                         22.
C****************************************************************         23.
C                                                                         24.
      SUBROUTINE  ROOTLS ( ROOT, XADJ, ADJNCY, MASK, NLVL, XLS, LS )      25.
C                                                                         26.
C****************************************************************         27.
C                                                                         28.
         INTEGER ADJNCY(1), LS(1), MASK(1), XLS(1)                        29.
         INTEGER XADJ(1), I, J, JSTOP, JSTRT, LBEGIN,                     30.
     1           CCSIZE, LVLEND, LVSIZE, NBR, NLVL,                       31.
     1           NODE, ROOT                                               32.
C                                                                         33.
C****************************************************************         34.
C                                                                         35.
C        -----------------------------------------------------            36.
C        INITIALIZATION ...                                               37.
C        -----------------------------------------------------            38.
         MASK(ROOT) = 0                                                   39.
         LS(1) = ROOT                                                     40.
         NLVL = 0                                                         41.
         LVLEND = 0                                                       42.
         CCSIZE = 1                                                       43.
C        -----------------------------------------------------            44.
C        LBEGIN IS THE POINTER TO THE BEGINNING OF THE CURRENT            45.
C        LEVEL, AND LVLEND POINTS TO THE END OF THIS LEVEL.               46.
C        -----------------------------------------------------            47.
  200    LBEGIN = LVLEND + 1                                              48.
         LVLEND = CCSIZE                                                  49.
         NLVL = NLVL + 1                                                  50.
         XLS(NLVL) = LBEGIN                                               51.
C        -----------------------------------------------------            52.
C        GENERATE THE NEXT LEVEL BY FINDING ALL THE MASKED                53.
C        NEIGHBORS OF NODES IN THE CURRENT LEVEL.                         54.
C        -----------------------------------------------------            55.
         DO 400 I = LBEGIN, LVLEND                                        56.
            NODE = LS(I)                                                  57.
            JSTRT = XADJ(NODE)                                            58.
            JSTOP = XADJ(NODE+1) - 1                                      59.
            IF (JSTOP .LT. JSTRT) GO TO 400                               60.
               DO 300 J = JSTRT, JSTOP                                    61.
                  NBR = ADJNCY(J)                                         62.
                  IF (MASK(NBR) .EQ. 0) GO TO 300                         63.
                     CCSIZE = CCSIZE + 1                                  64.
                     LS(CCSIZE) = NBR                                     65.
                     MASK(NBR) = 0                                        66.
  300          CONTINUE                                                   67.
  400    CONTINUE                                                         68.
C        -----------------------------------------------------            69.
C        COMPUTE THE CURRENT WIDTH.                                       70.
C        IF IT IS NONZERO, GENERATE THE NEXT LEVEL.                       71.
C        -----------------------------------------------------            72.
         LVSIZE = CCSIZE - LVLEND                                         73.
         IF (LVSIZE .GT. 0) GO TO 200                                     74.
C        -----------------------------------------------------            75.
C        RESET MASK TO ONE FOR THE NODES IN THE LEVEL STRUCTURE.          76.
C        -----------------------------------------------------            77.
         XLS(NLVL+1) = LVLEND + 1                                         78.
         DO 500 I = 1, CCSIZE                                             79.
            NODE = LS(I)                                                  80.
            MASK(NODE) = 1                                                81.
  500    CONTINUE                                                         82.
         RETURN                                                           83.
      END                                                                 84.

下面逐行进行解析，根据注释，首先是一个初始化过程，第一层就是 ROOT 节点本身，加入到 LS 中，并将其 MASK 值设为0。总层数 NLVL 初始化为 0，LVLEND 设为0 (这个变量的含义下面会提到)，CCSIZE 设为 1 (CCSIZE 为当前 LS 向量的长度)

1
2
3
4
5
6
7
8
C        -----------------------------------------------------            36.
C        INITIALIZATION ...                                               37.
C        -----------------------------------------------------            38.
         MASK(ROOT) = 0                                                   39.
         LS(1) = ROOT                                                     40.
         NLVL = 0                                                         41.
         LVLEND = 0                                                       42.
         CCSIZE = 1                                                       43.

下一步，我们获得当前层次在 LS 的初始位置 LBEGIN 和结束位置 LVLEND，由于第一层只有 ROOT 这个节点，因此目前 LBEGIN = LVLEND = 1。总层数 NLVL 加一，将本层次的初始位置加入到 XLS 中。

1
2
3
4
5
6
7
8
C        -----------------------------------------------------            44.
C        LBEGIN IS THE POINTER TO THE BEGINNING OF THE CURRENT            45.
C        LEVEL, AND LVLEND POINTS TO THE END OF THIS LEVEL.               46.
C        -----------------------------------------------------            47.
  200    LBEGIN = LVLEND + 1                                              48.
         LVLEND = CCSIZE                                                  49.
         NLVL = NLVL + 1                                                  50.
         XLS(NLVL) = LBEGIN                                               51.

对本层次的所有节点逐一查找相邻节点，其中如果某一个节点的相邻节点数目为 0 ，那么存在 JSTOP < JSTRT ，此时直接跳过这个节点（GO TO 400）。

对于下一层的向量节点，其必须满足 MASK 值不为 0 的条件，这样做有两个作用，一是找到的节点是我们需要的连通组分内的节点，二是避免 LS 向量中出现重复节点。对于满足 MASK 值不为 0 的相邻节点，CCSIZE 加1，将其加入到 LS 向量中，其 MASK 值设为 0。

因此，运行完这一步之后，LS 向量中加入了下一层次的节点，其长度 CCSIZE 同样更新为当前 LS 向量的长度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
C        -----------------------------------------------------            52.
C        GENERATE THE NEXT LEVEL BY FINDING ALL THE MASKED                53.
C        NEIGHBORS OF NODES IN THE CURRENT LEVEL.                         54.
C        -----------------------------------------------------            55.
         DO 400 I = LBEGIN, LVLEND                                        56.
            NODE = LS(I)                                                  57.
            JSTRT = XADJ(NODE)                                            58.
            JSTOP = XADJ(NODE+1) - 1                                      59.
            IF (JSTOP .LT. JSTRT) GO TO 400                               60.
               DO 300 J = JSTRT, JSTOP                                    61.
                  NBR = ADJNCY(J)                                         62.
                  IF (MASK(NBR) .EQ. 0) GO TO 300                         63.
                     CCSIZE = CCSIZE + 1                                  64.
                     LS(CCSIZE) = NBR                                     65.
                     MASK(NBR) = 0                                        66.
  300          CONTINUE                                                   67.
  400    CONTINUE                                                         68.

下面我们查看添加的节点数目是否为 0，如果为 0，那就说明找到头了；不为0，则 GO TO 200 ，接着循环，查找下一层次的节点。

这里 LVLEND 是上一次循环中 LS 向量的长度，CCSIZE 是当前 LS 向量的长度，二者一相减就是本次循环新增的节点数目。

1
2
3
4
5
6
C        -----------------------------------------------------            69.
C        COMPUTE THE CURRENT WIDTH.                                       70.
C        IF IT IS NONZERO, GENERATE THE NEXT LEVEL.                       71.
C        -----------------------------------------------------            72.
         LVSIZE = CCSIZE - LVLEND                                         73.
         IF (LVSIZE .GT. 0) GO TO 200                                     74.

因此，当循环结束后，我们就找到了 level structure 中的所有节点，因此此时 LS 向量就是最终的 LS 向量。

但是 XLS 向量还缺少一个元素，就是 LS 长度 + 1，即下面的 LVLEND + 1 (此时 LVLEND = CCSIZE, 也是 LS 向量长度)。

1
2
3
4
C        -----------------------------------------------------            75.
C        RESET MASK TO ONE FOR THE NODES IN THE LEVEL STRUCTURE.          76.
C        -----------------------------------------------------            77.
         XLS(NLVL+1) = LVLEND + 1                                         78.

最后，我们将 LS 向量中的所有节点的 MASK 重新设为 1，也就是说 “还原” MASK 向量。

1
2
3
4
5
6
       DO 500 I = 1, CCSIZE                                             79.
          NODE = LS(I)                                                  80.
          MASK(NODE) = 1                                                81.
500    CONTINUE                                                         82.
       RETURN                                                           83.
    END                                                                 84.

FNROOT

(FiNd ROOT) 这个子程序用于查找一个连通组分的 pseudo-peripheral node 。

其输入与上一个子程序相同，首先调用 ROOTLS 这个子程序，如果这个连通组分只有 ROOT 这一个节点或者只包含以 ROOT 为终点的一条链，那么 ROOT 就是要找的 pseudo-peripheral node ， LS 就是其相应的 rooted level structure ，程序终止 。

否则，程序会找到最后一个水平的度数最小的一个节点，针对这个节点生成其相应的新的 level structure ，我们会检查是否 ，如果不满足则程序终止，否则则重复上述流程。

当程序终止的时候，ROOT 就是我们想要找的 pseudo-peripheral node ，而数组对 (XLS, LS) 包含了相应的 rooted level structure 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
C****************************************************************          1.
C****************************************************************          2.
C**********    FNROOT ..... FIND PSEUDO-PERIPHERAL NODE       ***          3.
C****************************************************************          4.
C****************************************************************          5.
C                                                                          6.
C     PURPOSE - FNROOT IMPLEMENTS A MODIFIED VERSION OF THE                7.
C        SCHEME BY GIBBS, POOLE, AND STOCKMEYER TO FIND PSEUDO-            8.
C        PERIPHERAL NODES. IT DETERMINES SUCH A NODE FOR THE               9.
C        SECTION SUBGRAPH SPECIFIED BY MASK AND ROOT                      10.
C                                                                         11.
C    INPUT PARAMETERS -                                                   12.
C        (XADJ, ADJNCY) - ADJACENCY STRUCTURE PAIR FOR THE GRAPH.         13.
C        MASK - SPECIFIES A SECTION SUBGRAPH. NODES FOR WHICH             14.
C               MASK IS ZERO ARE IGNORED BY FNROOT.                       15.
C                                                                         16.
C     UPDATED PARAMETERS -                                                17.
C        ROOT - ON INPUT, IT (ALONG WITH MASK) DEFINES THE                18.
C               COMPONENT FOR WHICH A PSEUDO-PERIPHERAL NODE IS           19.
C               TO BE FOUND. ON OUTPUT, IT IS THE NODE OBTAINED           20.
C                                                                         21.
C     OUTPUT PARAMETERS -                                                 22.
C        NLVL - IS THE NUMBER OF LEVELS IN THE LEVEL STRUCTURE            23.
C               ROOTED AT THE NODE ROOT.                                  24.
C        (XLS, LS) - THE LEVEL STRUCTURE ARRAY PAIR CONTAINING            25.
C                    THE LEVEL STRUCTURE FOUND.                           26.
C                                                                         27.
C     PROGRAM SUBROUTINES -                                               28.
C        ROOTLS.                                                          29.
C                                                                         30.
C****************************************************************         31.
C                                                                         32.
      SUBROUTINE  FNROOT ( ROOT, XADJ, ADJNCY, MASK, NLVL, XLS, LS )      33.
C                                                                         34.
C****************************************************************         35.
C                                                                         36.
         INTEGER ADJNCY(1), LS(1), MASK(1), XLS(1)                        37.
         INTEGER XADJ(1), CCSIZE, J, JSTRT, K, KSTOP, KSTRT,              38.
     1           MINDEG, NABOR, NDEG, NLVL,NODE, NUNLVL,                  39.
     1           ROOT                                                     40.
C                                                                         41.
C****************************************************************         42.
C                                                                         43.
C        -----------------------------------------------------            44.
C        DETERMINE THE LEVEL STRUCTURE ROOTED AT ROOT.                    45.
C        -----------------------------------------------------            46.
         CALL ROOTLS( ROOT, XADJ, ADJNCY, MASK, NLVL, XLS, LS )           47.
         CCSIZE = XLS(NLVL+1) - 1                                         48.
         IF ( NLVL.EQ. 1 .OR. NLVL .EQ. CCSIZE ) RETURN                   49.
C        -----------------------------------------------------            50.
C        PICK A NODE WITH MINIMUM DEGREE FROM THE LAST LEVEL.             51.
C        -----------------------------------------------------            52.
  100    JSTRT = XLS(NLVL)                                                53.
         MINDEG = CCSIZE                                                  54.
         ROOT = LS(JSTRT)                                                 55.
         IF ( CCSIZE .EQ. JSTRT) GO TO 400                                56.
            DO 300 J = JSTRT, CCSIZE                                      57.
               NODE = LS(J)                                               58.
               NDEG = 0                                                   59.
               KSTRT = XADJ(NODE)                                         60.
               KSTOP = XADJ(NODE+1) - 1                                   61.
               DO 200 K = KSTRT, KSTOP                                    62.
                  NABOR = ADJNCY(K)                                       63.
                  IF ( MASK(NABOR) .GT. 0 ) NDEG = NDEG + 1               64.
  200          CONTINUE                                                   65.
               IF ( NDEG .GE. MINDEG ) GO TO 300                          66.
                  ROOT = NODE                                             67.
                  MINDEG = NDEG                                           68.
  300       CONTINUE                                                      69.
C        -----------------------------------------------------            70.
C        AND GENERATE ITS ROOTED LEVEL STRUCTURE                          71.
C        -----------------------------------------------------            72.
  400    CALL ROOTLS( ROOT, XADJ, ADJNCY, MASK, NUNLVL, XLS, LS )         73.
         IF (NUNLVL .LE. NLVL) RETURN                                     74.
            NLVL = NUNLVL                                                 75.
            IF ( NLVL .LT. CCSIZE ) GO TO 100                             76.
            RETURN                                                        77.
      END                                                                 78.

下面我们进行逐句解析

第一步，调用 ROOTLS 找到 ROOT 的 level structure，然后同样将 CCSIZE 设为 level structure 中节点数目。注意，由于我们是在一个连通组分中查找节点，因此所有可用的节点互相都是连通的，CCSIZE 就是这个连通组分的所有节点的数目，是一个固定值，也是所有可能的节点集合的节点数目的最大值。

然后我们做一个判断，如果总层数为1，或者总层数等于总节点数，则直接退出程序，输入节点 ROOT 就是我们要找的节点。总节点数为1，说明连通组分中只有 ROOT 这一个节点；总层数等于总节点数，说明每一层只有一个节点，也就是连通组分中只有一条链，而 ROOT 是这条链的一个起点/终点。因此，这两种情况下，输入节点 ROOT 就是我们要找的节点。

1
2
3
4
5
6
C        -----------------------------------------------------            44.
C        DETERMINE THE LEVEL STRUCTURE ROOTED AT ROOT.                    45.
C        -----------------------------------------------------            46.
         CALL ROOTLS( ROOT, XADJ, ADJNCY, MASK, NLVL, XLS, LS )           47.
         CCSIZE = XLS(NLVL+1) - 1                                         48.
         IF ( NLVL.EQ. 1 .OR. NLVL .EQ. CCSIZE ) RETURN                   49.

如果不是上面的两种特殊情况，那么我们接着往下走。

下面这一段脚本的意思，就是从 ROOT 的 level structure 的最后一层中，找到一个度数最小的节点，然后覆盖 ROOT 。

JSTRT 是最后一层的初始节点位置。采用一个中间变量 MINDEG，存储最小的度数，其初始值设置为可能的最大值 CCSIZE 。ROOT 首先设为最后一层的第一个节点。

如果 CCSIZE 等于 JSTRT ，那就说明最后一层只有一个节点，此时 ROOT 就是我们要找的节点，直接 GO TO 400 ，跳过下面这一段。

不然地话，采用 DO 300 循环便利最后一层的所有节点，对于每一个节点采用 DO 200 循环计算该节点的度数，存储在 NDEG 中。

如果 NDEG 小于 MINDEG ，则将 ROOT 赋值为 NOTE ，MINDEG 赋值为 NDEG 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
C        -----------------------------------------------------            50.
C        PICK A NODE WITH MINIMUM DEGREE FROM THE LAST LEVEL.             51.
C        -----------------------------------------------------            52.
  100    JSTRT = XLS(NLVL)                                                53.
         MINDEG = CCSIZE                                                  54.
         ROOT = LS(JSTRT)                                                 55.
         IF ( CCSIZE .EQ. JSTRT) GO TO 400                                56.
            DO 300 J = JSTRT, CCSIZE                                      57.
               NODE = LS(J)                                               58.
               NDEG = 0                                                   59.
               KSTRT = XADJ(NODE)                                         60.
               KSTOP = XADJ(NODE+1) - 1                                   61.
               DO 200 K = KSTRT, KSTOP                                    62.
                  NABOR = ADJNCY(K)                                       63.
                  IF ( MASK(NABOR) .GT. 0 ) NDEG = NDEG + 1               64.
  200          CONTINUE                                                   65.
               IF ( NDEG .GE. MINDEG ) GO TO 300                          66.
                  ROOT = NODE                                             67.
                  MINDEG = NDEG                                           68.
  300       CONTINUE                                                      69.

下面，我们对找到的最后一层度数最小的节点 ROOT ，再运行一遍 ROOTLS 子程序，唯一区别是总层数存储在变量 NUNLVL 中。

首先 NUNLVL 不可能小于 NLVL ，因为按照 level structure 的结构，NLVL - 1 就是初始节点到最后一层节点的最短路径长度， 因此反过来从最后一层的最小度数节点出发，其 level structure 的最后一层如果包含了初始节点，那么 NUNLVL = NLVL ；如果不包含初始节点，那么 NUNLVL > NLVL ，因此不存在 NUNLVL < NLVL 的可能。

因此这里的 IF (NUNLVL .LE. NLVL) RETURN 等价于IF (NUNLVL .EQ. NLVL) RETURN，也就是说如果 NUNLVL 等于 NLVL ，子程序终止 (RETURN) 。

否则的话，将 NUNLVL 的值赋值给 NLVL，判断 NLVL 是否小于 CCSIZE ，如果 NLVL小于 CCSIZE， GOTO 100，重新进入循环；如果 NLVL等于 CCSIZE (不可能大于)，那么说明总层数等于总位点数，说明连通组分中只有这一条链，而 ROOT 是这条链的起点/终点，此时 ROOT 就是要找的结果，程序终止。

1
2
3
4
5
6
7
8
9
C        -----------------------------------------------------            70.
C        AND GENERATE ITS ROOTED LEVEL STRUCTURE                          71.
C        -----------------------------------------------------            72.
  400    CALL ROOTLS( ROOT, XADJ, ADJNCY, MASK, NUNLVL, XLS, LS )         73.
         IF (NUNLVL .LE. NLVL) RETURN                                     74.
            NLVL = NUNLVL                                                 75.
            IF ( NLVL .LT. CCSIZE ) GO TO 100                             76.
            RETURN                                                        77.
      END                                                                 78.

参考文献

  1. George A, Liu J, Ng E. Computer solution of sparse linear systems[J]. Oak Ridge National Laboratory, 1994.
  • 版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！
  • Copyrights © 2019-2024 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信