Mercurial > octave
annotate src/DLD-FUNCTIONS/symrcm.cc @ 8377:25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 29 Oct 2008 16:52:10 +0100 |
| parents | eb7bdde776f2 |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 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" | |
|
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) |
| 88 { | |
| 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) | |
| 104 #define Q_empty(Q, N, qh, qt) ((qh) == (qt)) | |
| 105 | |
| 106 // A simple, array-based binary heap (used as a priority queue for nodes) | |
| 107 | |
| 108 // the left descendant of entry i | |
| 109 #define LEFT(i) (((i) << 1) + 1) // = (2*(i)+1) | |
| 110 // the right descendant of entry i | |
| 111 #define RIGHT(i) (((i) << 1) + 2) // = (2*(i)+2) | |
| 112 // the parent of entry i | |
| 113 #define PARENT(i) (((i) - 1) >> 1) // = floor(((i)-1)/2) | |
| 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) | |
| 129 smallest = l; | |
| 130 else | |
| 131 smallest = j; | |
| 132 | |
| 133 if (r < size && A[r].deg < A[smallest].deg) | |
| 134 smallest = r; | |
| 135 | |
| 136 if (smallest != j) | |
| 137 { | |
| 138 CMK_Node tmp = A[j]; | |
| 139 A[j] = A[smallest]; | |
| 140 A[smallest] = tmp; | |
| 141 j = smallest; | |
| 142 } | |
| 143 else | |
| 144 break; | |
| 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) | |
| 163 { | |
| 164 CMK_Node tmp = H[i]; | |
| 165 H[i] = H[p]; | |
| 166 H[p] = tmp; | |
| 167 | |
| 168 i = p; | |
| 169 } | |
| 170 else | |
| 171 break; | |
| 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) | |
| 190 #define H_empty(H, h) ((h) == 0) | |
| 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, |
| 197 const octave_idx_type *cidx, const octave_idx_type *ridx2, | |
| 198 const octave_idx_type *cidx2, octave_idx_type *D, | |
| 199 octave_idx_type start) | |
| 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)) | |
| 224 { | |
| 225 CMK_Node v = Q_deq (Q, N, qh); | |
| 226 | |
| 227 if (v.dist > x.dist || (v.id != x.id && v.deg > x.deg)) | |
| 228 x = v; | |
| 229 | |
| 230 octave_idx_type i = v.id; | |
| 231 | |
| 232 // add all unvisited neighbors to the queue | |
| 233 octave_idx_type j1 = cidx[i]; | |
| 234 octave_idx_type j2 = cidx2[i]; | |
| 235 while (j1 < cidx[i+1] || j2 < cidx2[i+1]) | |
| 236 { | |
| 237 OCTAVE_QUIT; | |
| 238 | |
| 239 if (j1 == cidx[i+1]) | |
| 240 { | |
| 241 octave_idx_type r2 = ridx2[j2++]; | |
| 242 if (! visit[r2]) | |
| 243 { | |
| 244 // the distance of node j is dist(i)+1 | |
| 245 w.id = r2; | |
| 246 w.deg = D[r2]; | |
| 247 w.dist = v.dist+1; | |
| 248 Q_enq (Q, N, qt, w); | |
| 249 visit[r2] = true; | |
| 250 | |
| 251 if (w.dist > level) | |
| 252 level = w.dist; | |
| 253 } | |
| 254 } | |
| 255 else if (j2 == cidx2[i+1]) | |
| 256 { | |
| 257 octave_idx_type r1 = ridx[j1++]; | |
| 258 if (! visit[r1]) | |
| 259 { | |
| 260 // the distance of node j is dist(i)+1 | |
| 261 w.id = r1; | |
| 262 w.deg = D[r1]; | |
| 263 w.dist = v.dist+1; | |
| 264 Q_enq (Q, N, qt, w); | |
| 265 visit[r1] = true; | |
| 266 | |
| 267 if (w.dist > level) | |
| 268 level = w.dist; | |
| 269 } | |
| 270 } | |
| 271 else | |
| 272 { | |
| 273 octave_idx_type r1 = ridx[j1]; | |
| 274 octave_idx_type r2 = ridx2[j2]; | |
| 275 if (r1 <= r2) | |
| 276 { | |
| 277 if (! visit[r1]) | |
| 278 { | |
| 279 w.id = r1; | |
| 280 w.deg = D[r1]; | |
| 281 w.dist = v.dist+1; | |
| 282 Q_enq (Q, N, qt, w); | |
| 283 visit[r1] = true; | |
| 284 | |
| 285 if (w.dist > level) | |
| 286 level = w.dist; | |
| 287 } | |
| 288 j1++; | |
| 289 if (r1 == r2) | |
| 290 j2++; | |
| 291 } | |
| 292 else | |
| 293 { | |
| 294 if (! visit[r2]) | |
| 295 { | |
| 296 w.id = r2; | |
| 297 w.deg = D[r2]; | |
| 298 w.dist = v.dist+1; | |
| 299 Q_enq (Q, N, qt, w); | |
| 300 visit[r2] = true; | |
| 301 | |
| 302 if (w.dist > level) | |
| 303 level = w.dist; | |
| 304 } | |
| 305 j2++; | |
| 306 } | |
| 307 } | |
| 308 } | |
| 309 } // finish of BFS | |
| 310 | |
| 311 if (max_dist < x.dist) | |
| 312 { | |
| 313 max_dist = x.dist; | |
| 314 | |
| 315 for (octave_idx_type i = 0; i < N; i++) | |
| 316 visit[i] = false; | |
| 317 | |
| 318 visit[x.id] = true; | |
| 319 x.dist = 0; | |
| 320 qt = qh = 0; | |
| 321 Q_enq (Q, N, qt, x); | |
| 322 } | |
| 323 else | |
| 324 break; | |
| 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, | |
| 335 const octave_idx_type *cidx, octave_idx_type *D) | |
| 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++) | |
| 345 { | |
| 346 OCTAVE_QUIT; | |
| 347 octave_idx_type k = ridx[i]; | |
| 348 // there is a non-zero element (k,j) | |
| 349 D[k]++; | |
| 350 if (D[k] > max_deg) | |
| 351 max_deg = D[k]; | |
| 352 // if there is no element (j,k) there is one in | |
| 353 // the symmetric matrix: | |
| 354 if (k != j) | |
| 355 { | |
| 356 bool found = false; | |
| 357 for (octave_idx_type l = cidx[k]; l < cidx[k + 1]; l++) | |
| 358 { | |
| 359 OCTAVE_QUIT; | |
| 360 | |
| 361 if (ridx[l] == j) | |
| 362 { | |
| 363 found = true; | |
| 364 break; | |
| 365 } | |
| 366 else if (ridx[l] > j) | |
| 367 break; | |
| 368 } | |
| 369 | |
| 370 if (! found) | |
| 371 { | |
| 372 // A(j,k) == 0 | |
| 373 D[j]++; | |
| 374 if (D[j] > max_deg) | |
| 375 max_deg = D[j]; | |
| 376 } | |
| 377 } | |
| 378 } | |
| 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, | |
| 387 const octave_idx_type *cidx, octave_idx_type *ridx2, | |
| 6959 | 388 octave_idx_type *cidx2) |
| 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 { | |
| 411 OCTAVE_QUIT; | |
| 412 octave_idx_type q = w [ridx[k]]++; | |
| 413 ridx2[q] = j; | |
| 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\ |
| 6608 | 421 Symmetric reverse Cuthill-McKee permutation of @var{S}.\n\ |
| 6610 | 422 Return a permutation vector @var{p} such that\n\ |
| 6608 | 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\ |
| 6608 | 425 or Cholesky 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\ | |
| 6610 | 432 E. Cuthill, J. McKee: Reducing the Bandwidth of Sparse Symmetric\n\ |
| 433 Matrices. Proceedings of the 24th ACM National Conference, 157-172\n\ | |
| 434 1969, Brandon Press, New Jersey.\n\ | |
| 6608 | 435 \n\ |
| 6610 | 436 Alan George, Joseph W. H. Liu: Computer Solution of Large Sparse\n\ |
| 437 Positive Definite Systems, Prentice Hall Series in Computational\n\ | |
| 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++) |
| 515 P(i) = i; | |
| 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. | |
| 538 | |
| 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++) |
| 6677 | 547 if (! visit[i]) |
| 6608 | 548 break; |
| 549 | |
| 550 // locate a probably better starting node | |
| 551 v.id = find_starting_node (N, ridx, cidx, ridx2, cidx2, D, i); | |
| 552 | |
| 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; |
| 6608 | 570 |
| 6677 | 571 while (! Q_empty (Q, N, qh, qt)) |
| 6608 | 572 { |
| 6959 | 573 v = Q_deq (Q, N, qh); |
| 6608 | 574 i = v.id; |
| 575 | |
| 576 c++; | |
| 577 | |
| 578 // for computing the inverse permutation P where | |
| 579 // A(inv(P),inv(P)) or P'*A*P is banded | |
| 580 // P(i) = c; | |
| 581 | |
| 582 // for computing permutation P where | |
| 583 // A(P(i),P(j)) or P*A*P' is banded | |
| 584 P(c) = i; | |
| 585 | |
| 586 // put all unvisited neighbors j of node i on the heap | |
| 587 s = 0; | |
| 588 octave_idx_type j1 = cidx[i]; | |
| 589 octave_idx_type j2 = cidx2[i]; | |
| 590 | |
| 591 OCTAVE_QUIT; | |
| 592 while (j1 < cidx[i+1] || j2 < cidx2[i+1]) | |
| 593 { | |
| 594 OCTAVE_QUIT; | |
| 595 if (j1 == cidx[i+1]) | |
| 596 { | |
| 597 octave_idx_type r2 = ridx2[j2++]; | |
| 6677 | 598 if (! visit[r2]) |
| 6608 | 599 { |
| 600 // the distance of node j is dist(i)+1 | |
| 601 w.id = r2; | |
| 602 w.deg = D[r2]; | |
| 603 w.dist = v.dist+1; | |
| 604 H_insert(S, s, w); | |
| 605 visit[r2] = true; | |
| 606 } | |
| 607 } | |
| 608 else if (j2 == cidx2[i+1]) | |
| 609 { | |
| 610 octave_idx_type r1 = ridx[j1++]; | |
| 6677 | 611 if (! visit[r1]) |
| 6608 | 612 { |
| 613 w.id = r1; | |
| 614 w.deg = D[r1]; | |
| 615 w.dist = v.dist+1; | |
| 616 H_insert(S, s, w); | |
| 617 visit[r1] = true; | |
| 618 } | |
| 619 } | |
| 620 else | |
| 621 { | |
| 622 octave_idx_type r1 = ridx[j1]; | |
| 623 octave_idx_type r2 = ridx2[j2]; | |
| 624 if (r1 <= r2) | |
| 625 { | |
| 6677 | 626 if (! visit[r1]) |
| 6608 | 627 { |
| 628 w.id = r1; | |
| 629 w.deg = D[r1]; | |
| 630 w.dist = v.dist+1; | |
| 631 H_insert(S, s, w); | |
| 632 visit[r1] = true; | |
| 633 } | |
| 634 j1++; | |
| 635 if (r1 == r2) | |
| 636 j2++; | |
| 637 } | |
| 638 else | |
| 639 { | |
| 6677 | 640 if (! visit[r2]) |
| 6608 | 641 { |
| 642 w.id = r2; | |
| 643 w.deg = D[r2]; | |
| 644 w.dist = v.dist+1; | |
| 645 H_insert(S, s, w); | |
| 646 visit[r2] = true; | |
| 647 } | |
| 648 j2++; | |
| 649 } | |
| 650 } | |
| 651 } | |
| 652 | |
| 653 // add the neighbors to the queue (sorted by node degree) | |
| 6959 | 654 while (! H_empty (S, s)) |
| 6608 | 655 { |
| 656 OCTAVE_QUIT; | |
| 657 | |
| 658 // locate a neighbor of i with minimal degree in O(log(N)) | |
| 659 v = H_remove_min(S, s, 1); | |
| 660 | |
| 661 // entered the BFS a new level? | |
| 662 if (v.dist > level) | |
| 663 { | |
| 664 // adjustment of bandwith: | |
| 665 // "[...] the minimum bandwidth that | |
| 666 // can be obtained [...] is the | |
| 667 // maximum number of nodes per level" | |
| 668 if (Bsub < level_N) | |
| 669 Bsub = level_N; | |
| 670 | |
| 671 level = v.dist; | |
| 672 // v is the first node on the new level | |
| 673 level_N = 1; | |
| 674 } | |
| 675 else | |
| 676 { | |
| 677 // there is no new level but another node on | |
| 678 // this level: | |
| 679 level_N++; | |
| 680 } | |
| 681 | |
| 682 // enqueue v in O(1) | |
| 6959 | 683 Q_enq (Q, N, qt, v); |
| 6608 | 684 } |
| 685 | |
| 686 // synchronize the bandwidth with level_N once again: | |
| 687 if (Bsub < level_N) | |
| 688 Bsub = level_N; | |
| 689 } | |
| 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) | |
| 694 B = Bsub; | |
| 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 } |