稀疏矩阵算法最小度数算法四之最小度数算法主程序

本章节介绍 the minimum degree algorithm ,我个人将其翻译为最小度数算法。这里我们看最小度数算法的主程序部分。

最小度数算法的实施

在最小度数算法的实施过程中,不可区分的节点会视为一个 supernode 。注意此时我们有两种不同情况下的 supernode,一个是消元过程中已经消元的连接起来的节点组成的 supernode ,一个是上面提到的不可区分的节点组成的 supernode

下图就是对于这两种 supernode 的说明,它是对图 5.4.5 中的图是如何通过两种商图存储的,其中的阴影的双层圆圈表示已经消元的节点组成的supernodes ,而空白的双层圆圈表示不可区分的节点组成的 supernodes

1

我们现在会讨论如何实施这个算法的脚本,其输入为图的结构 (XADJ, ADJNCY) 和方程数目 NEQNS ,输出为 PERMINVP

这个子程序要求一些临时向量来实施商图模型和不可区分的节点。当前消元图中所有节点的度数存储在数组 DEG 中,其中已经被消除的节点的 DEG 值会设为 -1 。

在商图模型中,连通的已经消除的节点会合并为一个 supernode ,此时我们需要挑出一个节点作为 supernode 的代表。如果 就是这样一个连通组分,我们通常会选择 中最后消除的节点作为 的代表,其它节点不再使用。

同样地,对于未消除的不可区分的节点,也会挑出一个节点作为代表,其它节点不再使用。

中间向量 MARKER 用于标记在相邻结构中可以忽略的节点,这样的节点的 MARKER 值会设为 -1 ,这个向量也会临时用于生成 reachable setsQSIZEQLINK 数组用于指定不可区分的节点,如果节点 作为代表,那么在这个 supernode 中的节点数目为 QSIZE(i) ,这些节点为

下图说明了 QSIZE, QLINKMARKER 的使用,节点 形成了一个由不可区分的节点组成的 supernode ,其代表为节点 2 。因此节点 5 和 8 的 MARKER 值设置为 -1;节点 形成了一个已经消元的 supernode ,其代表为节点 9 ,因此节点 3 和 6 的 MARKER 值设置为 -1

1

这里总共有 5 个子程序,分别为 GENQMD, QMDRCH, QMDQT, QMDUPD, QMDMRG ,它们的控制关系见下图。

1

GENQMD

(General Quotient Minimum Degree algorithm),这个子程序的目的是对一个不连通的图找到其最小度数排序。其输入为图的结构 (XADJ, ADJNCY) 和方程数目 NEQNS ,输出为 PERMINVP 。运算过程中, (XADJ, ADJNCY) 会被修改,用于存储商图的结构。

这个子程序一开始会初始化中间数组 QSIZE, QLINK, MARKERDEG 向量,之后进行到主循环中。在主循环中,子程序首先会通过 threshold searching 找到一个度数最小的节点,其中有两个变量 THRESHMINDEG,任何当前的度数等于 THRESH 的节点均是度数最小的节点。变量 MINDEG 保存大于 THRESH 的值的最小的度数,并用于更新 THRESH 的值。

当找到了一个最小度数的节点 NODE ,之后通过调用子程序 QMDRCH 来通过已经消除的节点获得 NODE 的 reachable set 。这个集合会存储在向量 RCHSET 中,其大小为 RCHSZE。通过 QLINK 获得 NODE 的不可区分的节点,并且进行编号和消除。

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
C----- SUBROUTINE GENQMD
C**************************************************************** 1.
C**************************************************************** 2.
C********** GENQMD ..... QUOT MIN DEGREE ORDERING ********* 3.
C**************************************************************** 4.
C**************************************************************** 5.
C 6.
C PURPOSE - THIS ROUTINE IMPLEMENTS THE MINIMUM DEGREE 7.
C ALGORITHM. IT MAKES USE OF THE IMPLICIT REPRESENT- 8.
C ATION OF THE ELIMINATION GRAPHS BY QUOTIENT GRAPHS, 9.
C AND THE NOTION OF INDISTINGUISHABLE NODES. 10.
C CAUTION - THE ADJACENCY VECTOR ADJNCY WILL BE 11.
C DESTROYED. 12.
C 13.
C INPUT PARAMETERS - 14.
C NEQNS - NUMBER OF EQUATIONS. 15.
C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 16.
C 17.
C OUTPUT PARAMETERS - 18.
C PERM - THE MINIMUM DEGREE ORDERING. 19.
C INVP - THE INVERSE OF PERM. 20.
C 21.
C WORKING PARAMETERS - 22.
C DEG - THE DEGREE VECTOR. DEG(I) IS NEGATIVE MEANS 23.
C NODE I HAS BEEN NUMBERED. 24.
C MARKER - A MARKER VECTOR, WHERE MARKER(I) IS 25.
C NEGATIVE MEANS NODE I HAS BEEN MERGED WITH 26.
C ANOTHER NODE AND THUS CAN BE IGNORED. 27.
C RCHSET - VECTOR USED FOR THE REACHABLE SET. 28.
C NBRHD - VECTOR USED FOR THE NEIGHBORHOOD SET. 29.
C QSIZE - VECTOR USED TO STORE THE SIZE OF 30.
C INDISTINGUISHABLE SUPERNODES. 31.
C QLINK - VECTOR TO STORE INDISTINGUISHABLE NODES, 32.
C I, QLINK(I), QLINK(QLINK(I)) ... ARE THE 33.
C MEMBERS OF THE SUPERNODE REPRESENTED BY I. 34.
C 35.
C PROGRAM SUBROUTINES - 36.
C QMDRCH, QMDQT, QMDUPD. 37.
C 38.
C**************************************************************** 39.
C 40.
C 41.
SUBROUTINE GENQMD ( NEQNS, XADJ, ADJNCY, PERM, INVP, DEG, 42.
1 MARKER, RCHSET, NBRHD, QSIZE, QLINK, 43.
1 NOFSUB ) 44.
C 45.
C**************************************************************** 46.
C 47.
INTEGER ADJNCY(1), PERM(1), INVP(1), DEG(1), MARKER(1), 48.
1 RCHSET(1), NBRHD(1), QSIZE(1), QLINK(1) 49.
INTEGER XADJ(1), INODE, IP, IRCH, J, MINDEG, NDEG, 50.
1 NEQNS, NHDSZE, NODE, NOFSUB, NP, NUM, NUMP1, 51.
1 NXNODE, RCHSZE, SEARCH, THRESH 52.
C 53.
C**************************************************************** 54.
C 55.
C ----------------------------------------------------- 56.
C INITIALIZE DEGREE VECTOR AND OTHER WORKING VARIABLES. 57.
C ----------------------------------------------------- 58.
MINDEG = NEQNS 59.
NOFSUB = 0 60.
DO 100 NODE = 1, NEQNS 61.
PERM(NODE) = NODE 62.
INVP(NODE) = NODE 63.
MARKER(NODE) = 0 64.
QSIZE(NODE) = 1 65.
QLINK(NODE) = 0 66.
NDEG = XADJ(NODE+1) - XADJ(NODE) 67.
DEG(NODE) = NDEG 68.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 69.
100 CONTINUE 70.
NUM = 0 71.
C ----------------------------------------------------- 72.
C PERFORM THRESHOLD SEARCH TO GET A NODE OF MIN DEGREE. 73.
C VARIABLE SEARCH POINTS TO WHERE SEARCH SHOULD START. 74.
C ----------------------------------------------------- 75.
200 SEARCH = 1 76.
THRESH = MINDEG 77.
MINDEG = NEQNS 78.
300 NUMP1 = NUM + 1 79.
IF ( NUMP1 .GT. SEARCH ) SEARCH = NUMP1 80.
DO 400 J = SEARCH, NEQNS 81.
NODE = PERM(J) 82.
IF ( MARKER(NODE) .LT. 0 ) GOTO 400 83.
NDEG = DEG(NODE) 84.
IF ( NDEG .LE. THRESH ) GO TO 500 85.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 86.
400 CONTINUE 87.
GO TO 200 88.
C --------------------------------------------------- 89.
C NODE HAS MINIMUM DEGREE. FIND ITS REACHABLE SETS BY 90.
C CALLING QMDRCH. 91.
C --------------------------------------------------- 92.
500 SEARCH = J 93.
NOFSUB = NOFSUB + DEG(NODE) 94.
MARKER(NODE) = 1 95.
CALL QMDRCH (NODE, XADJ, ADJNCY, DEG, MARKER, 96.
1 RCHSZE, RCHSET, NHDSZE, NBRHD ) 97.
C ------------------------------------------------ 98.
C ELIMINATE ALL NODES INDISTINGUISHABLE FROM NODE. 99.
C THEY ARE GIVEN BY NODE, QLINK(NODE), .... 100.
C ------------------------------------------------ 101.
NXNODE = NODE 102.
600 NUM = NUM + 1 103.
NP = INVP(NXNODE) 104.
IP = PERM(NUM) 105.
PERM(NP) = IP 106.
INVP(IP) = NP 107.
PERM(NUM) = NXNODE 108.
INVP(NXNODE) = NUM 109.
DEG(NXNODE) = - 1 110.
NXNODE = QLINK(NXNODE) 111.
IF (NXNODE .GT. 0) GOTO 600 112.
C 113.
IF ( RCHSZE .LE. 0 ) GO TO 800 114.
C ------------------------------------------------ 115.
C UPDATE THE DEGREES OF THE NODES IN THE REACHABLE 116.
C SET AND IDENTIFY INDISTINGUISHABLE NODES. 117.
C ------------------------------------------------ 118.
CALL QMDUPD ( XADJ, ADJNCY, RCHSZE, RCHSET, DEG, 119.
1 QSIZE, QLINK, MARKER, RCHSET(RCHSZE+1), 120.
1 NBRHD(NHDSZE+1) ) 121.
C ------------------------------------------- 122.
C RESET MARKER VALUE OF NODES IN REACH SET. 123.
C UPDATE THRESHOLD VALUE FOR CYCLIC SEARCH. 124.
C ALSO CALL QMDQT TO FORM NEW QUOTIENT GRAPH. 125.
C ------------------------------------------- 126.
MARKER(NODE) = 0 127.
DO 700 IRCH = 1, RCHSZE 128.
INODE = RCHSET(IRCH) 129.
IF ( MARKER(INODE) .LT. 0 ) GOTO 700 130.
MARKER(INODE) = 0 131.
NDEG = DEG(INODE) 132.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 133.
IF ( NDEG .GT. THRESH ) GOTO 700 134.
MINDEG = THRESH 135.
THRESH = NDEG 136.
SEARCH = INVP(INODE) 137.
700 CONTINUE 138.
IF ( NHDSZE .GT. 0 ) CALL QMDQT ( NODE, XADJ, 139.
1 ADJNCY, MARKER, RCHSZE, RCHSET, NBRHD ) 140.
800 IF ( NUM .LT. NEQNS ) GO TO 300 141.
RETURN 142.
END 143.

下面开始逐行解析脚本,首先初始化 DEG 向量和其它变量

MARKER 值全部设为 0,QSIZE 值和 QLINK 值全部为 1和 0,DEG 为原始图中的度数,MINDEG 是原始图中的最小度数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
C        -----------------------------------------------------            56.
C INITIALIZE DEGREE VECTOR AND OTHER WORKING VARIABLES. 57.
C ----------------------------------------------------- 58.
MINDEG = NEQNS 59.
NOFSUB = 0 60.
DO 100 NODE = 1, NEQNS 61.
PERM(NODE) = NODE 62.
INVP(NODE) = NODE 63.
MARKER(NODE) = 0 64.
QSIZE(NODE) = 1 65.
QLINK(NODE) = 0 66.
NDEG = XADJ(NODE+1) - XADJ(NODE) 67.
DEG(NODE) = NDEG 68.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 69.
100 CONTINUE 70.
NUM = 0 71.

这里采用 threshold search 来查找拥有最小度数的一个节点,SEARCH 变量指向查找的起始位置。

将之前 MINDEG 的值赋值给 THRESH ,MINDEG 重新设置为 NEQNS ,通过DO 400循环遍历需要搜索的所有节点,如果找到拥有最小度数的节点,则跳出循环,GO TO 500 。而 MINDEG 是跳出循环前大于 THRESH 的最小度数值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
C        -----------------------------------------------------            72.
C PERFORM THRESHOLD SEARCH TO GET A NODE OF MIN DEGREE. 73.
C VARIABLE SEARCH POINTS TO WHERE SEARCH SHOULD START. 74.
C ----------------------------------------------------- 75.
200 SEARCH = 1 76.
THRESH = MINDEG 77.
MINDEG = NEQNS 78.
300 NUMP1 = NUM + 1 79.
IF ( NUMP1 .GT. SEARCH ) SEARCH = NUMP1 80.
DO 400 J = SEARCH, NEQNS 81.
NODE = PERM(J) 82.
IF ( MARKER(NODE) .LT. 0 ) GOTO 400 83.
NDEG = DEG(NODE) 84.
IF ( NDEG .LE. THRESH ) GO TO 500 85.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 86.
400 CONTINUE 87.
GO TO 200 88.

对于最小度数节点,将其 MARKER 值设为1,通过调用 QMDRCH 子程序找到它通过已经消除的节点的 REACHABLE SETS,存储在 RCHSET 向量中(MARKER 值均修改为1),其向量长度为 RCHSZE ;与 NODE 节点相邻的已经消元的 supernodes 存储在 NBRHD 向量中(MARKER 值均修改为 -1),其向量长度为 NHDSZE 。

1
2
3
4
5
6
7
8
9
C           ---------------------------------------------------           89.
C NODE HAS MINIMUM DEGREE. FIND ITS REACHABLE SETS BY 90.
C CALLING QMDRCH. 91.
C --------------------------------------------------- 92.
500 SEARCH = J 93.
NOFSUB = NOFSUB + DEG(NODE) 94.
MARKER(NODE) = 1 95.
CALL QMDRCH (NODE, XADJ, ADJNCY, DEG, MARKER, 96.
1 RCHSZE, RCHSET, NHDSZE, NBRHD ) 97.

消除与 NODE 节点所有不可区分的节点,这些节点为 NODE, QLINK(NODE), … 等。

将 NODE 赋值给 NXNODE ,NUM 加 1,应该是 NXNODE 排序后的新位置。

NP 是 NXNODE 节点的原始位置,IP 是 NUM 位置现在的节点。

下面四句互换 NUM 和 NP 位置的节点,PERM(NP) = IPINVP(IP) = NP 将 IP 节点放置到 NP 位置。PERM(NUM) = NXNODEINVP(NXNODE) = NUM 将 NXNODE 放置到 NUM 位置。

排好序的 NXNODE 节点就可以消除了,消除方式是将其 DEG 值设为 -1 。

下面我们获取与 NXNODE 节点不可区分的节点,一并消除。NXNODE = QLINK(NXNODE) 此命令将 NXNODE 更新为下一个不可区分的节点。IF (NXNODE .GT. 0) GO TO 600 ,如果其不为0,则再进行循环中进行编号消除。

将所有不可区分的节点消除完毕之后,如果 researchable sets 长度为0,则 GO TO 800,跳过度数更新这一步,直接进入下一个节点的排序。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
C           ------------------------------------------------              98.
C ELIMINATE ALL NODES INDISTINGUISHABLE FROM NODE. 99.
C THEY ARE GIVEN BY NODE, QLINK(NODE), .... 100.
C ------------------------------------------------ 101.
NXNODE = NODE 102.
600 NUM = NUM + 1 103.
NP = INVP(NXNODE) 104.
IP = PERM(NUM) 105.
PERM(NP) = IP 106.
INVP(IP) = NP 107.
PERM(NUM) = NXNODE 108.
INVP(NXNODE) = NUM 109.
DEG(NXNODE) = - 1 110.
NXNODE = QLINK(NXNODE) 111.
IF (NXNODE .GT. 0) GOTO 600 112.
C 113.
IF ( RCHSZE .LE. 0 ) GO TO 800 114.

如果 researchable sets 长度不为0,我们则对其 researchable sets 中节点的度数进行更新,并查找不可区分的节点。这里我们通过调用子程序 QMDUPD 来实现这一点(此时 researchable sets 节点的 MARKER 值均为1)。最后两个参数 RCHSET(RCHSZE+1)NBRHD(NHDSZE+1) 利用 RCHSET 和 NBRHD 中尚未利用的空间作为临时变量。

运行结束后,确定新的不可区分的节点,修改 RCHSET 中的 DEG 值。RCHSET 中除了被合并到不可区分的 supernode 的节点的 MARKER 值改为-1之外,其它 MARKER 值改为2

1
2
3
4
5
6
7
C              ------------------------------------------------          115.
C UPDATE THE DEGREES OF THE NODES IN THE REACHABLE 116.
C SET AND IDENTIFY INDISTINGUISHABLE NODES. 117.
C ------------------------------------------------ 118.
CALL QMDUPD ( XADJ, ADJNCY, RCHSZE, RCHSET, DEG, 119.
1 QSIZE, QLINK, MARKER, RCHSET(RCHSZE+1), 120.
1 NBRHD(NHDSZE+1) ) 121.

将要消元的节点 NODE 的 MARKER 值修改为 0,并重置 RCHSET 中节点的 MARKER 值,遍历 RCHSET 中的节点 INODE,如果其 MARKER 值小于0,说明其已经合并到其它节点中了,则跳过这个节点,不然的话将其 MARKER 值重置为 0。节点 INODE 的 DEG 值赋值给 NDEG,如果 NDEG 小于 MINDEG,则更新 MINDEG。

这一步完成后,所有节点的 MAKER 值均已经正确设置(NODE 的 MARKER 值为0,RCHSET 除了被合并的节点之外 MARKER 值也为 0,NBRHD 的 MARKER 值仍为 -1 ,因此这些节点之后会合并到 NODE 中)。

如果 NDEG 小于等于 THRESH,则将MINDEG 改为旧的 THRESH 值,将 THRESH 改为 NDEG 值,SEARCH 改为 INODE 所处的位置。此时说明 INODE 就是下一个消元的节点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
C              -------------------------------------------               122.
C RESET MARKER VALUE OF NODES IN REACH SET. 123.
C UPDATE THRESHOLD VALUE FOR CYCLIC SEARCH. 124.
C ALSO CALL QMDQT TO FORM NEW QUOTIENT GRAPH. 125.
C ------------------------------------------- 126.
MARKER(NODE) = 0 127.
DO 700 IRCH = 1, RCHSZE 128.
INODE = RCHSET(IRCH) 129.
IF ( MARKER(INODE) .LT. 0 ) GOTO 700 130.
MARKER(INODE) = 0 131.
NDEG = DEG(INODE) 132.
IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 133.
IF ( NDEG .GT. THRESH ) GOTO 700 134.
MINDEG = THRESH 135.
THRESH = NDEG 136.
SEARCH = INVP(INODE) 137.
700 CONTINUE 138.

如果 NHDSZE 大于0,说明 NODE 需要和其它消元节点进行合并,调用 QMDQT 来形成新的商图。

1
2
          IF ( NHDSZE .GT. 0 )  CALL  QMDQT ( NODE, XADJ,           139.
1 ADJNCY, MARKER, RCHSZE, RCHSET, NBRHD ) 140.

最后,如果 NUM 小于 NEQNS,GO TO 300 ,重新进入循环。

1
2
3
800    IF ( NUM .LT. NEQNS )  GO TO 300                                141.
RETURN 142.
END 143.

讨论

这里用的所谓的 threshold search 查找度数最小的节点的算法写得并不很合理,这一步实际上就是从尚未消元的节点中,找到一个度数最小的节点,内容其实很简单,完全可以直接用一遍循环得到,原76至95行可以修改如下。

遍历过程中如果存在节点度数小于等于上次的最小度数 THRESH ,则这个节点就是本次消元的节点,GO TO 500(THRESH 不变) 。否则的话遍历结束后,MINDEG 和 NODE 就时最小度数和相应节点,最后再将 MINDEG 赋值给 THRESH 。这样这一块最多遍历一遍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
          THRESH = MINDEG                                              
300 MINDEG = NEQNS
NUMP1 = NUM + 1
DO 400 J = NUMP1, NEQNS
INODE = PERM(J)
IF ( MARKER(INODE) .LT. 0 ) GOTO 400
NDEG = DEG(INODE)
IF ( NDEG .LE. THRESH ) THEN
NODE = INODE
GO TO 500
END IF
IF ( NDEG .LT. MINDEG ) THEN
MINDEG = NDEG
NODE = INODE
ENDIF
400 CONTINUE
THRESH = MINDEG

500 NOFSUB = NOFSUB + DEG(NODE)
...

原脚本中采用 SEARCH,MINDEG 和 THRESH 三个变量来查找最小位点,那么查找消元节点的过程如下。第一次查找到的节点假设为 ,MINDEG 为原始 1到 j-1 节点中的原始最小度数,SEARCH 设为 J 。节点 1 和 j 互换后,MINDEG 就是PERM 中 2 到 j 位置节点中的原始最小度数。更新度数后,如果 RCHSET 中节点度数小于 MINDEG ,更新 MINDEG 。 假设我们最后重置 RCHSET 节点的 MARKER 值时均没有找到下一个消元的节点(没有更新 THRESH ,注意此时仍为第一次的最小度数), 之后我们查找第二个消元节点,GO TO 300 ,遍历 PERM 中第 J 到 最后一个位置的所有节点,如果没有度数小于等于 THRESH (由于 THRESH 仍是上一次的最小度数,很可能就是找不到的),遍历完 400 循环之后, MINDEG 值为最小度数的下限GO TO 200 , 更新 THRESH 值为 MINDEG,再次进入 400 循环 ,这一次大概率可以找到需要的节点 NODE ,退出循环,GO TO 500 。但是还是有一个小概率遍历完还是没找到,但是遍历完之后, MINDEG 就会更新为所有剩余节点的实际的最小度数,GO TO 200 ,再经过一次 400 循环 ,这一次绝对可以找到需要的节点 NODE,然后跳出循环。

因此这里查找最小度数过程中可能经过多次循环(最佳一次,最差三次),而且可读性差,逻辑性差。

参考文献

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

请我喝杯茶吧~

支付宝
微信