PaStiX Handbook  6.3.2
order_find_supernodes.c
Go to the documentation of this file.
1 /**
2  *
3  * @file order_find_supernodes.c
4  *
5  * PaStiX order function to find supernodes out of a given permutation.
6  *
7  * @copyright 2004-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
9  *
10  * @version 6.3.2
11  * @author Pascal Henon
12  * @author Mathieu Faverge
13  * @author Pierre Ramet
14  * @author Tony Delarue
15  * @date 2023-07-21
16  *
17  **/
18 #include "common.h"
19 #include "graph/graph.h"
20 #include "order/order_internal.h"
21
22 /**
23  *******************************************************************************
24  *
25  * @ingroup order_dev
26  *
27  * @brief Computes the size of each subtree.
28  *
29  *******************************************************************************
30  *
31  * @param[in] n
32  * The number of nodes.
33  *
34  * @param[in] father
35  * Array of size n.
36  * List of father to each node. If node i is a root then father[i] = i.
37  *
38  * @param[in] iperm
39  * Array of size n. The inverse permutation vector.
40  *
41  * @param[out] T
42  * Array of size n.
43  * On exit, size of each subtree.
44  *
45  *******************************************************************************/
46 static inline void
48  const pastix_int_t *father,
49  const pastix_int_t *iperm,
50  pastix_int_t *T)
51 {
52  /********************************************/
53  /* Compute the size of each subtree given */
54  /* the number of the father of each node */
55  /********************************************/
56  pastix_int_t k, i;
57
58  /* TODO pas la peine d'utiliser un tas; il suffit de parcourir iperm pour assurer
59  de toujours traiter un fils avant son pere */
60
61  memset( T, 0, sizeof(pastix_int_t) * n );
62
63  for(k=0;k<n;k++)
64  {
65  i = iperm[k];
66  T[i]++;
67  if(i!=father[i])
68  T[father[i]] += T[i];
69  }
70
71 #if defined(PASTIX_DEBUG_ORDERING)
72  {
73  pastix_int_t sum = 0;
74
75  for(i=0;i<n;i++) {
76  if(father[i] == i)
77  sum += T[i];
78  }
79  if(sum != n)
80  pastix_print_error( "compute_subtree_size: sum of the subtree = %ld n = %ld", (long)sum, (long)n );
81  assert(n == sum);
82  }
83 #endif
84 }
85
86 /**
87  *******************************************************************************
88  *
89  * @ingroup order_dev
90  *
91  * @brief Computes the post order of the elimination tree given on
92  * entry.
93  *
94  *******************************************************************************
95  *
96  * @param[in] n
97  * The number of nodes.
98  *
99  * @param[in] father
100  * Array of size n.
101  * List of father to each node. If node i is a root then father[i] = i.
102  *
103  * @param[inout] perm
104  * Array of size n. The permutation vector.
105  * On exit, the postorder permutation vector.
106  *
107  * @param[inout] invp
108  * Array of size n. The inverse permutation vector.
109  * On exit, the postorder inverse permutation vector.
110  *
111  * @param[inout] T
112  * Workspace of size n.
113  *
114  *******************************************************************************/
115 static inline void
117  const pastix_int_t *father,
118  pastix_int_t *perm,
119  pastix_int_t *invp,
120  pastix_int_t *T)
121 {
122  pastix_int_t i;
123  pastix_int_t j, k, t;
124
125  /*
126  * First compute the number of node in the subtree rooted in node i
127  */
128  compute_subtree_size(n, father, invp, T);
129
130  /*
131  * When multiple roots are present we have to compute the start index of
132  * each root
133  */
134  t = 0;
135  for(k=0;k<n;k++) {
136  i = invp[k];
137  if(father[i] == i) {
138  /* This is a root */
139  j = T[i];
140  T[i] += t;
141  t += j;
142  }
143  }
144
145 #if defined(PASTIX_DEBUG_ORDERING)
146  for(i=0;i<n;i++)
147  assert(T[i] <= T[father[i]]);
148 #endif
149
150  for(k=n-1;k>=0;k--)
151  {
152  i = invp[k];
153  perm[i] = T[father[i]]; /* We MUST HAVE father[i] == i for a root ! */
154  T[father[i]] -= T[i];
155  T[i] = perm[i]-1;
156  assert(perm[father[i]] >= perm[i]);
157  }
158
159  /*
160  * We need to retrieve 1 for the C numbering compliance
161  */
162  for(i=0;i<n;i++)
163  perm[i]--;
164
165 #if defined(PASTIX_DEBUG_ORDERING)
166  /* Check the permutation vector */
167  for(i=0;i<n;i++)
168  {
169  assert(perm[i] >= 0);
170  assert(perm[i] < n);
171  }
172
173  memset( invp, 0, sizeof(pastix_int_t) * n );
174  for(i=0;i<n;i++) {
175  invp[perm[i]]++;
176  }
177
178  k = 0;
179  for(i=0;i<n;i++) {
180  if(invp[i] != 1) { k++; }
181  }
182  if ( k > 0 ) {
183  pastix_print_error( "Number of errors in perm vector in postorder %ld", (long)k );
184  }
185  assert(k==0);
186 #endif
187
188  /*
189  * Compute the invp vector
190  */
191  for(i=0; i<n; i++) {
192  invp[perm[i]] = i;
193  }
194 }
195
196 /**
197  *******************************************************************************
198  *
199  * @ingroup order_dev
200  *
201  * @brief Compute the elimination tree of a matrix A
202  * (without computing the symbolic factorization) associated with a reordering
203  * of the matrix.
204  *
205  *******************************************************************************
206  *
207  * @param[in] n
208  * The number of vertices in the given CSC
209  *
210  * @param[in] ia
211  * Array of size n+1.
212  * The column pointer tabular for the given CSC. (In C numbering)
213  *
214  * @param[in] ja
215  * Array of size ia[n].
216  * The list of edges in the given CSC (In C numbering)
217  *
218  * @param[in] perm
219  * Array of size n. The permutation vector.
220  *
221  * @param[in] invp
222  * Array of size n. The inverse permutation vector.
223  *
224  * @param[inout] father
225  * Array of size n.
226  * On entry, an allocated array of size n.
227  * On exit, father[i] = father of ith node on the eliminination
228  * tree. If node i is a root then father[i] = i.
229  *
230  *******************************************************************************/
231 static inline void
233  const pastix_int_t *ia,
234  const pastix_int_t *ja,
235  const pastix_int_t *perm,
236  const pastix_int_t *invp,
237  pastix_int_t *father)
238 {
239  pastix_int_t i, j, k;
240  pastix_int_t node;
241  pastix_int_t vroot;
242
243  /* Optim */
244  pastix_int_t flag, ind;
245  pastix_int_t *jrev = NULL;
246  pastix_int_t *jf = NULL;
247
248  MALLOC_INTERN(jrev, n, pastix_int_t);
249  for(i=0;i<n;i++)
250  jrev[i] = -1;
251
252  MALLOC_INTERN(jf, n, pastix_int_t);
253  memset( jf, 0, sizeof(pastix_int_t) * n );
254
255  for(i=0;i<n;i++)
256  father[i] = -1;
257
258  for(i=0;i<n;i++)
259  {
260  ind = 0;
261  node = invp[i];
262  for(j=ia[node];j<ia[node+1];j++)
263  {
264
265  k = ja[j];
266  if(perm[k] < perm[node])
267  {
268  flag = 1;
269  vroot = k;
270  while(father[vroot] != -1 && father[vroot] != node)
271  {
272  if(jrev[vroot] >= 0)
273  {
274  flag = 0;
275  break;
276  }
277  jrev[vroot] = ind;
278  jf[ind] = vroot;
279  ind++;
280
281  vroot = father[vroot];
282  }
283  if(flag == 1)
284  father[vroot] = node;
285  }
286  }
287  /* reinit jrev */
288  for(j=0;j<ind;j++)
289  jrev[jf[j]]=-1;
290  }
291
292  memFree_null(jrev);
293  memFree_null(jf);
294
295  for(i=0;i<n;i++)
296  if(father[i] == -1)
297  father[i] = i;
298
299 #if defined(PASTIX_DEBUG_ORDERING)
300  /* Check to see if a father has a lower rank in the permutation array than one of its sons */
301  for(i=0;i<n;i++)
302  {
303  if(perm[i] > perm[father[i]])
304  {
305  fprintf(stderr, "Node %ld perm=%ld Father %ld perm=%ld \n",
306  (long)i, (long)perm[i], (long)father[i], (long)perm[father[i]]);
307  assert(perm[i] <= perm[father[i]]);
308  }
309  }
310 #endif
311 }
312
313
314 /**
315  *******************************************************************************
316  *
317  * @ingroup pastix_order
318  *
319  * @brief Computes the set of supernodes for a given permutation.
320  *
321  * The permutation of the matrix given on entry is modified to obtain a
322  * postorder of the elimination tree: this does not affect the fill-in properties
323  * of the initial ordering.
324  *
325  * WARNING: The matrix pattern is assumed to be symmetric.
326  *
327  * NOTE: This function can take on entry the lower triangular part or the whole
328  * matrix A.
329  *
330  *******************************************************************************
331  *
332  * @param[in] graph
333  * The graph associated with the order structuree on which we need to
334  * find the supernodes.
335  *
336  * @param[inout] ordeptr
337  * Pointer to a pastix_order_t structure, that will be further initialized by
338  * the routine.
339  *
340  * On entry:
341  * - orderptr->permtab: the original permutation vector for the
342  * elimination tree.
343  * - orderptr->permtab: the original inverse permutation vector for the
344  * elimination tree.
345  *
346  * On exit:
347  * - orderptr->cblknbr: The number of supernodes found.
348  * - orderptr->rangtab: Contains the first element of each
349  * supernode. rangtab[i] is the first element of the ith
350  * supernode.
351  * - orderptr->permtab: the permutation vector for the
352  * postorder of the nodes in the elimination tree.
353  * - orderptr->permtab: the inverse permutation vector for the
354  * postorder of the nodes in the elimination tree.
355  * - orderptr->treetab: treetab[i] is the number of the father of
356  * supernode i in the supernodal elimination tree.
357  *
358  *******************************************************************************/
359 void
360 orderFindSupernodes( const pastix_graph_t *graph,
361  pastix_order_t * const ordeptr )
362 {
363  pastix_int_t *father = NULL; /* father[i] is the father of node i in he elimination tree of A */
364  pastix_int_t *T; /* T[j] is the number of node in the subtree rooted in node j in
365  the elimination tree of A */
366  pastix_int_t *S = NULL; /* S[i] is the number of sons for node i in the elimination tree */
367  pastix_int_t *isleaf = NULL;
368  pastix_int_t *prev_rownz = NULL;
369  pastix_int_t *treetab = NULL;
370  pastix_int_t i, j, k;
371  pastix_int_t pi, pj;
372
373  pastix_int_t snodenbr;
374  pastix_int_t *perm = ordeptr->permtab;
375  pastix_int_t *invp = ordeptr->peritab;
376  pastix_int_t n = graph->n;
377  pastix_int_t *ia = graph->colptr;
378  pastix_int_t *ja = graph->rowptr;
379
380  assert( graph->colptr[0] == 0 );
381
382  /* Free the rangtab/treetab array if existing */
383  if ( ordeptr->rangtab != NULL ) {
384  memFree_null( ordeptr->rangtab );
385  }
386  if ( ordeptr->treetab != NULL ) {
387  memFree_null( ordeptr->treetab );
388  }
389
390  MALLOC_INTERN(S, n, pastix_int_t);
391  MALLOC_INTERN(father, n, pastix_int_t);
392
393 #if defined(PASTIX_DEBUG_ORDERING)
394  /* Check that the permutation vector is 0 based */
395  for(i=0;i<n;i++) {
396  assert(perm[i] >= 0);
397  assert(perm[i] < n);
398  assert(invp[i] >= 0);
399  assert(invp[i] < n);
400  }
401
402  memset( S, 0, sizeof(pastix_int_t) * n );
403  for(i=0;i<n;i++)
404  S[perm[i]]++;
405
406  k = 0;
407  for(i=0;i<n;i++) {
408  if(S[i] != 1) {
409  k++;
410  }
411  }
412
413  if(k>0) {
414  pastix_print_error( "perm array is not valid, number of error = %ld", (long)k );
415  }
416  assert(k==0);
417 #endif
418
419  /*
420  * Compute the elimination tree of A
421  */
422  compute_elimination_tree(n, ia, ja, perm, invp, father);
423
424  /*
425  * Compute the postorder of the elimination tree
426  * Warning: This operation modifies perm and invp
427  */
428  MALLOC_INTERN(T, n+1, pastix_int_t);
429  compute_post_order(n, father, perm, invp, T);
430
431  /*
432  * Compute the number of descendant of each node i in the elimination tree
433  */
434  compute_subtree_size(n, father, invp, T);
435
436  MALLOC_INTERN(isleaf, n, pastix_int_t);
437  MALLOC_INTERN(prev_rownz, n, pastix_int_t);
438  memset(isleaf, 0, sizeof(pastix_int_t) * n );
439  memset(prev_rownz, 0, sizeof(pastix_int_t) * n );
440
441  for(j=0;j<n;j++)
442  {
443  pj = invp[j];
444  for(i=ia[pj];i<ia[pj+1];i++)
445  {
446  pi = perm[ja[i]];
447  if(pi > j)
448  {
449  k = prev_rownz[pi];
450  if(k < j - T[pj]+1 ) {
451  isleaf[j] = 1;
452  }
453
454  prev_rownz[pi] = j;
455  }
456  }
457  }
458  memFree(prev_rownz);
459
460  /*
461  * Compute the number of sons of each node in the elimination tree
462  * The snodetab/rangtab is computed in the workspace T.
463  */
464  memset( S, 0, sizeof(pastix_int_t) * n );
465  for(i=0;i<n;i++) {
466  if(father[i] != i) {
467  S[father[i]]++;
468  }
469  }
470
471  for(i=0;i<n;i++) {
472  if(S[i] != 1) {
473  isleaf[perm[i]] = 1;
474  }
475  }
476
477  snodenbr = 0;
478  for(i=0;i<n;i++) {
479  if(isleaf[i] == 1)
480  {
481  T[snodenbr] = i;
482  snodenbr++;
483  }
484  }
485  T[snodenbr] = n;
486
487  memFree(isleaf);
488
489  /*
490  * If the treetab is required, we compute it before to free the local data.
491  */
492  assert( snodenbr > 0 );
493  MALLOC_INTERN( ordeptr->treetab, snodenbr, pastix_int_t );
494  treetab = ordeptr->treetab;
495  {
497
498  /* Node to supernode conversion vector */
499  for(i=0;i<snodenbr;i++) {
500  for(j=T[i];j<T[i+1];j++) {
501  S[j] = i;
502  }
503  }
504
505  /* Fill the treetab info */
506  for(i=0;i<snodenbr;i++)
507  {
508  k=snodenbr;
509  for(j=T[i];j<T[i+1];j++)
510  {
514  }
515  }
516  treetab[i] = k;
517  if(k == snodenbr)
518  {
519  treetab[i] = -1; /* This is a root */
520  }
521  assert((treetab[i] == -1) || (treetab[i] >= i));
522  }
523  }
524
525  memFree(father);
526  memFree(S);
527
528  /*
529  * Save result to ordeptr (switch former rangtab and new one allocated to
530  * restricted array)
531  */
532  ordeptr->cblknbr = snodenbr;
533  MALLOC_INTERN( ordeptr->rangtab, snodenbr+1, pastix_int_t );
534  memcpy( ordeptr->rangtab, T, (snodenbr+1)*sizeof(pastix_int_t) );
535  memFree(T);
536
537  /* Check that the order is 0 based */
538  assert(ordeptr->rangtab[0] == 0);
539 }
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
static void compute_elimination_tree(pastix_int_t n, const pastix_int_t *ia, const pastix_int_t *ja, const pastix_int_t *perm, const pastix_int_t *invp, pastix_int_t *father)
Compute the elimination tree of a matrix A (without computing the symbolic factorization) associated ...
static void compute_post_order(pastix_int_t n, const pastix_int_t *father, pastix_int_t *perm, pastix_int_t *invp, pastix_int_t *T)
Computes the post order of the elimination tree given on entry.
static void compute_subtree_size(pastix_int_t n, const pastix_int_t *father, const pastix_int_t *iperm, pastix_int_t *T)
Computes the size of each subtree.
pastix_int_t * treetab
Definition: order.h:54
pastix_int_t * permtab
Definition: order.h:51
pastix_int_t * peritab
Definition: order.h:52
pastix_int_t cblknbr
Definition: order.h:50
pastix_int_t * rangtab
Definition: order.h:53
void orderFindSupernodes(const pastix_graph_t *graph, pastix_order_t *const ordeptr)
Computes the set of supernodes for a given permutation.
Order structure.
Definition: order.h:47