Mercurial > octave
annotate src/DLD-FUNCTIONS/symrcm.cc @ 11553:01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Use same variable names in error() strings and in documentation.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sun, 16 Jan 2011 22:13:23 -0800 |
parents | fd0a3ac60b0e |
children | 12df7854fa7c |
rev | line source |
---|---|
6608 | 1 /* |
2 | |
11523 | 3 Copyright (C) 2007-2011 Michael Weitzel |
6608 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
6608 | 11 |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
6608 | 20 |
21 */ | |
22 | |
23 /* | |
24 An implementation of the Reverse Cuthill-McKee algorithm (symrcm) | |
25 | |
26 The implementation of this algorithm is based in the descriptions found in | |
27 | |
28 @INPROCEEDINGS{, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
29 author = {E. Cuthill and J. McKee}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
30 title = {Reducing the Bandwidth of Sparse Symmetric Matrices}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
31 booktitle = {Proceedings of the 24th ACM National Conference}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
32 publisher = {Brandon Press}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
33 pages = {157 -- 172}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
34 location = {New Jersey}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
35 year = {1969} |
6608 | 36 } |
37 | |
38 @BOOK{, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
39 author = {Alan George and Joseph W. H. Liu}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
40 title = {Computer Solution of Large Sparse Positive Definite Systems}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
41 publisher = {Prentice Hall Series in Computational Mathematics}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
42 ISBN = {0-13-165274-5}, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
43 year = {1981} |
6608 | 44 } |
45 | |
46 The algorithm represents a heuristic approach to the NP-complete minimum | |
47 bandwidth problem. | |
48 | |
6610 | 49 Written by Michael Weitzel <michael.weitzel@@uni-siegen.de> |
50 <weitzel@@ldknet.org> | |
6608 | 51 */ |
52 | |
53 #ifdef HAVE_CONFIG_H | |
54 #include <config.h> | |
55 #endif | |
56 | |
57 #include "ov.h" | |
58 #include "defun-dld.h" | |
59 #include "error.h" | |
60 #include "gripes.h" | |
61 #include "utils.h" | |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
62 #include "oct-locbuf.h" |
6608 | 63 |
64 #include "ov-re-mat.h" | |
65 #include "ov-re-sparse.h" | |
66 #include "ov-cx-sparse.h" | |
67 #include "oct-sparse.h" | |
68 | |
69 // A node struct for the Cuthill-McKee algorithm | |
70 struct CMK_Node | |
71 { | |
72 // the node's id (matrix row index) | |
73 octave_idx_type id; | |
74 // the node's degree | |
75 octave_idx_type deg; | |
76 // minimal distance to the root of the spanning tree | |
77 octave_idx_type dist; | |
78 }; | |
79 | |
80 // A simple queue. | |
81 // Queues Q have a fixed maximum size N (rows,cols of the matrix) and are | |
82 // stored in an array. qh and qt point to queue head and tail. | |
83 | |
84 // Enqueue operation (adds a node "o" at the tail) | |
6959 | 85 |
6608 | 86 inline static void |
6959 | 87 Q_enq (CMK_Node *Q, octave_idx_type N, octave_idx_type& qt, const CMK_Node& o) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
88 { |
6959 | 89 Q[qt] = o; |
90 qt = (qt + 1) % (N + 1); | |
91 } | |
6608 | 92 |
93 // Dequeue operation (removes a node from the head) | |
6959 | 94 |
6608 | 95 inline static CMK_Node |
6959 | 96 Q_deq (CMK_Node * Q, octave_idx_type N, octave_idx_type& qh) |
97 { | |
98 CMK_Node r = Q[qh]; | |
99 qh = (qh + 1) % (N + 1); | |
100 return r; | |
101 } | |
6608 | 102 |
103 // Predicate (queue empty) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
104 #define Q_empty(Q, N, qh, qt) ((qh) == (qt)) |
6608 | 105 |
106 // A simple, array-based binary heap (used as a priority queue for nodes) | |
107 | |
108 // the left descendant of entry i | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
109 #define LEFT(i) (((i) << 1) + 1) // = (2*(i)+1) |
6608 | 110 // the right descendant of entry i |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
111 #define RIGHT(i) (((i) << 1) + 2) // = (2*(i)+2) |
6608 | 112 // the parent of entry i |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
113 #define PARENT(i) (((i) - 1) >> 1) // = floor(((i)-1)/2) |
6608 | 114 |
115 // Builds a min-heap (the root contains the smallest element). A is an array | |
116 // with the graph's nodes, i is a starting position, size is the length of A. | |
6959 | 117 |
6608 | 118 static void |
6959 | 119 H_heapify_min (CMK_Node *A, octave_idx_type i, octave_idx_type size) |
120 { | |
121 octave_idx_type j = i; | |
122 for (;;) | |
123 { | |
124 octave_idx_type l = LEFT(j); | |
125 octave_idx_type r = RIGHT(j); | |
126 | |
127 octave_idx_type smallest; | |
128 if (l < size && A[l].deg < A[j].deg) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
129 smallest = l; |
6959 | 130 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
131 smallest = j; |
6959 | 132 |
133 if (r < size && A[r].deg < A[smallest].deg) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
134 smallest = r; |
6959 | 135 |
136 if (smallest != j) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
137 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
138 CMK_Node tmp = A[j]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
139 A[j] = A[smallest]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
140 A[smallest] = tmp; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
141 j = smallest; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
142 } |
6959 | 143 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
144 break; |
6959 | 145 } |
146 } | |
6608 | 147 |
148 // Heap operation insert. Running time is O(log(n)) | |
6959 | 149 |
6608 | 150 static void |
6959 | 151 H_insert (CMK_Node *H, octave_idx_type& h, const CMK_Node& o) |
152 { | |
153 octave_idx_type i = h++; | |
154 | |
155 H[i] = o; | |
156 | |
157 if (i == 0) | |
158 return; | |
159 do | |
160 { | |
161 octave_idx_type p = PARENT(i); | |
162 if (H[i].deg < H[p].deg) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
163 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
164 CMK_Node tmp = H[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
165 H[i] = H[p]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
166 H[p] = tmp; |
6959 | 167 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
168 i = p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
169 } |
6959 | 170 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
171 break; |
6959 | 172 } |
173 while (i > 0); | |
174 } | |
6608 | 175 |
176 // Heap operation remove-min. Removes the smalles element in O(1) and | |
177 // reorganizes the heap optionally in O(log(n)) | |
6959 | 178 |
6608 | 179 inline static CMK_Node |
6959 | 180 H_remove_min (CMK_Node *H, octave_idx_type& h, int reorg/*=1*/) |
181 { | |
182 CMK_Node r = H[0]; | |
183 H[0] = H[--h]; | |
184 if (reorg) | |
185 H_heapify_min(H, 0, h); | |
186 return r; | |
187 } | |
6608 | 188 |
189 // Predicate (heap empty) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
190 #define H_empty(H, h) ((h) == 0) |
6608 | 191 |
192 // Helper function for the Cuthill-McKee algorithm. Tries to determine a | |
193 // pseudo-peripheral node of the graph as starting node. | |
6959 | 194 |
6608 | 195 static octave_idx_type |
6959 | 196 find_starting_node (octave_idx_type N, const octave_idx_type *ridx, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
197 const octave_idx_type *cidx, const octave_idx_type *ridx2, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
198 const octave_idx_type *cidx2, octave_idx_type *D, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
199 octave_idx_type start) |
6959 | 200 { |
201 CMK_Node w; | |
202 | |
203 OCTAVE_LOCAL_BUFFER (CMK_Node, Q, N+1); | |
204 boolNDArray btmp (dim_vector (1, N), false); | |
205 bool *visit = btmp.fortran_vec (); | |
206 | |
207 octave_idx_type qh = 0; | |
208 octave_idx_type qt = 0; | |
209 CMK_Node x; | |
210 x.id = start; | |
211 x.deg = D[start]; | |
212 x.dist = 0; | |
213 Q_enq (Q, N, qt, x); | |
214 visit[start] = true; | |
215 | |
216 // distance level | |
217 octave_idx_type level = 0; | |
218 // current largest "eccentricity" | |
219 octave_idx_type max_dist = 0; | |
220 | |
221 for (;;) | |
222 { | |
223 while (! Q_empty (Q, N, qh, qt)) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
224 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
225 CMK_Node v = Q_deq (Q, N, qh); |
6959 | 226 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
227 if (v.dist > x.dist || (v.id != x.id && v.deg > x.deg)) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
228 x = v; |
6959 | 229 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
230 octave_idx_type i = v.id; |
6959 | 231 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
232 // add all unvisited neighbors to the queue |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
233 octave_idx_type j1 = cidx[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
234 octave_idx_type j2 = cidx2[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
235 while (j1 < cidx[i+1] || j2 < cidx2[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
236 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
237 OCTAVE_QUIT; |
6959 | 238 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
239 if (j1 == cidx[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
240 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
241 octave_idx_type r2 = ridx2[j2++]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
242 if (! visit[r2]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
243 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
244 // the distance of node j is dist(i)+1 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
245 w.id = r2; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
246 w.deg = D[r2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
247 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
248 Q_enq (Q, N, qt, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
249 visit[r2] = true; |
6959 | 250 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
251 if (w.dist > level) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
252 level = w.dist; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
253 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
254 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
255 else if (j2 == cidx2[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
256 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
257 octave_idx_type r1 = ridx[j1++]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
258 if (! visit[r1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
259 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
260 // the distance of node j is dist(i)+1 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
261 w.id = r1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
262 w.deg = D[r1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
263 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
264 Q_enq (Q, N, qt, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
265 visit[r1] = true; |
6959 | 266 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
267 if (w.dist > level) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
268 level = w.dist; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
269 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
270 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
271 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
272 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
273 octave_idx_type r1 = ridx[j1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
274 octave_idx_type r2 = ridx2[j2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
275 if (r1 <= r2) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
276 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
277 if (! visit[r1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
278 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
279 w.id = r1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
280 w.deg = D[r1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
281 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
282 Q_enq (Q, N, qt, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
283 visit[r1] = true; |
6959 | 284 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
285 if (w.dist > level) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
286 level = w.dist; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
287 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
288 j1++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
289 if (r1 == r2) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
290 j2++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
291 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
292 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
293 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
294 if (! visit[r2]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
295 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
296 w.id = r2; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
297 w.deg = D[r2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
298 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
299 Q_enq (Q, N, qt, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
300 visit[r2] = true; |
6959 | 301 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
302 if (w.dist > level) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
303 level = w.dist; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
304 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
305 j2++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
306 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
307 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
308 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
309 } // finish of BFS |
6959 | 310 |
311 if (max_dist < x.dist) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
312 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
313 max_dist = x.dist; |
6959 | 314 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
315 for (octave_idx_type i = 0; i < N; i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
316 visit[i] = false; |
6959 | 317 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
318 visit[x.id] = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
319 x.dist = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
320 qt = qh = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
321 Q_enq (Q, N, qt, x); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
322 } |
6959 | 323 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
324 break; |
6959 | 325 } |
326 return x.id; | |
327 } | |
6608 | 328 |
329 // Calculates the node's degrees. This means counting the non-zero elements | |
330 // in the symmetric matrix' rows. This works for non-symmetric matrices | |
331 // as well. | |
6959 | 332 |
333 static octave_idx_type | |
334 calc_degrees (octave_idx_type N, const octave_idx_type *ridx, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
335 const octave_idx_type *cidx, octave_idx_type *D) |
6959 | 336 { |
337 octave_idx_type max_deg = 0; | |
338 | |
339 for (octave_idx_type i = 0; i < N; i++) | |
340 D[i] = 0; | |
341 | |
342 for (octave_idx_type j = 0; j < N; j++) | |
343 { | |
344 for (octave_idx_type i = cidx[j]; i < cidx[j+1]; i++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
345 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
346 OCTAVE_QUIT; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
347 octave_idx_type k = ridx[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
348 // there is a non-zero element (k,j) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
349 D[k]++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
350 if (D[k] > max_deg) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
351 max_deg = D[k]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
352 // if there is no element (j,k) there is one in |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
353 // the symmetric matrix: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
354 if (k != j) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
355 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
356 bool found = false; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
357 for (octave_idx_type l = cidx[k]; l < cidx[k + 1]; l++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
358 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
359 OCTAVE_QUIT; |
6959 | 360 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
361 if (ridx[l] == j) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
362 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
363 found = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
364 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
365 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
366 else if (ridx[l] > j) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
367 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
368 } |
6959 | 369 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
370 if (! found) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
371 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
372 // A(j,k) == 0 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
373 D[j]++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
374 if (D[j] > max_deg) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
375 max_deg = D[j]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
376 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
377 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
378 } |
6959 | 379 } |
380 return max_deg; | |
381 } | |
6608 | 382 |
383 // Transpose of the structure of a square sparse matrix | |
6959 | 384 |
6608 | 385 static void |
386 transpose (octave_idx_type N, const octave_idx_type *ridx, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
387 const octave_idx_type *cidx, octave_idx_type *ridx2, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
388 octave_idx_type *cidx2) |
6959 | 389 { |
390 octave_idx_type nz = cidx[N]; | |
391 | |
392 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, N + 1); | |
393 for (octave_idx_type i = 0; i < N; i++) | |
394 w[i] = 0; | |
395 for (octave_idx_type i = 0; i < nz; i++) | |
396 w[ridx[i]]++; | |
397 nz = 0; | |
398 for (octave_idx_type i = 0; i < N; i++) | |
399 { | |
400 OCTAVE_QUIT; | |
401 cidx2[i] = nz; | |
402 nz += w[i]; | |
403 w[i] = cidx2[i]; | |
404 } | |
405 cidx2[N] = nz; | |
406 w[N] = nz; | |
407 | |
408 for (octave_idx_type j = 0; j < N; j++) | |
409 for (octave_idx_type k = cidx[j]; k < cidx[j + 1]; k++) | |
410 { | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
411 OCTAVE_QUIT; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
412 octave_idx_type q = w [ridx[k]]++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
413 ridx2[q] = j; |
6959 | 414 } |
415 } | |
6608 | 416 |
417 // An implementation of the Cuthill-McKee algorithm. | |
418 DEFUN_DLD (symrcm, args, , | |
419 "-*- texinfo -*-\n\ | |
7650 | 420 @deftypefn {Loadable Function} {@var{p} =} symrcm (@var{S})\n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
421 Return the symmetric reverse Cuthill-McKee permutation of @var{S}.\n\ |
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
422 @var{p} is a permutation vector such that\n\ |
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
423 @code{@var{S}(@var{p}, @var{p})} tends to have its diagonal elements\n\ |
6610 | 424 closer to the diagonal than @var{S}. This is a good preordering for LU\n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
425 or Cholesky@tie{}factorization of matrices that come from 'long, skinny'\n\ |
6610 | 426 problems. It works for both symmetric and asymmetric @var{S}.\n\ |
6608 | 427 \n\ |
428 The algorithm represents a heuristic approach to the NP-complete\n\ | |
6610 | 429 bandwidth minimization problem. The implementation is based in the\n\ |
6608 | 430 descriptions found in\n\ |
431 \n\ | |
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10154
diff
changeset
|
432 E. Cuthill, J. McKee. @cite{Reducing the Bandwidth of Sparse Symmetric\n\ |
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10154
diff
changeset
|
433 Matrices}. Proceedings of the 24th ACM National Conference, 157--172\n\ |
6610 | 434 1969, Brandon Press, New Jersey.\n\ |
6608 | 435 \n\ |
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10154
diff
changeset
|
436 A. George, J.W.H. Liu. @cite{Computer Solution of Large Sparse\n\ |
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10154
diff
changeset
|
437 Positive Definite Systems}, Prentice Hall Series in Computational\n\ |
6610 | 438 Mathematics, ISBN 0-13-165274-5, 1981.\n\ |
6608 | 439 \n\ |
440 @seealso{colperm, colamd, symamd}\n\ | |
441 @end deftypefn") | |
442 { | |
443 octave_value retval; | |
6959 | 444 int nargin = args.length (); |
6608 | 445 |
446 if (nargin != 1) | |
447 { | |
448 print_usage (); | |
449 return retval; | |
450 } | |
451 | |
6959 | 452 octave_value arg = args(0); |
6608 | 453 |
454 // the parameter of the matrix is converted into a sparse matrix | |
455 //(if necessary) | |
456 octave_idx_type *cidx; | |
457 octave_idx_type *ridx; | |
458 SparseMatrix Ar; | |
459 SparseComplexMatrix Ac; | |
460 | |
461 if (arg.is_real_type ()) | |
462 { | |
6959 | 463 Ar = arg.sparse_matrix_value (); |
6608 | 464 // Note cidx/ridx are const, so use xridx and xcidx... |
465 cidx = Ar.xcidx (); | |
466 ridx = Ar.xridx (); | |
467 } | |
468 else | |
469 { | |
6959 | 470 Ac = arg.sparse_complex_matrix_value (); |
6608 | 471 cidx = Ac.xcidx (); |
472 ridx = Ac.xridx (); | |
473 } | |
474 | |
475 if (error_state) | |
476 return retval; | |
477 | |
478 octave_idx_type nr = arg.rows (); | |
479 octave_idx_type nc = arg.columns (); | |
480 | |
481 if (nr != nc) | |
482 { | |
6959 | 483 gripe_square_matrix_required ("symrcm"); |
6608 | 484 return retval; |
485 } | |
486 | |
487 if (nr == 0 && nc == 0) | |
6959 | 488 return octave_value (NDArray (dim_vector (1, 0))); |
489 | |
6608 | 490 // sizes of the heaps |
491 octave_idx_type s = 0; | |
6959 | 492 |
6608 | 493 // head- and tail-indices for the queue |
6959 | 494 octave_idx_type qt = 0, qh = 0; |
6608 | 495 CMK_Node v, w; |
496 // dimension of the matrix | |
6959 | 497 octave_idx_type N = nr; |
498 | |
6608 | 499 OCTAVE_LOCAL_BUFFER (octave_idx_type, cidx2, N + 1); |
500 OCTAVE_LOCAL_BUFFER (octave_idx_type, ridx2, cidx[N]); | |
501 transpose (N, ridx, cidx, ridx2, cidx2); | |
502 | |
503 // the permutation vector | |
504 NDArray P (dim_vector (1, N)); | |
505 | |
506 // compute the node degrees | |
6959 | 507 OCTAVE_LOCAL_BUFFER (octave_idx_type, D, N); |
508 octave_idx_type max_deg = calc_degrees (N, ridx, cidx, D); | |
6608 | 509 |
510 // if none of the nodes has a degree > 0 (a matrix of zeros) | |
511 // the return value corresponds to the identity permutation | |
512 if (max_deg == 0) | |
513 { | |
6959 | 514 for (octave_idx_type i = 0; i < N; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
515 P(i) = i; |
6959 | 516 return octave_value (P); |
6608 | 517 } |
518 | |
519 // a heap for the a node's neighbors. The number of neighbors is | |
520 // limited by the maximum degree max_deg: | |
6959 | 521 OCTAVE_LOCAL_BUFFER (CMK_Node, S, max_deg); |
6608 | 522 |
523 // a queue for the BFS. The array is always one element larger than | |
524 // the number of entries that are stored. | |
6959 | 525 OCTAVE_LOCAL_BUFFER (CMK_Node, Q, N+1); |
6608 | 526 |
527 // a counter (for building the permutation) | |
6959 | 528 octave_idx_type c = -1; |
6608 | 529 |
6959 | 530 // upper bound for the bandwidth (=quality of solution) |
6608 | 531 // initialize the bandwidth of the graph with 0. B contains the |
532 // the maximum of the theoretical lower limits of the subgraphs | |
533 // bandwidths. | |
6959 | 534 octave_idx_type B = 0; |
6608 | 535 |
536 // mark all nodes as unvisited; with the exception of the nodes | |
537 // that have degree==0 and build a CC of the graph. | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
538 |
6608 | 539 boolNDArray btmp (dim_vector (1, N), false); |
540 bool *visit = btmp.fortran_vec (); | |
541 | |
542 do | |
543 { | |
544 // locate an unvisited starting node of the graph | |
6959 | 545 octave_idx_type i; |
6608 | 546 for (i = 0; i < N; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
547 if (! visit[i]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
548 break; |
6608 | 549 |
550 // locate a probably better starting node | |
551 v.id = find_starting_node (N, ridx, cidx, ridx2, cidx2, D, i); | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
552 |
6608 | 553 // mark the node as visited and enqueue it (a starting node |
554 // for the BFS). Since the node will be a root of a spanning | |
555 // tree, its dist is 0. | |
556 v.deg = D[v.id]; | |
557 v.dist = 0; | |
558 visit[v.id] = true; | |
6959 | 559 Q_enq (Q, N, qt, v); |
6608 | 560 |
6959 | 561 // lower bound for the bandwidth of a subgraph |
6608 | 562 // keep a "level" in the spanning tree (= min. distance to the |
563 // root) for determining the bandwidth of the computed | |
564 // permutation P | |
6959 | 565 octave_idx_type Bsub = 0; |
6608 | 566 // min. dist. to the root is 0 |
6959 | 567 octave_idx_type level = 0; |
6608 | 568 // the root is the first/only node on level 0 |
6959 | 569 octave_idx_type level_N = 1; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
570 |
6677 | 571 while (! Q_empty (Q, N, qh, qt)) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
572 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
573 v = Q_deq (Q, N, qh); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
574 i = v.id; |
6608 | 575 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
576 c++; |
6608 | 577 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
578 // for computing the inverse permutation P where |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
579 // A(inv(P),inv(P)) or P'*A*P is banded |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
580 // P(i) = c; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
581 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
582 // for computing permutation P where |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
583 // A(P(i),P(j)) or P*A*P' is banded |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
584 P(c) = i; |
6608 | 585 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
586 // put all unvisited neighbors j of node i on the heap |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
587 s = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
588 octave_idx_type j1 = cidx[i]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
589 octave_idx_type j2 = cidx2[i]; |
6608 | 590 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
591 OCTAVE_QUIT; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
592 while (j1 < cidx[i+1] || j2 < cidx2[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
593 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
594 OCTAVE_QUIT; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
595 if (j1 == cidx[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
596 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
597 octave_idx_type r2 = ridx2[j2++]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
598 if (! visit[r2]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
599 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
600 // the distance of node j is dist(i)+1 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
601 w.id = r2; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
602 w.deg = D[r2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
603 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
604 H_insert(S, s, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
605 visit[r2] = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
606 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
607 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
608 else if (j2 == cidx2[i+1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
609 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
610 octave_idx_type r1 = ridx[j1++]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
611 if (! visit[r1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
612 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
613 w.id = r1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
614 w.deg = D[r1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
615 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
616 H_insert(S, s, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
617 visit[r1] = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
618 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
619 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
620 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
621 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
622 octave_idx_type r1 = ridx[j1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
623 octave_idx_type r2 = ridx2[j2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
624 if (r1 <= r2) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
625 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
626 if (! visit[r1]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
627 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
628 w.id = r1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
629 w.deg = D[r1]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
630 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
631 H_insert(S, s, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
632 visit[r1] = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
633 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
634 j1++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
635 if (r1 == r2) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
636 j2++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
637 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
638 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
639 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
640 if (! visit[r2]) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
641 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
642 w.id = r2; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
643 w.deg = D[r2]; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
644 w.dist = v.dist+1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
645 H_insert(S, s, w); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
646 visit[r2] = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
647 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
648 j2++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
649 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
650 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
651 } |
6608 | 652 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
653 // add the neighbors to the queue (sorted by node degree) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
654 while (! H_empty (S, s)) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
655 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
656 OCTAVE_QUIT; |
6608 | 657 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
658 // locate a neighbor of i with minimal degree in O(log(N)) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
659 v = H_remove_min(S, s, 1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
660 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
661 // entered the BFS a new level? |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
662 if (v.dist > level) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
663 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
664 // adjustment of bandwith: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
665 // "[...] the minimum bandwidth that |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
666 // can be obtained [...] is the |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
667 // maximum number of nodes per level" |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
668 if (Bsub < level_N) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
669 Bsub = level_N; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
670 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
671 level = v.dist; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
672 // v is the first node on the new level |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
673 level_N = 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
674 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
675 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
676 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
677 // there is no new level but another node on |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
678 // this level: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
679 level_N++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
680 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
681 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
682 // enqueue v in O(1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
683 Q_enq (Q, N, qt, v); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
684 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
685 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
686 // synchronize the bandwidth with level_N once again: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
687 if (Bsub < level_N) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
688 Bsub = level_N; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
689 } |
6608 | 690 // finish of BFS. If there are still unvisited nodes in the graph |
691 // then it is split into CCs. The computed bandwidth is the maximum | |
692 // of all subgraphs. Update: | |
693 if (Bsub > B) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9066
diff
changeset
|
694 B = Bsub; |
6608 | 695 } |
696 // are there any nodes left? | |
697 while (c+1 < N); | |
698 | |
699 // compute the reverse-ordering | |
700 s = N / 2 - 1; | |
6959 | 701 for (octave_idx_type i = 0, j = N - 1; i <= s; i++, j--) |
6608 | 702 { |
703 double tmp = P.elem(i); | |
704 P.elem(i) = P.elem(j); | |
705 P.elem(j) = tmp; | |
706 } | |
707 | |
708 // increment all indices, since Octave is not C | |
6959 | 709 return octave_value (P+1); |
6608 | 710 } |