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