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