1 /*****************************************************************************
2 * *
3 * Sparse-graph-specific auxiliary source file for version 2.7 of nauty. *
4 * *
5 * Copyright (2004-2016) Brendan McKay. All rights reserved. *
6 * Subject to waivers and disclaimers in nauty.h. *
7 * *
8 * CHANGE HISTORY *
9 * 26-Oct-04 : initial creation *
10 * 23-Nov-06 : dispatch uses targetcell_sg, not bestcell_sg *
11 * 8-Dec-06 : add adjacencies_sg() *
12 * 10-Nov-09 : remove types shortish and permutation *
13 * 14-Nov-09 : added copy_sg() *
14 * 11-May-10 : use sorttemplates.c for sorting procedures *
15 * 19-May-10 : add two *_tr procedures for traces. *
16 * 21-May-10 : fixes for nde,v fields becoming size_t *
17 * 23-May-10 : add sparsenauty() *
18 * 15-Jan-12 : add TLS_ATTR attributes *
19 * 17-Dec-15 : extend sortlists_sg() to sort weights *
20 * : add weights to copy_sg() and updatecan_sg() *
21 * 11-Mar-16 : add cleanup_sg(). This can be used in the cleanup *
22 * field of the dispatch vector to sort the lists of the *
23 * canonical graph, but isn't there by default. *
24 * 15-Oct-19 : fix static declaration of snwork[] *
25 * *
26 *****************************************************************************/
27
28 #define TMP
29
30 /* #define ONE_WORD_SETS not sure about this! See notes.txt. */
31 #include "nausparse.h"
32
33 /* macros for hash-codes: */
34 #define MASH(l,i) ((((l) ^ 065435) + (i)) & 077777)
35 /* : expression whose long value depends only on long l and int/long i.
36 Anything goes, preferably non-commutative. */
37
38 #define CLEANUP(l) ((int)((l) % 077777))
39 /* : expression whose value depends on long l and is less than 077777
40 when converted to int then short. Anything goes. */
41
42 #if MAXM==1
43 #define M 1
44 #else
45 #define M m
46 #endif
47
48 #define ACCUM(x,y) x = (((x) + (y)) & 077777)
49
50 static const int fuzz1[] = {037541,061532,005257,026416};
51 static const int fuzz2[] = {006532,070236,035523,062437};
52
53 #define FUZZ1(x) ((x) ^ fuzz1[(x)&3])
54 #define FUZZ2(x) ((x) ^ fuzz2[(x)&3])
55
56 /* aproto: header new_nauty_protos.h */
57
58 dispatchvec dispatch_sparse =
59 {isautom_sg,testcanlab_sg,updatecan_sg,refine_sg,refine_sg,cheapautom_sg,
60 targetcell_sg,nausparse_freedyn,nausparse_check,init_sg,NULL};
61
62 #if !MAXN
63 DYNALLSTAT(short,vmark1,vmark1_sz);
64 DYNALLSTAT(short,vmark2,vmark2_sz);
65 DYNALLSTAT(int,work1,work1_sz);
66 DYNALLSTAT(int,work2,work2_sz);
67 DYNALLSTAT(int,work3,work3_sz);
68 DYNALLSTAT(int,work4,work4_sz);
69 DYNALLSTAT(set,snwork,snwork_sz);
70 #else
71 static TLS_ATTR short vmark1[MAXN];
72 static TLS_ATTR short vmark2[MAXN];
73 static TLS_ATTR int work1[MAXN];
74 static TLS_ATTR int work2[MAXN];
75 static TLS_ATTR int work3[MAXN];
76 static TLS_ATTR int work4[MAXN];
77 static TLS_ATTR set snwork[2*500*MAXM];
78 #endif
79
80 static TLS_ATTR short vmark1_val = 32000;
81 #define MARK1(i) vmark1[i] = vmark1_val
82 #define UNMARK1(i) vmark1[i] = 0
83 #define ISMARKED1(i) (vmark1[i] == vmark1_val)
84 #define ISNOTMARKED1(i) (vmark1[i] != vmark1_val)
85
86 static TLS_ATTR short vmark2_val = 32000;
87 #define MARK2(i) vmark2[i] = vmark2_val
88 #define UNMARK2(i) vmark2[i] = 0
89 #define ISMARKED2(i) (vmark2[i] == vmark2_val)
90 #define ISNOTMARKED2(i) (vmark2[i] != vmark2_val)
91
92 #if !MAXN
93 #define RESETMARKS1 {if (vmark1_val++ >= 32000) \
94 {size_t ij; for (ij=0;ij<vmark1_sz;++ij) vmark1[ij]=0; vmark1_val=1;}}
95 #define PREPAREMARKS1(nn) preparemarks1(nn)
96 #define RESETMARKS2 {if (vmark2_val++ >= 32000) \
97 {size_t ij; for (ij=0;ij<vmark2_sz;++ij) vmark2[ij]=0; vmark2_val=1;}}
98 #define PREPAREMARKS2(nn) preparemarks2(nn)
99 #else
100 #define RESETMARKS1 {if (vmark1_val++ >= 32000) \
101 {size_t ij; for (ij=0;ij<MAXN;++ij) vmark1[ij]=0; vmark1_val=1;}}
102 #define PREPAREMARKS1(nn)
103 #define RESETMARKS2 {if (vmark2_val++ >= 32000) \
104 {size_t ij; for (ij=0;ij<MAXN;++ij) vmark2[ij]=0; vmark2_val=1;}}
105 #define PREPAREMARKS2(nn)
106 #endif
107
108 /*****************************************************************************
109 * *
110 * preparemarks1(N) and preparemarks2(N) *
111 * make vmarks array large enough to mark 0..N-1 and such that *
112 * the next RESETMARKS command will work correctly *
113 * *
114 *****************************************************************************/
115
116 #if !MAXN
117 static void
preparemarks1(size_t nn)118 preparemarks1(size_t nn)
119 {
120 size_t oldsize;
121 short *oldpos;
122
123 oldsize = vmark1_sz;
124 oldpos = vmark1;
125 DYNALLOC1(short,vmark1,vmark1_sz,nn,"preparemarks");
126 if (vmark1_sz != oldsize || vmark1 != oldpos) vmark1_val = 32000;
127 }
128 #endif
129
130 #if !MAXN
131 static void
preparemarks2(size_t nn)132 preparemarks2(size_t nn)
133 {
134 size_t oldsize;
135 short *oldpos;
136
137 oldsize = vmark2_sz;
138 oldpos = vmark2;
139 DYNALLOC1(short,vmark2,vmark2_sz,nn,"preparemarks");
140 if (vmark2_sz != oldsize || vmark2 != oldpos) vmark2_val = 32000;
141 }
142 #endif
143
144 /*****************************************************************************
145 * *
146 * isautom_sg(g,perm,digraph,m,n) = TRUE iff perm is an automorphism of g *
147 * (i.e., g^perm = g). Symmetry is assumed unless digraph = TRUE. *
148 * *
149 *****************************************************************************/
150
151 boolean
isautom_sg(graph * g,int * p,boolean digraph,int m,int n)152 isautom_sg(graph *g, int *p, boolean digraph, int m, int n)
153 {
154 int *d,*e;
155 size_t *v;
156 int i,pi,di;
157 size_t vi,vpi,j;
158
159 SG_VDE(g,v,d,e);
160
161 PREPAREMARKS1(n);
162
163 for (i = 0; i < n; ++i)
164 if (p[i] != i || digraph)
165 {
166 pi = p[i];
167 di = d[i];
168 if (d[pi] != di) return FALSE;
169
170 vi = v[i];
171 vpi = v[pi];
172 RESETMARKS1;
173 for (j = 0; j < di; ++j) MARK1(p[e[vi+j]]);
174 for (j = 0; j < di; ++j) if (ISNOTMARKED1(e[vpi+j])) return FALSE;
175 }
176
177 return TRUE;
178 }
179
180 /*****************************************************************************
181 * *
182 * aresame_sg(g1,g2) = TRUE iff g1 and g2 are identical as labelled digraphs *
183 * *
184 *****************************************************************************/
185
186 boolean
aresame_sg(sparsegraph * g1,sparsegraph * g2)187 aresame_sg(sparsegraph *g1, sparsegraph *g2)
188 {
189 int *d1,*e1;
190 int *d2,*e2;
191 int n,i,di;
192 size_t vi,*v1,*v2,j;
193
194 n = g1->nv;
195 if (g2->nv != n || g2->nde != g1->nde) return FALSE;
196
197 SG_VDE(g1,v1,d1,e1);
198 SG_VDE(g2,v2,d2,e2);
199
200 PREPAREMARKS1(n);
201
202 for (i = 0; i < n; ++i)
203 {
204 di = d1[i];
205 if (d2[i] != di) return FALSE;
206
207 vi = v1[i];
208 RESETMARKS1;
209 for (j = 0; j < di; ++j) MARK1(e1[vi+j]);
210 vi = v2[i];
211 for (j = 0; j < di; ++j) if (ISNOTMARKED1(e2[vi+j])) return FALSE;
212 }
213
214 return TRUE;
215 }
216
217 /*****************************************************************************
218 * *
219 * testcanlab_sg(g,canong,lab,samerows,m,n) compares g^lab to canong, *
220 * using an ordering which is immaterial since it's only used here. The *
221 * value returned is -1,0,1 if g^lab <,=,> canong. *samerows is set to *
222 * the number of rows (0..n) of canong which are the same as those of g^lab. *
223 * *
224 *****************************************************************************/
225
226 int
testcanlab_sg(graph * g,graph * canong,int * lab,int * samerows,int m,int n)227 testcanlab_sg(graph *g, graph *canong, int *lab, int *samerows, int m, int n)
228 {
229 int *d,*e;
230 int *cd,*ce;
231 int i,k,di,dli;
232 size_t j,vi,vli,*v,*cv;
233 int mina;
234
235 SG_VDE(g,v,d,e);
236 SG_VDE(canong,cv,cd,ce);
237
238 #if !MAXN
239 DYNALLOC1(int,work1,work1_sz,n,"testcanlab_sg");
240 #endif
241 #define INVLAB work1
242
243 PREPAREMARKS1(n);
244
245 for (i = 0; i < n; ++i) INVLAB[lab[i]] = i;
246
247 for (i = 0; i < n; ++i)
248 {
249 /* compare g[lab[i]]^INVLAB to canong[i] */
250 vi = cv[i];
251 di = cd[i];
252 vli = v[lab[i]];
253 dli = d[lab[i]];
254
255 if (di != dli)
256 {
257 *samerows = i;
258 if (di < dli) return -1;
259 return 1;
260 }
261
262 RESETMARKS1;
263 mina = n;
264 for (j = 0; j < di; ++j) MARK1(ce[vi+j]);
265 for (j = 0; j < di; ++j)
266 {
267 k = INVLAB[e[vli+j]];
268 if (ISMARKED1(k)) UNMARK1(k);
269 else if (k < mina) mina = k;
270 }
271 if (mina != n)
272 {
273 *samerows = i;
274 for (j = 0; j < di; ++j)
275 {
276 k = ce[vi+j];
277 if (ISMARKED1(k) && k < mina) return -1;
278 }
279 return 1;
280 }
281 }
282
283 *samerows = n;
284 return 0;
285 }
286
287 /*****************************************************************************
288 * *
289 * updatecan_sg(g,canong,lab,samerows,m,n) sets canong = g^lab, assuming *
290 * the first samerows vertices of canong are ok already. Also assumes *
291 * contiguity and ample space in canong. *
292 * *
293 *****************************************************************************/
294
295 void
updatecan_sg(graph * g,graph * canong,int * lab,int samerows,int m,int n)296 updatecan_sg(graph *g, graph *canong, int *lab, int samerows, int m, int n)
297 {
298 int *d,*e;
299 int *cd,*ce;
300 int i,dli;
301 size_t *v,*cv,vli,j,k;
302 sg_weight *wt,*cwt;
303
304 SWG_VDE(g,v,d,e,wt);
305 SWG_VDE(canong,cv,cd,ce,cwt);
306
307 #if !MAXN
308 DYNALLOC1(int,work1,work1_sz,n,"testcanlab_sg");
309 #endif
310 #define INVLAB work1
311
312 ((sparsegraph*)canong)->nv = n;
313 ((sparsegraph*)canong)->nde = ((sparsegraph*)g)->nde;
314
315 for (i = 0; i < n; ++i) INVLAB[lab[i]] = i;
316
317 if (samerows == 0) k = 0;
318 else k = cv[samerows-1]+cd[samerows-1];
319
320 for (i = samerows; i < n; ++i)
321 {
322 cv[i] = k;
323 cd[i] = dli = d[lab[i]];
324 vli = v[lab[i]];
325 if (wt)
326 {
327 for (j = 0; j < dli; ++j)
328 {
329 ce[k] = INVLAB[e[vli+j]];
330 cwt[k] = wt[vli+j];
331 ++k;
332 }
333 }
334 else
335 for (j = 0; j < dli; ++j) ce[k++] = INVLAB[e[vli+j]];
336 }
337 }
338
339 /*****************************************************************************
340 * *
341 * comparelab_tr(g,lab1,invlab1,lab2,invlab2,cls,col) compares *
342 * g^lab1 to g^lab2 and returns -1,0,1 according to the comparison. *
343 * invlab1[] and invlab2[] are assumed to hold inverses of lab1,lab2. *
344 * *
345 *****************************************************************************/
346
347 int
comparelab_tr(sparsegraph * g,int * lab1,int * invlab1,int * lab2,int * invlab2,int * cls,int * col)348 comparelab_tr(sparsegraph *g,
349 int *lab1, int *invlab1, int *lab2, int *invlab2, int *cls, int *col)
350 {
351 int d1,*e1,d2,*e2;
352 int i,j,k,n,c,end;
353 int mina;
354
355 n = g->nv;
356 PREPAREMARKS1(n);
357
358 for (c=0; c<n; c+=cls[c])
359 {
360 if (cls[c] == 1)
361 {
362 end = c+cls[c];
363 for (i = c; i < end; ++i)
364 {
365 e1 = g->e + g->v[lab1[i]];
366 d1 = g->d[lab1[i]];
367 e2 = g->e + g->v[lab2[i]];
368 d2 = g->d[lab2[i]];
369 if (d1 < d2) return -1;
370 else if (d1 > d2) return 1;
371
372 RESETMARKS1;
373 mina = n;
374 for (j = 0; j < d1; ++j) MARK1(col[invlab1[e1[j]]]);
375
376 for (j = 0; j < d1; ++j)
377 {
378 k = col[invlab2[e2[j]]];
379 if (ISMARKED1(k)) UNMARK1(k);
380 else if (k < mina) mina = k;
381 }
382 if (mina != n)
383 {
384 for (j = 0; j < d1; ++j)
385 {
386 k = col[invlab1[e1[j]]];
387 if (ISMARKED1(k) && k < mina) return -1;
388 }
389 return 1;
390 }
391 }
392 }
393 }
394
395 return 0;
396 }
397
398 /*****************************************************************************
399 * *
400 * testcanlab_tr(g,canong,lab,invlab,samerows) compares g^lab to canong, *
401 * using an ordering which is immaterial since it's only used here. The *
402 * value returned is -1,0,1 if g^lab <,=,> canong. *samerows is set to *
403 * the number of rows (0..n) of canong which are the same as those of g^lab. *
404 * invlab[] is assumed to hold the inverse of lab[] *
405 * *
406 *****************************************************************************/
407
408 int
testcanlab_tr(sparsegraph * g,sparsegraph * canong,int * lab,int * invlab,int * samerows)409 testcanlab_tr(sparsegraph *g, sparsegraph *canong,
410 int *lab, int *invlab, int *samerows)
411 {
412 int *d,*e;
413 int *cd,*ce;
414 int i,di,dli;
415 int k,n;
416 size_t *v,*cv,vi,vli,j;
417 int mina;
418
419 SG_VDE(g,v,d,e);
420 SG_VDE(canong,cv,cd,ce);
421 n = g->nv;
422
423 PREPAREMARKS1(n);
424
425 for (i = 0; i < n; ++i)
426 {
427 /* compare g[lab[i]]^invlab to canong[i] */
428 vi = cv[i];
429 di = cd[i];
430 vli = v[lab[i]];
431 dli = d[lab[i]];
432
433 if (di != dli)
434 {
435 *samerows = i;
436 if (di < dli) return -1;
437 return 1;
438 }
439
440 RESETMARKS1;
441 mina = n;
442 for (j = 0; j < di; ++j) MARK1(ce[vi+j]);
443
444 for (j = 0; j < di; ++j)
445 {
446 k = invlab[e[vli+j]];
447 if (ISMARKED1(k)) UNMARK1(k);
448 else if (k < mina) mina = k;
449 }
450 if (mina != n)
451 {
452 *samerows = i;
453 for (j = 0; j < di; ++j)
454 {
455 k = ce[vi+j];
456 if (ISMARKED1(k) && k < mina) return -1;
457 }
458 return 1;
459 }
460 }
461
462 *samerows = n;
463 return 0;
464 }
465
466 /*****************************************************************************
467 * *
468 * updatecan_tr(g,canong,lab,invlab,samerows) sets canong = g^lab, *
469 * assuming the first samerows vertices of canong are ok already. *
470 * Also assumes contiguity and ample space in canong. *
471 * Assumes invlab[] holds the inverse of lab[] *
472 * *
473 *****************************************************************************/
474
475 void
updatecan_tr(sparsegraph * g,sparsegraph * canong,int * lab,int * invlab,int samerows)476 updatecan_tr(sparsegraph *g, sparsegraph *canong,
477 int *lab, int *invlab, int samerows)
478 {
479 int *d,*e;
480 int *cd,*ce;
481 int i,dli,n;
482 size_t *v,*cv,vli,j,k;
483
484 SG_VDE(g,v,d,e);
485 SG_VDE(canong,cv,cd,ce);
486 n = g->nv;
487
488 PREPAREMARKS1(n);
489
490 canong->nv = n;
491 canong->nde = g->nde;
492
493 if (samerows == 0) k = 0;
494 else k = cv[samerows-1]+cd[samerows-1];
495
496 for (i = samerows; i < n; ++i)
497 {
498 cv[i] = k;
499 cd[i] = dli = d[lab[i]];
500 vli = v[lab[i]];
501 for (j = 0; j < dli; ++j) ce[k++] = invlab[e[vli+j]];
502 }
503 }
504
505 #define SORT_OF_SORT 3
506 #define SORT_NAME sortindirect
507 #define SORT_TYPE1 int
508 #define SORT_TYPE2 int
509 #include "sorttemplates.c"
510
511 #define SORT_OF_SORT 1
512 #define SORT_NAME sortints
513 #define SORT_TYPE1 int
514 #include "sorttemplates.c"
515
516 #define SORT_OF_SORT 2
517 #define SORT_NAME sortweights
518 #define SORT_TYPE1 int
519 #define SORT_TYPE2 sg_weight
520 #include "sorttemplates.c"
521
522 /*****************************************************************************
523 * *
524 * init_sg(graph *gin, graph **gout, graph *hin, graph **hout, *
525 * int *lab, int *ptn, set *active, optionblk *options, *
526 * int *status, int m, int n) *
527 * Initialise routine for dispatch vector. This one just makes sure *
528 * that *hin has enough space and sets fields for n=0. *
529 * *
530 *****************************************************************************/
531
532 void
init_sg(graph * gin,graph ** gout,graph * hin,graph ** hout,int * lab,int * ptn,set * active,struct optionstruct * options,int * status,int m,int n)533 init_sg(graph *gin, graph **gout, graph *hin, graph **hout, int *lab,
534 int *ptn, set *active, struct optionstruct *options, int *status,
535 int m, int n)
536 {
537 sparsegraph *sg,*sh;
538
539 if (options->getcanon)
540 {
541 sg = (sparsegraph*)gin;
542 sh = (sparsegraph*)hin;
543 SG_ALLOC(*sh,sg->nv,sg->nde,"init_sg");
544 sh->nv = sg->nv;
545 sh->nde = sg->nde;
546 }
547 *status = 0;
548 }
549
550 /*****************************************************************************
551 * *
552 * cleanup_sg(graph *gin, graph **gout, graph *hin, graph **hout, *
553 * int *lab, int *ptn, optionblk *options, *
554 * statsblk *stats, int m, int n) *
555 * Cleanup routine for dispatch vector. This one sorts the adjacency *
556 * lists for the canonical labelling. *
557 * *
558 *****************************************************************************/
559
560 void
cleanup_sg(graph * gin,graph ** gout,graph * hin,graph ** hout,int * lab,int * ptn,optionblk * options,statsblk * stats,int m,int n)561 cleanup_sg(graph *gin, graph **gout, graph *hin, graph **hout, int *lab,
562 int *ptn, optionblk *options, statsblk *stats, int m, int n)
563 {
564 sparsegraph *sh;
565
566 if (options->getcanon
567 && (stats->errstatus == 0 || stats->errstatus == NAUABORTED))
568 {
569 sh = (sparsegraph*)hin;
570 sortlists_sg(sh);
571 }
572 }
573
574 /*****************************************************************************
575 * *
576 * distvals(sparsegraph *sg, int v0, int *dist, int n) sets dist[i] *
577 * to the distance from v0 to i, for each i, or to n if there is no such *
578 * distance. work4[] is used as a queue. *
579 * *
580 *****************************************************************************/
581
582 void
distvals(sparsegraph * g,int v0,int * dist,int n)583 distvals(sparsegraph *g, int v0, int *dist, int n)
584 {
585 int *d,*e;
586 int i,head,tail;
587 int di,k;
588 size_t *v,vi,j;
589
590 SG_VDE(g,v,d,e);
591 #if !MAXN
592 DYNALLOC1(int,work4,work4_sz,n,"distvals");
593 #endif
594 #define QUEUE work4
595
596 for (i = 0; i < n; ++i) dist[i] = n;
597
598 QUEUE[0] = v0;
599 dist[v0] = 0;
600
601 head = 0;
602 tail = 1;
603 while (tail < n && head < tail)
604 {
605 i = QUEUE[head++];
606 vi = v[i];
607 di = d[i];
608 for (j = 0; j < di; ++j)
609 {
610 k = e[vi+j];
611 if (dist[k] == n)
612 {
613 dist[k] = dist[i] + 1;
614 QUEUE[tail++] = k;
615 }
616 }
617 }
618 }
619
620 /*****************************************************************************
621 * *
622 * refine_sg(g,lab,ptn,level,numcells,count,active,code,m,n) performs a *
623 * refinement operation on the partition at the specified level of the *
624 * partition nest (lab,ptn). *numcells is assumed to contain the number of *
625 * cells on input, and is updated. The initial set of active cells (alpha *
626 * in the paper) is specified in the set active. Precisely, x is in active *
627 * iff the cell starting at index x in lab is active. *
628 * The resulting partition is equitable if active is correct (see the paper *
629 * and the Guide). *
630 * *code is set to a value which depends on the fine detail of the *
631 * algorithm, but which is independent of the labelling of the graph. *
632 * count is used for work space. *
633 * *
634 *****************************************************************************/
635
636 void
refine_sg(graph * g,int * lab,int * ptn,int level,int * numcells,int * count,set * active,int * code,int m,int n)637 refine_sg(graph *g, int *lab, int *ptn, int level, int *numcells,
638 int *count, set *active, int *code, int m, int n)
639 {
640 int i,j,k,l,v1,v2,v3,isplit;
641 int w1,w2,w3;
642 long longcode;
643 int *d,*e;
644 int size,bigsize,bigpos;
645 int nactive,hitcells;
646 int lj,di,splitv;
647 boolean trivsplit;
648 size_t *v,vi,ii;
649
650 SG_VDE(g,v,d,e);
651
652 #if !MAXN
653 DYNALLOC1(int,work1,work1_sz,n,"refine_sg");
654 DYNALLOC1(int,work2,work2_sz,n,"refine_sg");
655 DYNALLOC1(int,work3,work3_sz,n,"refine_sg");
656 DYNALLOC1(int,work4,work4_sz,n,"refine_sg");
657 #endif
658 #define CELLSTART work1
659 #define ACTIVE work2
660 #define HITS work3
661 #define HITCELL work4
662
663 PREPAREMARKS1(n);
664 PREPAREMARKS2(n);
665
666 longcode = *numcells;
667
668 /* Set ACTIVE[0..nactive-1] = queue of active cell starts */
669
670 nactive = 0;
671 for (i = -1; (i = nextelement(active,m,i)) >= 0;)
672 ACTIVE[nactive++] = i;
673
674 if (nactive == 0)
675 {
676 *code = CLEANUP(longcode);
677 return;
678 }
679
680 /* Set CELLSTART[i] = starting point in lab[] of nontrivial cell
681 containing i, or n if i is a singleton */
682
683 for (i = 0; i < n; )
684 {
685 /* Just here, i is a cell starting position */
686 if (ptn[i] <= level)
687 {
688 CELLSTART[lab[i]] = n;
689 ++i;
690 }
691 else
692 {
693 j = i;
694 do
695 {
696 CELLSTART[lab[i]] = j;
697 } while (ptn[i++] > level);
698 }
699 }
700
701 if (level <= 2 && nactive == 1 && ptn[ACTIVE[0]] <= level
702 && *numcells <= n/8)
703 {
704 isplit = ACTIVE[--nactive];
705 DELELEMENT(active,isplit);
706
707 distvals((sparsegraph*)g,lab[isplit],HITS,n);
708
709 for (v1 = 0; v1 < n; )
710 {
711 if (ptn[v1] <= level)
712 {
713 ++v1;
714 continue;
715 }
716
717 longcode = MASH(longcode,v1);
718 w1 = HITS[lab[v1]];
719
720 v2 = v1+1;
721 while (ptn[v2-1] > level && HITS[lab[v2]] == w1) ++v2;
722
723 if (ptn[v2-1] <= level)
724 {
725 v1 = v2;
726 continue;
727 }
728
729 w2 = NAUTY_INFINITY;
730 v3 = j = v2;
731
732 do
733 {
734 lj = lab[j];
735 w3 = HITS[lj];
736 if (w3 == w1)
737 {
738 lab[j] = lab[v3];
739 lab[v3] = lab[v2];
740 lab[v2] = lj;
741 ++v2;
742 ++v3;
743 }
744 else if (w3 == w2)
745 {
746 lab[j] = lab[v3];
747 lab[v3] = lj;
748 ++v3;
749 }
750 else if (w3 < w1)
751 {
752 lab[j] = lab[v2];
753 lab[v2] = lab[v1];
754 lab[v1] = lj;
755 v3 = v2 + 1;
756 v2 = v1 + 1;
757 w2 = w1;
758 w1 = w3;
759 }
760 else if (w3 < w2)
761 {
762 lab[j] = lab[v2];
763 lab[v2] = lj;
764 v3 = v2 + 1;
765 w2 = w3;
766 }
767 } while (ptn[j++] > level);
768
769 longcode = MASH(longcode,w2);
770 longcode = MASH(longcode,v2);
771 if (j != v2) /* At least two fragments
772 * v1..v2-1 = w1; v2..v3-1 = w2 */
773 {
774 if (v2 == v1+1)
775 CELLSTART[lab[v1]] = n;
776
777 if (v3 == v2+1)
778 CELLSTART[lab[v2]] = n;
779 else
780 for (k = v2; k < v3; ++k)
781 CELLSTART[lab[k]] = v2;
782 ++*numcells;
783 ptn[v2-1] = level;
784
785 if (j == v3)
786 {
787 /* Two fragments only */
788 if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1))
789 {
790 ADDELEMENT(active,v1);
791 ACTIVE[nactive++] = v1;
792 }
793 else
794 {
795 ADDELEMENT(active,v2);
796 ACTIVE[nactive++] = v2;
797 }
798 }
799 else
800 {
801 /* Extra fragments: v3..j-1 > w2 */
802 sortindirect(lab+v3,HITS,j-v3);
803 ACTIVE[nactive++] = v2;
804 ADDELEMENT(active,v2);
805 if (v2-v1 >= v3-v2)
806 {
807 bigpos = -1;
808 bigsize = v2-v1;
809 }
810 else
811 {
812 bigpos = nactive-1;
813 bigsize = v3-v2;
814 }
815 for (k = v3-1; k < j-1;)
816 {
817 ptn[k] = level;
818 longcode = MASH(longcode,k);
819 ++*numcells;
820 l = k+1;
821 ADDELEMENT(active,l);
822 ACTIVE[nactive++] = l;
823 w3 = HITS[lab[l]];
824 for (k = l; k < j-1
825 && HITS[lab[k+1]] == w3; ++k)
826 CELLSTART[lab[k+1]] = l;
827 size = k-l+1;
828 if (size == 1)
829 CELLSTART[lab[l]] = n;
830 else
831 {
832 CELLSTART[lab[l]] = l;
833 if (size > bigsize)
834 {
835 bigsize = size;
836 bigpos = nactive-1;
837 }
838 }
839 }
840
841 if (bigpos >= 0 && !ISELEMENT(active,v1))
842 {
843 longcode = MASH(longcode,bigpos);
844 DELELEMENT(active,ACTIVE[bigpos]);
845 ADDELEMENT(active,v1);
846 ACTIVE[bigpos] = v1;
847 }
848 }
849 }
850 v1 = j;
851 }
852 }
853
854 /* Iterate until complete */
855 while (nactive > 0 && *numcells < n)
856 {
857 for (i = 0; i < nactive && i < 10; ++i)
858 if (ptn[ACTIVE[i]] <= level) break;
859
860 if (i < nactive && i < 10)
861 {
862 trivsplit = TRUE;
863 isplit = ACTIVE[i];
864 ACTIVE[i] = ACTIVE[--nactive];
865 }
866 else
867 {
868 isplit = ACTIVE[--nactive];
869 trivsplit = ptn[isplit] <= level;
870 }
871
872 DELELEMENT(active,isplit);
873 longcode = MASH(longcode,isplit);
874
875 if (trivsplit)
876 {
877 RESETMARKS1;
878 RESETMARKS2;
879 hitcells = 0;
880 splitv = lab[isplit];
881 vi = v[splitv];
882 di = d[splitv];
883 for (ii = 0; ii < di; ++ii)
884 {
885 j = e[vi+ii];
886 MARK2(j);
887 k = CELLSTART[j];
888 if (k != n && ISNOTMARKED1(k))
889 {
890 MARK1(k);
891 HITCELL[hitcells++] = k;
892 }
893 }
894
895 if (hitcells > 1) sortints(HITCELL,hitcells);
896 longcode = MASH(longcode,hitcells);
897
898 /* divide cells according to which vertices are hit */
899
900 for (i = 0; i < hitcells; ++i)
901 {
902 j = v1 = v2 = HITCELL[i];
903 longcode = MASH(longcode,v2);
904 k = 0;
905 do
906 {
907 lj = lab[j];
908 if (ISMARKED2(lj))
909 HITS[k++] = lj;
910 else
911 lab[v2++] = lj;
912 } while (ptn[j++] > level);
913
914 longcode = MASH(longcode,k);
915 v3 = v2;
916 while (--k >= 0)
917 {
918 j = HITS[k];
919 CELLSTART[j] = v2;
920 lab[v3++] = j;
921 }
922
923 if (v2 != v3 && v2 != v1)
924 {
925 ++*numcells;
926 if (v2 == v1+1) CELLSTART[lab[v1]] = n;
927 if (v3 == v2+1) CELLSTART[lab[v2]] = n;
928 ptn[v2-1] = level;
929 longcode = MASH(longcode,v2);
930 if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1))
931 {
932 ADDELEMENT(active,v1);
933 ACTIVE[nactive++] = v1;
934 }
935 else
936 {
937 ADDELEMENT(active,v2);
938 ACTIVE[nactive++] = v2;
939 }
940 }
941 }
942 }
943 else /* non-trivial splitting */
944 {
945 /* isplit is the start of the splitting cell.
946 Set HITS[i] = hits of i for i in non-trivial cells,
947 HITCELL[0..hitcells-1] = starts of hit non-trivial cells */
948
949 RESETMARKS1;
950 hitcells = 0;
951 do
952 {
953 vi = v[lab[isplit]];
954 di = d[lab[isplit]];
955 for (ii = 0; ii < di; ++ii)
956 {
957 j = e[vi+ii];
958 k = CELLSTART[j];
959 if (k != n)
960 {
961 if (ISNOTMARKED1(k))
962 {
963 MARK1(k);
964 HITCELL[hitcells++] = k;
965 do
966 {
967 HITS[lab[k]] = 0;
968 } while (ptn[k++] > level);
969 }
970 ++HITS[j];
971 }
972 }
973 } while (ptn[isplit++] > level);
974
975 if (hitcells > 1) sortints(HITCELL,hitcells);
976
977 /* divide cells according to hit counts */
978
979 longcode = MASH(longcode,hitcells);
980 for (i = 0; i < hitcells; ++i)
981 {
982 v1 = HITCELL[i];
983 w1 = HITS[lab[v1]];
984 longcode = MASH(longcode,v1);
985
986 v2 = v1+1;
987 while (ptn[v2-1] > level && HITS[lab[v2]] == w1) ++v2;
988
989 if (ptn[v2-1] <= level) continue;
990 w2 = NAUTY_INFINITY;
991 v3 = j = v2;
992
993 do
994 {
995 lj = lab[j];
996 w3 = HITS[lj];
997 if (w3 == w1)
998 {
999 lab[j] = lab[v3];
1000 lab[v3] = lab[v2];
1001 lab[v2] = lj;
1002 ++v2;
1003 ++v3;
1004 }
1005 else if (w3 == w2)
1006 {
1007 lab[j] = lab[v3];
1008 lab[v3] = lj;
1009 ++v3;
1010 }
1011 else if (w3 < w1)
1012 {
1013 lab[j] = lab[v2];
1014 lab[v2] = lab[v1];
1015 lab[v1] = lj;
1016 v3 = v2 + 1;
1017 v2 = v1 + 1;
1018 w2 = w1;
1019 w1 = w3;
1020 }
1021 else if (w3 < w2)
1022 {
1023 lab[j] = lab[v2];
1024 lab[v2] = lj;
1025 v3 = v2 + 1;
1026 w2 = w3;
1027 }
1028 } while (ptn[j++] > level);
1029
1030 longcode = MASH(longcode,w1);
1031 longcode = MASH(longcode,v2);
1032 if (j != v2) /* At least two fragments
1033 * v1..v2-1 = w1; v2..v3-1 = w2 */
1034 {
1035 if (v2 == v1+1)
1036 CELLSTART[lab[v1]] = n;
1037
1038 if (v3 == v2+1)
1039 CELLSTART[lab[v2]] = n;
1040 else
1041 for (k = v2; k < v3; ++k)
1042 CELLSTART[lab[k]] = v2;
1043 ++*numcells;
1044 ptn[v2-1] = level;
1045
1046 if (j == v3)
1047 {
1048 /* Two fragments only */
1049 if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1))
1050 {
1051 ADDELEMENT(active,v1);
1052 ACTIVE[nactive++] = v1;
1053 }
1054 else
1055 {
1056 ADDELEMENT(active,v2);
1057 ACTIVE[nactive++] = v2;
1058 }
1059 }
1060 else
1061 {
1062 /* Extra fragments: v3..j-1 > w2 */
1063 longcode = MASH(longcode,v3);
1064 sortindirect(lab+v3,HITS,j-v3);
1065 ACTIVE[nactive++] = v2;
1066 ADDELEMENT(active,v2);
1067 if (v2-v1 >= v3-v2)
1068 {
1069 bigpos = -1;
1070 bigsize = v2-v1;
1071 }
1072 else
1073 {
1074 bigpos = nactive-1;
1075 bigsize = v3-v2;
1076 longcode = MASH(longcode,bigsize);
1077 }
1078 for (k = v3-1; k < j-1;)
1079 {
1080 ptn[k] = level;
1081 ++*numcells;
1082 l = k+1;
1083 ADDELEMENT(active,l);
1084 ACTIVE[nactive++] = l;
1085 w3 = HITS[lab[l]];
1086 longcode = MASH(longcode,w3);
1087 for (k = l; k < j-1
1088 && HITS[lab[k+1]] == w3; ++k)
1089 CELLSTART[lab[k+1]] = l;
1090 size = k-l+1;
1091 if (size == 1)
1092 CELLSTART[lab[l]] = n;
1093 else
1094 {
1095 CELLSTART[lab[l]] = l;
1096 if (size > bigsize)
1097 {
1098 bigsize = size;
1099 bigpos = nactive-1;
1100 }
1101 }
1102 }
1103
1104 if (bigpos >= 0 && !ISELEMENT(active,v1))
1105 {
1106 DELELEMENT(active,ACTIVE[bigpos]);
1107 ADDELEMENT(active,v1);
1108 ACTIVE[bigpos] = v1;
1109 }
1110 }
1111 }
1112 }
1113 }
1114 }
1115
1116 longcode = MASH(longcode,*numcells);
1117 *code = CLEANUP(longcode);
1118 }
1119
1120 /*****************************************************************************
1121 * *
1122 * cheapautom_sg(ptn,level,digraph,n) returns TRUE if the partition at the *
1123 * specified level in the partition nest (lab,ptn) {lab is not needed here} *
1124 * satisfies a simple sufficient condition for its cells to be the orbits of *
1125 * some subgroup of the automorphism group. Otherwise it returns FALSE. *
1126 * It always returns FALSE if digraph!=FALSE. *
1127 * *
1128 * nauty assumes that this function will always return TRUE for any *
1129 * partition finer than one for which it returns TRUE. *
1130 * *
1131 *****************************************************************************/
1132
1133 boolean
cheapautom_sg(int * ptn,int level,boolean digraph,int n)1134 cheapautom_sg(int *ptn, int level, boolean digraph, int n)
1135 {
1136 int i,k,nnt;
1137
1138 if (digraph) return FALSE;
1139
1140 k = n;
1141 nnt = 0;
1142 for (i = 0; i < n; ++i)
1143 {
1144 --k;
1145 if (ptn[i] > level)
1146 {
1147 ++nnt;
1148 while (ptn[++i] > level) {}
1149 }
1150 }
1151
1152 return (k <= nnt + 1 || k <= 4);
1153 }
1154
1155 /*****************************************************************************
1156 * *
1157 * bestcell_sg(g,lab,ptn,level,tc_level,m,n) returns the index in lab of *
1158 * the start of the "best non-singleton cell" for fixing. If there is no *
1159 * non-singleton cell it returns n. *
1160 * This implementation finds the first cell which is non-trivially joined *
1161 * to the greatest number of other cells, assuming equitability. *
1162 * This is not good for digraphs! *
1163 * *
1164 *****************************************************************************/
1165
1166 static int
bestcell_sg(graph * g,int * lab,int * ptn,int level,int tc_level,int m,int n)1167 bestcell_sg(graph *g, int *lab, int *ptn, int level,
1168 int tc_level, int m, int n)
1169 {
1170 int nnt;
1171 int *d,*e;
1172 int i,k,di;
1173 int *work1b;
1174 int maxcnt;
1175 size_t *v,vi,j;
1176
1177 SG_VDE(g,v,d,e);
1178
1179 #if !MAXN
1180 DYNALLOC1(int,work1,work1_sz,n,"bestcell_sg");
1181 DYNALLOC1(int,work2,work2_sz,n,"bestcell_sg");
1182 DYNALLOC1(int,work3,work3_sz,n,"bestcell_sg");
1183 DYNALLOC1(int,work4,work4_sz,n,"bestcell_sg");
1184 #endif
1185 work1b = work1 + (n/2);
1186 #define START work1
1187 #define SIZE work1b
1188 #define NNTCELL work2
1189 #define HITS work3
1190 #define COUNT work4
1191
1192 /* find non-singleton cells: put starts in START[0..nnt-1],
1193 sizes in SIZE[0..nnt-1].
1194 Also NNTCELL[i] = n if {i} is a singelton, else index of
1195 nontriv cell containing i. */
1196
1197 i = nnt = 0;
1198
1199 while (i < n)
1200 {
1201 if (ptn[i] > level)
1202 {
1203 START[nnt] = i;
1204 j = i;
1205 do
1206 NNTCELL[lab[j]] = nnt;
1207 while (ptn[j++] > level);
1208 SIZE[nnt] = j-i;
1209 ++nnt;
1210 i = j;
1211 }
1212 else
1213 {
1214 NNTCELL[lab[i]] = n;
1215 ++i;
1216 }
1217 }
1218
1219 if (nnt == 0) return n;
1220
1221 /* set COUNT[i] to # non-trivial neighbours of n.s. cell i */
1222
1223 for (i = 0; i < nnt; ++i) HITS[i] = COUNT[i] = 0;
1224
1225 for (i = 0; i < nnt; ++i)
1226 {
1227 vi = v[lab[START[i]]];
1228 di = d[lab[START[i]]];
1229
1230 for (j = 0; j < di; ++j)
1231 {
1232 k = NNTCELL[e[vi+j]];
1233 if (k != n) ++HITS[k];
1234 }
1235 for (j = 0; j < di; ++j)
1236 {
1237 k = NNTCELL[e[vi+j]];
1238 if (k != n)
1239 {
1240 if (HITS[k] > 0 && HITS[k] < SIZE[k]) ++COUNT[i];
1241 HITS[k] = 0;
1242 }
1243 }
1244 }
1245
1246 /* find first greatest bucket value */
1247
1248 j = 0;
1249 maxcnt = COUNT[0];
1250 for (i = 1; i < nnt; ++i)
1251 if (COUNT[i] > maxcnt)
1252 {
1253 j = i;
1254 maxcnt = COUNT[i];
1255 }
1256
1257 return (int)START[j];
1258 }
1259 /*****************************************************************************
1260 * *
1261 * targetcell_sg(g,lab,ptn,level,tc_level,digraph,hint,m,n) returns the *
1262 * index in lab of the next cell to split. *
1263 * hint is a suggestion for the answer, which is obeyed if it is valid. *
1264 * Otherwise we use bestcell() up to tc_level and the first non-trivial *
1265 * cell after that. *
1266 * *
1267 *****************************************************************************/
1268
1269 int
targetcell_sg(graph * g,int * lab,int * ptn,int level,int tc_level,boolean digraph,int hint,int m,int n)1270 targetcell_sg(graph *g, int *lab, int *ptn, int level, int tc_level,
1271 boolean digraph, int hint, int m, int n)
1272 {
1273 int i;
1274
1275 if (hint >= 0 && ptn[hint] > level &&
1276 (hint == 0 || ptn[hint-1] <= level))
1277 return hint;
1278 else if (level <= tc_level)
1279 return bestcell_sg(g,lab,ptn,level,tc_level,m,n);
1280 else
1281 {
1282 for (i = 0; i < n && ptn[i] <= level; ++i) {}
1283 return (i == n ? 0 : i);
1284 }
1285 }
1286
1287 /*****************************************************************************
1288 * *
1289 * sortlists_sg(g) sorts the adjacency lists into numerical order *
1290 * *
1291 *****************************************************************************/
1292
1293 void
sortlists_sg(sparsegraph * g)1294 sortlists_sg(sparsegraph *g)
1295 {
1296 int *d,*e;
1297 int n,i;
1298 size_t *v;
1299 sg_weight *wt;
1300
1301 SWG_VDE(g,v,d,e,wt);
1302 n = g->nv;
1303
1304 if (wt)
1305 {
1306 for (i = 0; i < n; ++i)
1307 if (d[i] > 1) sortweights(e+v[i],wt+v[i],d[i]);
1308 }
1309 else
1310 {
1311 for (i = 0; i < n; ++i)
1312 if (d[i] > 1) sortints(e+v[i],d[i]);
1313 }
1314 }
1315
1316 /*****************************************************************************
1317 * *
1318 * put_sg(f,sg,digraph,linelength) writes the sparse graph to file f using *
1319 * at most linelength characters per line. If digraph then all directed *
1320 * edges are written; else one v-w for w>=v is written. *
1321 * *
1322 *****************************************************************************/
1323
1324 void
put_sg(FILE * f,sparsegraph * sg,boolean digraph,int linelength)1325 put_sg(FILE *f, sparsegraph *sg, boolean digraph, int linelength)
1326 {
1327 int *d,*e;
1328 int n,di;
1329 int i,curlen,slen;
1330 size_t *v,vi,j;
1331 char s[12];
1332
1333 SG_VDE(sg,v,d,e);
1334 n = sg->nv;
1335
1336 for (i = 0; i < n; ++i)
1337 {
1338 vi = v[i];
1339 di = d[i];
1340 if (di == 0) continue;
1341 slen = itos(i+labelorg,s);
1342 putstring(f,s);
1343 putstring(f," :");
1344 curlen = slen + 2;
1345
1346 for (j = 0; j < di; ++j)
1347 {
1348 if (!digraph && e[vi+j] < i) continue;
1349 slen = itos(e[vi+j]+labelorg,s);
1350 if (linelength && curlen + slen + 1 >= linelength)
1351 {
1352 putstring(f,"\n ");
1353 curlen = 2;
1354 }
1355 PUTC(' ',f);
1356 putstring(f,s);
1357 curlen += slen + 1;
1358 }
1359 PUTC('\n',f);
1360 }
1361 }
1362
1363 /*****************************************************************************
1364 * *
1365 * sg_to_nauty(sg,g,reqm,&m) creates a nauty-format graph from a sparse *
1366 * graph. reqm is the required m value (computed from n if reqm=0), and *
1367 * m is the actual value used. g is dynamically generated if NULL is given. *
1368 * A pointer to g is returned. *
1369 * *
1370 *****************************************************************************/
1371
1372 graph*
sg_to_nauty(sparsegraph * sg,graph * g,int reqm,int * pm)1373 sg_to_nauty(sparsegraph *sg, graph *g, int reqm, int *pm)
1374 {
1375 int *d,*e;
1376 int m,n,i,di;
1377 size_t *v,vi,j;
1378 set *gi;
1379
1380 SG_VDE(sg,v,d,e);
1381 n = sg->nv;
1382 if (reqm != 0 && reqm*WORDSIZE < n)
1383 {
1384 fprintf(ERRFILE,"sg_to_nauty: reqm is impossible\n");
1385 exit(1);
1386 }
1387
1388 if (reqm != 0) m = reqm;
1389 else m = (n+WORDSIZE-1)/WORDSIZE;
1390
1391 *pm = m;
1392
1393 if (g == NULL)
1394 {
1395 if ((g = (graph*)ALLOCS(n,m*sizeof(graph))) == NULL)
1396 {
1397 fprintf(ERRFILE,"sg_to_nauty: malloc failed\n");
1398 exit(1);
1399 }
1400 }
1401
1402 for (i = 0, gi = g; i < n; ++i, gi += m)
1403 {
1404 vi = v[i];
1405 di = d[i];
1406 EMPTYSET(gi,m);
1407 for (j = 0; j < di; ++j) ADDELEMENT(gi,e[vi+j]);
1408 }
1409
1410 return g;
1411 }
1412
1413 /*****************************************************************************
1414 * *
1415 * copy_sg(sg1,sg2) makes a copy of sg1 into sg2. *
1416 * If sg2 is not NULL, it is assumed that the vlen,dlen,elen fields are *
1417 * correct and v,d,e are dynamically allocated (or NULL); they are *
1418 * reallocated if necessary. If sg2==NULL, a new structure is allocated. *
1419 * A pointer to the copy is returned. *
1420 * The new graph e component is the same, no compression is done. *
1421 * *
1422 *****************************************************************************/
1423
1424 sparsegraph*
copy_sg(sparsegraph * sg1,sparsegraph * sg2)1425 copy_sg(sparsegraph *sg1, sparsegraph *sg2)
1426 {
1427 int *d1,*e1,*d2,*e2;
1428 int i,n;
1429 size_t *v1,*v2,k;
1430 sg_weight *wt1,*wt2;
1431
1432 if (!sg2)
1433 {
1434 if ((sg2 = (sparsegraph*)ALLOCS(1,sizeof(sparsegraph))) == NULL)
1435 {
1436 fprintf(ERRFILE,"copy_sg: malloc failed\n");
1437 exit(1);
1438 }
1439 SG_INIT(*sg2);
1440 }
1441
1442 SWG_VDE(sg1,v1,d1,e1,wt1);
1443
1444 n = sg1->nv;
1445
1446 k = 0;
1447 for (i = 0; i < n; ++i)
1448 if (v1[i]+d1[i]>k) k = v1[i] + d1[i];
1449
1450 if (wt1)
1451 SWG_ALLOC(*sg2,n,k,"copy_sg malloc");
1452 else
1453 {
1454 SG_ALLOC(*sg2,n,k,"copy_sg malloc");
1455 DYNFREE(sg2->w,sg2->wlen);
1456 }
1457 SWG_VDE(sg2,v2,d2,e2,wt2);
1458
1459 sg2->nv = n;
1460 sg2->nde = sg1->nde;
1461 memcpy(v2,v1,n*sizeof(size_t));
1462 memcpy(d2,d1,n*sizeof(int));
1463 memcpy(e2,e1,k*sizeof(int));
1464 if (wt1) memcpy(wt2,wt1,k*sizeof(sg_weight));
1465
1466 return sg2;
1467 }
1468
1469 /*****************************************************************************
1470 * *
1471 * nauty_to_sg(g,sg,m,n) creates a sparse graph from a nauty format graph *
1472 * If sg is not NULL, it is assumed that the vlen,dlen,elen fields are *
1473 * correct and v,d,e are dynamically allocated (or NULL); they are *
1474 * reallocated if necessary. If sg==NULL, a new structure is allocated. *
1475 * A pointer to the sparse graph is returned. *
1476 * *
1477 *****************************************************************************/
1478
1479 sparsegraph*
nauty_to_sg(graph * g,sparsegraph * sg,int m,int n)1480 nauty_to_sg(graph *g, sparsegraph *sg, int m, int n)
1481 {
1482 int *d,*e;
1483 int i,k;
1484 set *gi;
1485 size_t j,*v,nde;
1486
1487 if (!sg)
1488 {
1489 if ((sg = (sparsegraph*)ALLOCS(1,sizeof(sparsegraph))) == NULL)
1490 {
1491 fprintf(ERRFILE,"nauty_to_sg: malloc failed\n");
1492 exit(1);
1493 }
1494 SG_INIT(*sg);
1495 }
1496
1497 nde = 0;
1498 for (gi = g + (size_t)m*(size_t)n; --gi >= g; )
1499 if (*gi != 0) nde += POPCOUNT(*gi);
1500
1501 sg->nv = n;
1502 sg->nde = nde;
1503
1504 SG_ALLOC(*sg,n,nde,"nauty_to_sg");
1505
1506 SG_VDE(sg,v,d,e);
1507
1508 j = 0;
1509 for (i = 0, gi = g; i < n; ++i, gi += m)
1510 {
1511 v[i] = j;
1512 for (k = -1; (k = nextelement(gi,m,k)) >= 0; )
1513 e[j++] = k;
1514 d[i] = j - v[i];
1515 }
1516
1517 return sg;
1518 }
1519
1520 /*****************************************************************************
1521 * *
1522 * distances_sg() assigns to each vertex v a value depending on the number *
1523 * of vertices at each distance from v, and what cells they lie in. *
1524 * If we find any cell which is split in this manner, we don't try any *
1525 * further cells. *
1526 * *
1527 *****************************************************************************/
1528
1529 void
distances_sg(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1530 distances_sg(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1531 int *invar, int invararg, boolean digraph, int m, int n)
1532 {
1533 int *d,*e;
1534 int i,k,dlim,wt;
1535 int di;
1536 int cell1,cell2,iv,liv,kcode;
1537 int head,tail;
1538 long longcode;
1539 size_t *v,vi,j;
1540 boolean success;
1541
1542 SG_VDE(g,v,d,e);
1543
1544 #if !MAXN
1545 DYNALLOC1(int,work1,work1_sz,n,"distances_sg");
1546 DYNALLOC1(int,work4,work4_sz,n,"distances_sg");
1547 DYNALLOC1(int,work3,work3_sz,n,"distances_sg");
1548 #endif
1549 #define CELLCODE work1
1550 #define QUEUE work4
1551 #define DIST work3
1552
1553 for (i = n; --i >= 0;) invar[i] = 0;
1554
1555 wt = 1;
1556 for (i = 0; i < n; ++i)
1557 {
1558 CELLCODE[lab[i]] = FUZZ1(wt);
1559 if (ptn[i] <= level) ++wt;
1560 }
1561
1562 if (invararg > n || invararg == 0) dlim = n;
1563 else dlim = invararg+1;
1564
1565 success = FALSE;
1566 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
1567 {
1568 for (cell2 = cell1; ptn[cell2] > level; ++cell2) {}
1569 if (cell2 == cell1) continue;
1570
1571 for (iv = cell1; iv <= cell2; ++iv)
1572 {
1573 liv = lab[iv];
1574 QUEUE[0] = liv;
1575 DIST[liv] = 0;
1576 RESETMARKS1;
1577 MARK1(liv);
1578 longcode = 0;
1579 head = 0;
1580 tail = 1;
1581
1582 while (tail < n && head < tail)
1583 {
1584 i = QUEUE[head++];
1585 if (DIST[i] >= dlim) break;
1586 vi = v[i];
1587 di = d[i];
1588
1589 for (j = 0; j < di; ++j)
1590 {
1591 k = e[vi+j];
1592 if (ISNOTMARKED1(k))
1593 {
1594 MARK1(k);
1595 DIST[k] = DIST[i] + 1;
1596 kcode = DIST[k]+CELLCODE[k];
1597 ACCUM(longcode,FUZZ1(kcode));
1598 QUEUE[tail++] = k;
1599 }
1600 }
1601 }
1602 invar[liv] = CLEANUP(longcode);
1603 if (invar[liv] != invar[lab[cell1]]) success = TRUE;
1604 }
1605 if (success) break;
1606 }
1607 }
1608
1609 /*****************************************************************************
1610 * *
1611 * adjacencies_sg() assigns to each vertex v a code depending on which cells *
1612 * it is joined to and from, and how many times. It is intended to provide *
1613 * better partitioning that the normal refinement routine for digraphs. *
1614 * It will not help with undirected graphs in nauty at all. *
1615 * *
1616 *****************************************************************************/
1617
1618 void
adjacencies_sg(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1619 adjacencies_sg(graph *g, int *lab, int *ptn, int level, int numcells,
1620 int tvpos, int *invar, int invararg, boolean digraph,
1621 int m, int n)
1622 {
1623 int *d,*e;
1624 int vwt,wwt;
1625 int *ei,di,i;
1626 size_t *v,j;
1627
1628 SG_VDE(g,v,d,e);
1629
1630 #if !MAXN
1631 DYNALLOC1(int,work2,work2_sz,n,"adjacencies_sg");
1632 #endif
1633
1634 vwt = 1;
1635 for (i = 0; i < n; ++i)
1636 {
1637 work2[lab[i]] = vwt;
1638 if (ptn[i] <= level) ++vwt;
1639 invar[i] = 0;
1640 }
1641
1642 for (i = 0; i < n; ++i)
1643 {
1644 vwt = FUZZ1(work2[i]);
1645 wwt = 0;
1646 di = d[i];
1647 ei = e + v[i];
1648 for (j = 0; j < di; ++j)
1649 {
1650 ACCUM(wwt,FUZZ2(work2[ei[j]]));
1651 ACCUM(invar[ei[j]],vwt);
1652 }
1653 ACCUM(invar[i],wwt);
1654 }
1655 }
1656
1657 /*****************************************************************************
1658 * *
1659 * sparsenauty(g,lab,ptn,orbits,&options,&stats,h) *
1660 * is a slightly simplified interface to nauty(). It allocates enough *
1661 * workspace for 500 automorphisms and checks that the sparsegraph dispatch *
1662 * vector is in use. *
1663 * *
1664 *****************************************************************************/
1665
1666 void
sparsenauty(sparsegraph * g,int * lab,int * ptn,int * orbits,optionblk * options,statsblk * stats,sparsegraph * h)1667 sparsenauty(sparsegraph *g, int *lab, int *ptn, int *orbits,
1668 optionblk *options, statsblk *stats, sparsegraph *h)
1669 {
1670 int m,n;
1671
1672 if (options->dispatch != &dispatch_sparse)
1673 {
1674 fprintf(ERRFILE,"Error: sparsenauty() needs standard options block\n");
1675 exit(1);
1676 }
1677
1678 n = g->nv;
1679 m = SETWORDSNEEDED(n);
1680
1681 #if !MAXN
1682 /* Don't increase 2*500*m in the following without also increasing
1683 the static decalaration of snwork[] above. */
1684 DYNALLOC1(set,snwork,snwork_sz,2*500*m,"densenauty malloc");
1685 #endif
1686
1687 nauty((graph*)g,lab,ptn,NULL,orbits,options,stats,
1688 snwork,2*500*m,m,n,(graph*)h);
1689 }
1690
1691 /*****************************************************************************
1692 * *
1693 * nausparse_check() checks that this file is compiled compatibly with the *
1694 * given parameters. If not, call exit(1). *
1695 * *
1696 *****************************************************************************/
1697
1698 void
nausparse_check(int wordsize,int m,int n,int version)1699 nausparse_check(int wordsize, int m, int n, int version)
1700 {
1701 if (wordsize != WORDSIZE)
1702 {
1703 fprintf(ERRFILE,"Error: WORDSIZE mismatch in nausparse.c\n");
1704 exit(1);
1705 }
1706
1707 #if MAXN
1708 if (m > MAXM)
1709 {
1710 fprintf(ERRFILE,"Error: MAXM inadequate in nausparse.c\n");
1711 exit(1);
1712 }
1713
1714 if (n > MAXN)
1715 {
1716 fprintf(ERRFILE,"Error: MAXN inadequate in nausparse.c\n");
1717 exit(1);
1718 }
1719 #endif
1720
1721 if (version < NAUTYREQUIRED)
1722 {
1723 fprintf(ERRFILE,"Error: nausparse.c version mismatch\n");
1724 exit(1);
1725 }
1726 }
1727
1728 /*****************************************************************************
1729 * *
1730 * nausparse_freedyn() - free the dynamic memory in this module *
1731 * *
1732 *****************************************************************************/
1733
1734 void
nausparse_freedyn(void)1735 nausparse_freedyn(void)
1736 {
1737 #if !MAXN
1738 DYNFREE(vmark1,vmark1_sz);
1739 DYNFREE(vmark2,vmark2_sz);
1740 DYNFREE(work1,work1_sz);
1741 DYNFREE(work2,work2_sz);
1742 DYNFREE(work3,work3_sz);
1743 DYNFREE(work4,work4_sz);
1744 DYNFREE(snwork,snwork_sz);
1745 #endif
1746 }
1747