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