/***************************************************************************** * * * Sparse-graph-specific auxiliary source file for version 2.7 of nauty. * * * * Copyright (2004-2016) Brendan McKay. All rights reserved. * * Subject to waivers and disclaimers in nauty.h. * * * * CHANGE HISTORY * * 26-Oct-04 : initial creation * * 23-Nov-06 : dispatch uses targetcell_sg, not bestcell_sg * * 8-Dec-06 : add adjacencies_sg() * * 10-Nov-09 : remove types shortish and permutation * * 14-Nov-09 : added copy_sg() * * 11-May-10 : use sorttemplates.c for sorting procedures * * 19-May-10 : add two *_tr procedures for traces. * * 21-May-10 : fixes for nde,v fields becoming size_t * * 23-May-10 : add sparsenauty() * * 15-Jan-12 : add TLS_ATTR attributes * * 17-Dec-15 : extend sortlists_sg() to sort weights * * : add weights to copy_sg() and updatecan_sg() * * 11-Mar-16 : add cleanup_sg(). This can be used in the cleanup * * field of the dispatch vector to sort the lists of the * * canonical graph, but isn't there by default. * * 15-Oct-19 : fix static declaration of snwork[] * * * *****************************************************************************/ #define TMP /* #define ONE_WORD_SETS not sure about this! See notes.txt. */ #include "nausparse.h" /* macros for hash-codes: */ #define MASH(l,i) ((((l) ^ 065435) + (i)) & 077777) /* : expression whose long value depends only on long l and int/long i. Anything goes, preferably non-commutative. */ #define CLEANUP(l) ((int)((l) % 077777)) /* : expression whose value depends on long l and is less than 077777 when converted to int then short. Anything goes. */ #if MAXM==1 #define M 1 #else #define M m #endif #define ACCUM(x,y) x = (((x) + (y)) & 077777) static const int fuzz1[] = {037541,061532,005257,026416}; static const int fuzz2[] = {006532,070236,035523,062437}; #define FUZZ1(x) ((x) ^ fuzz1[(x)&3]) #define FUZZ2(x) ((x) ^ fuzz2[(x)&3]) /* aproto: header new_nauty_protos.h */ dispatchvec dispatch_sparse = {isautom_sg,testcanlab_sg,updatecan_sg,refine_sg,refine_sg,cheapautom_sg, targetcell_sg,nausparse_freedyn,nausparse_check,init_sg,NULL}; #if !MAXN DYNALLSTAT(short,vmark1,vmark1_sz); DYNALLSTAT(short,vmark2,vmark2_sz); DYNALLSTAT(int,work1,work1_sz); DYNALLSTAT(int,work2,work2_sz); DYNALLSTAT(int,work3,work3_sz); DYNALLSTAT(int,work4,work4_sz); DYNALLSTAT(set,snwork,snwork_sz); #else static TLS_ATTR short vmark1[MAXN]; static TLS_ATTR short vmark2[MAXN]; static TLS_ATTR int work1[MAXN]; static TLS_ATTR int work2[MAXN]; static TLS_ATTR int work3[MAXN]; static TLS_ATTR int work4[MAXN]; static TLS_ATTR set snwork[2*60*MAXM]; #endif static TLS_ATTR short vmark1_val = 32000; #define MARK1(i) vmark1[i] = vmark1_val #define UNMARK1(i) vmark1[i] = 0 #define ISMARKED1(i) (vmark1[i] == vmark1_val) #define ISNOTMARKED1(i) (vmark1[i] != vmark1_val) static TLS_ATTR short vmark2_val = 32000; #define MARK2(i) vmark2[i] = vmark2_val #define UNMARK2(i) vmark2[i] = 0 #define ISMARKED2(i) (vmark2[i] == vmark2_val) #define ISNOTMARKED2(i) (vmark2[i] != vmark2_val) #if !MAXN #define RESETMARKS1 {if (vmark1_val++ >= 32000) \ {size_t ij; for (ij=0;ij= 32000) \ {size_t ij; for (ij=0;ij= 32000) \ {size_t ij; for (ij=0;ij= 32000) \ {size_t ij; for (ij=0;ijnv; if (g2->nv != n || g2->nde != g1->nde) return FALSE; SG_VDE(g1,v1,d1,e1); SG_VDE(g2,v2,d2,e2); PREPAREMARKS1(n); for (i = 0; i < n; ++i) { di = d1[i]; if (d2[i] != di) return FALSE; vi = v1[i]; RESETMARKS1; for (j = 0; j < di; ++j) MARK1(e1[vi+j]); vi = v2[i]; for (j = 0; j < di; ++j) if (ISNOTMARKED1(e2[vi+j])) return FALSE; } return TRUE; } /***************************************************************************** * * * testcanlab_sg(g,canong,lab,samerows,m,n) compares g^lab to canong, * * using an ordering which is immaterial since it's only used here. The * * value returned is -1,0,1 if g^lab <,=,> canong. *samerows is set to * * the number of rows (0..n) of canong which are the same as those of g^lab. * * * *****************************************************************************/ int testcanlab_sg(graph *g, graph *canong, int *lab, int *samerows, int m, int n) { int *d,*e; int *cd,*ce; int i,k,di,dli; size_t j,vi,vli,*v,*cv; int mina; SG_VDE(g,v,d,e); SG_VDE(canong,cv,cd,ce); #if !MAXN DYNALLOC1(int,work1,work1_sz,n,"testcanlab_sg"); #endif #define INVLAB work1 PREPAREMARKS1(n); for (i = 0; i < n; ++i) INVLAB[lab[i]] = i; for (i = 0; i < n; ++i) { /* compare g[lab[i]]^INVLAB to canong[i] */ vi = cv[i]; di = cd[i]; vli = v[lab[i]]; dli = d[lab[i]]; if (di != dli) { *samerows = i; if (di < dli) return -1; return 1; } RESETMARKS1; mina = n; for (j = 0; j < di; ++j) MARK1(ce[vi+j]); for (j = 0; j < di; ++j) { k = INVLAB[e[vli+j]]; if (ISMARKED1(k)) UNMARK1(k); else if (k < mina) mina = k; } if (mina != n) { *samerows = i; for (j = 0; j < di; ++j) { k = ce[vi+j]; if (ISMARKED1(k) && k < mina) return -1; } return 1; } } *samerows = n; return 0; } /***************************************************************************** * * * updatecan_sg(g,canong,lab,samerows,m,n) sets canong = g^lab, assuming * * the first samerows vertices of canong are ok already. Also assumes * * contiguity and ample space in canong. * * * *****************************************************************************/ void updatecan_sg(graph *g, graph *canong, int *lab, int samerows, int m, int n) { int *d,*e; int *cd,*ce; int i,dli; size_t *v,*cv,vli,j,k; sg_weight *wt,*cwt; SWG_VDE(g,v,d,e,wt); SWG_VDE(canong,cv,cd,ce,cwt); #if !MAXN DYNALLOC1(int,work1,work1_sz,n,"testcanlab_sg"); #endif #define INVLAB work1 ((sparsegraph*)canong)->nv = n; ((sparsegraph*)canong)->nde = ((sparsegraph*)g)->nde; for (i = 0; i < n; ++i) INVLAB[lab[i]] = i; if (samerows == 0) k = 0; else k = cv[samerows-1]+cd[samerows-1]; for (i = samerows; i < n; ++i) { cv[i] = k; cd[i] = dli = d[lab[i]]; vli = v[lab[i]]; if (wt) { for (j = 0; j < dli; ++j) { ce[k] = INVLAB[e[vli+j]]; cwt[k] = wt[vli+j]; ++k; } } else for (j = 0; j < dli; ++j) ce[k++] = INVLAB[e[vli+j]]; } } /***************************************************************************** * * * comparelab_tr(g,lab1,invlab1,lab2,invlab2,cls,col) compares * * g^lab1 to g^lab2 and returns -1,0,1 according to the comparison. * * invlab1[] and invlab2[] are assumed to hold inverses of lab1,lab2. * * * *****************************************************************************/ int comparelab_tr(sparsegraph *g, int *lab1, int *invlab1, int *lab2, int *invlab2, int *cls, int *col) { int d1,*e1,d2,*e2; int i,j,k,n,c,end; int mina; n = g->nv; PREPAREMARKS1(n); for (c=0; ce + g->v[lab1[i]]; d1 = g->d[lab1[i]]; e2 = g->e + g->v[lab2[i]]; d2 = g->d[lab2[i]]; if (d1 < d2) return -1; else if (d1 > d2) return 1; RESETMARKS1; mina = n; for (j = 0; j < d1; ++j) MARK1(col[invlab1[e1[j]]]); for (j = 0; j < d1; ++j) { k = col[invlab2[e2[j]]]; if (ISMARKED1(k)) UNMARK1(k); else if (k < mina) mina = k; } if (mina != n) { for (j = 0; j < d1; ++j) { k = col[invlab1[e1[j]]]; if (ISMARKED1(k) && k < mina) return -1; } return 1; } } } } return 0; } /***************************************************************************** * * * testcanlab_tr(g,canong,lab,invlab,samerows) compares g^lab to canong, * * using an ordering which is immaterial since it's only used here. The * * value returned is -1,0,1 if g^lab <,=,> canong. *samerows is set to * * the number of rows (0..n) of canong which are the same as those of g^lab. * * invlab[] is assumed to hold the inverse of lab[] * * * *****************************************************************************/ int testcanlab_tr(sparsegraph *g, sparsegraph *canong, int *lab, int *invlab, int *samerows) { int *d,*e; int *cd,*ce; int i,di,dli; int k,n; size_t *v,*cv,vi,vli,j; int mina; SG_VDE(g,v,d,e); SG_VDE(canong,cv,cd,ce); n = g->nv; PREPAREMARKS1(n); for (i = 0; i < n; ++i) { /* compare g[lab[i]]^invlab to canong[i] */ vi = cv[i]; di = cd[i]; vli = v[lab[i]]; dli = d[lab[i]]; if (di != dli) { *samerows = i; if (di < dli) return -1; return 1; } RESETMARKS1; mina = n; for (j = 0; j < di; ++j) MARK1(ce[vi+j]); for (j = 0; j < di; ++j) { k = invlab[e[vli+j]]; if (ISMARKED1(k)) UNMARK1(k); else if (k < mina) mina = k; } if (mina != n) { *samerows = i; for (j = 0; j < di; ++j) { k = ce[vi+j]; if (ISMARKED1(k) && k < mina) return -1; } return 1; } } *samerows = n; return 0; } /***************************************************************************** * * * updatecan_tr(g,canong,lab,invlab,samerows) sets canong = g^lab, * * assuming the first samerows vertices of canong are ok already. * * Also assumes contiguity and ample space in canong. * * Assumes invlab[] holds the inverse of lab[] * * * *****************************************************************************/ void updatecan_tr(sparsegraph *g, sparsegraph *canong, int *lab, int *invlab, int samerows) { int *d,*e; int *cd,*ce; int i,dli,n; size_t *v,*cv,vli,j,k; SG_VDE(g,v,d,e); SG_VDE(canong,cv,cd,ce); n = g->nv; PREPAREMARKS1(n); canong->nv = n; canong->nde = g->nde; if (samerows == 0) k = 0; else k = cv[samerows-1]+cd[samerows-1]; for (i = samerows; i < n; ++i) { cv[i] = k; cd[i] = dli = d[lab[i]]; vli = v[lab[i]]; for (j = 0; j < dli; ++j) ce[k++] = invlab[e[vli+j]]; } } #define SORT_OF_SORT 3 #define SORT_NAME sortindirect #define SORT_TYPE1 int #define SORT_TYPE2 int #include "sorttemplates.c" #define SORT_OF_SORT 1 #define SORT_NAME sortints #define SORT_TYPE1 int #include "sorttemplates.c" #define SORT_OF_SORT 2 #define SORT_NAME sortweights #define SORT_TYPE1 int #define SORT_TYPE2 sg_weight #include "sorttemplates.c" /***************************************************************************** * * * init_sg(graph *gin, graph **gout, graph *hin, graph **hout, * * int *lab, int *ptn, set *active, optionblk *options, * * int *status, int m, int n) * * Initialise routine for dispatch vector. This one just makes sure * * that *hin has enough space and sets fields for n=0. * * * *****************************************************************************/ 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) { sparsegraph *sg,*sh; if (options->getcanon) { sg = (sparsegraph*)gin; sh = (sparsegraph*)hin; SG_ALLOC(*sh,sg->nv,sg->nde,"init_sg"); sh->nv = sg->nv; sh->nde = sg->nde; } *status = 0; } /***************************************************************************** * * * cleanup_sg(graph *gin, graph **gout, graph *hin, graph **hout, * * int *lab, int *ptn, optionblk *options, * * statsblk *stats, int m, int n) * * Cleanup routine for dispatch vector. This one sorts the adjacency * * lists for the canonical labelling. * * * *****************************************************************************/ void cleanup_sg(graph *gin, graph **gout, graph *hin, graph **hout, int *lab, int *ptn, optionblk *options, statsblk *stats, int m, int n) { sparsegraph *sh; if (options->getcanon && (stats->errstatus == 0 || stats->errstatus == NAUABORTED)) { sh = (sparsegraph*)hin; sortlists_sg(sh); } } /***************************************************************************** * * * distvals(sparsegraph *sg, int v0, int *dist, int n) sets dist[i] * * to the distance from v0 to i, for each i, or to n if there is no such * * distance. work4[] is used as a queue. * * * *****************************************************************************/ void distvals(sparsegraph *g, int v0, int *dist, int n) { int *d,*e; int i,head,tail; int di,k; size_t *v,vi,j; SG_VDE(g,v,d,e); #if !MAXN DYNALLOC1(int,work4,work4_sz,n,"distvals"); #endif #define QUEUE work4 for (i = 0; i < n; ++i) dist[i] = n; QUEUE[0] = v0; dist[v0] = 0; head = 0; tail = 1; while (tail < n && head < tail) { i = QUEUE[head++]; vi = v[i]; di = d[i]; for (j = 0; j < di; ++j) { k = e[vi+j]; if (dist[k] == n) { dist[k] = dist[i] + 1; QUEUE[tail++] = k; } } } } /***************************************************************************** * * * refine_sg(g,lab,ptn,level,numcells,count,active,code,m,n) performs a * * refinement operation on the partition at the specified level of the * * partition nest (lab,ptn). *numcells is assumed to contain the number of * * cells on input, and is updated. The initial set of active cells (alpha * * in the paper) is specified in the set active. Precisely, x is in active * * iff the cell starting at index x in lab is active. * * The resulting partition is equitable if active is correct (see the paper * * and the Guide). * * *code is set to a value which depends on the fine detail of the * * algorithm, but which is independent of the labelling of the graph. * * count is used for work space. * * * *****************************************************************************/ void refine_sg(graph *g, int *lab, int *ptn, int level, int *numcells, int *count, set *active, int *code, int m, int n) { int i,j,k,l,v1,v2,v3,isplit; int w1,w2,w3; long longcode; int *d,*e; int size,bigsize,bigpos; int nactive,hitcells; int lj,di,splitv; boolean trivsplit; size_t *v,vi,ii; SG_VDE(g,v,d,e); #if !MAXN DYNALLOC1(int,work1,work1_sz,n,"refine_sg"); DYNALLOC1(int,work2,work2_sz,n,"refine_sg"); DYNALLOC1(int,work3,work3_sz,n,"refine_sg"); DYNALLOC1(int,work4,work4_sz,n,"refine_sg"); #endif #define CELLSTART work1 #define ACTIVE work2 #define HITS work3 #define HITCELL work4 PREPAREMARKS1(n); PREPAREMARKS2(n); longcode = *numcells; /* Set ACTIVE[0..nactive-1] = queue of active cell starts */ nactive = 0; for (i = -1; (i = nextelement(active,m,i)) >= 0;) ACTIVE[nactive++] = i; if (nactive == 0) { *code = CLEANUP(longcode); return; } /* Set CELLSTART[i] = starting point in lab[] of nontrivial cell containing i, or n if i is a singleton */ for (i = 0; i < n; ) { /* Just here, i is a cell starting position */ if (ptn[i] <= level) { CELLSTART[lab[i]] = n; ++i; } else { j = i; do { CELLSTART[lab[i]] = j; } while (ptn[i++] > level); } } if (level <= 2 && nactive == 1 && ptn[ACTIVE[0]] <= level && *numcells <= n/8) { isplit = ACTIVE[--nactive]; DELELEMENT(active,isplit); distvals((sparsegraph*)g,lab[isplit],HITS,n); for (v1 = 0; v1 < n; ) { if (ptn[v1] <= level) { ++v1; continue; } longcode = MASH(longcode,v1); w1 = HITS[lab[v1]]; v2 = v1+1; while (ptn[v2-1] > level && HITS[lab[v2]] == w1) ++v2; if (ptn[v2-1] <= level) { v1 = v2; continue; } w2 = NAUTY_INFINITY; v3 = j = v2; do { lj = lab[j]; w3 = HITS[lj]; if (w3 == w1) { lab[j] = lab[v3]; lab[v3] = lab[v2]; lab[v2] = lj; ++v2; ++v3; } else if (w3 == w2) { lab[j] = lab[v3]; lab[v3] = lj; ++v3; } else if (w3 < w1) { lab[j] = lab[v2]; lab[v2] = lab[v1]; lab[v1] = lj; v3 = v2 + 1; v2 = v1 + 1; w2 = w1; w1 = w3; } else if (w3 < w2) { lab[j] = lab[v2]; lab[v2] = lj; v3 = v2 + 1; w2 = w3; } } while (ptn[j++] > level); longcode = MASH(longcode,w2); longcode = MASH(longcode,v2); if (j != v2) /* At least two fragments * v1..v2-1 = w1; v2..v3-1 = w2 */ { if (v2 == v1+1) CELLSTART[lab[v1]] = n; if (v3 == v2+1) CELLSTART[lab[v2]] = n; else for (k = v2; k < v3; ++k) CELLSTART[lab[k]] = v2; ++*numcells; ptn[v2-1] = level; if (j == v3) { /* Two fragments only */ if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1)) { ADDELEMENT(active,v1); ACTIVE[nactive++] = v1; } else { ADDELEMENT(active,v2); ACTIVE[nactive++] = v2; } } else { /* Extra fragments: v3..j-1 > w2 */ sortindirect(lab+v3,HITS,j-v3); ACTIVE[nactive++] = v2; ADDELEMENT(active,v2); if (v2-v1 >= v3-v2) { bigpos = -1; bigsize = v2-v1; } else { bigpos = nactive-1; bigsize = v3-v2; } for (k = v3-1; k < j-1;) { ptn[k] = level; longcode = MASH(longcode,k); ++*numcells; l = k+1; ADDELEMENT(active,l); ACTIVE[nactive++] = l; w3 = HITS[lab[l]]; for (k = l; k < j-1 && HITS[lab[k+1]] == w3; ++k) CELLSTART[lab[k+1]] = l; size = k-l+1; if (size == 1) CELLSTART[lab[l]] = n; else { CELLSTART[lab[l]] = l; if (size > bigsize) { bigsize = size; bigpos = nactive-1; } } } if (bigpos >= 0 && !ISELEMENT(active,v1)) { longcode = MASH(longcode,bigpos); DELELEMENT(active,ACTIVE[bigpos]); ADDELEMENT(active,v1); ACTIVE[bigpos] = v1; } } } v1 = j; } } /* Iterate until complete */ while (nactive > 0 && *numcells < n) { for (i = 0; i < nactive && i < 10; ++i) if (ptn[ACTIVE[i]] <= level) break; if (i < nactive && i < 10) { trivsplit = TRUE; isplit = ACTIVE[i]; ACTIVE[i] = ACTIVE[--nactive]; } else { isplit = ACTIVE[--nactive]; trivsplit = ptn[isplit] <= level; } DELELEMENT(active,isplit); longcode = MASH(longcode,isplit); if (trivsplit) { RESETMARKS1; RESETMARKS2; hitcells = 0; splitv = lab[isplit]; vi = v[splitv]; di = d[splitv]; for (ii = 0; ii < di; ++ii) { j = e[vi+ii]; MARK2(j); k = CELLSTART[j]; if (k != n && ISNOTMARKED1(k)) { MARK1(k); HITCELL[hitcells++] = k; } } if (hitcells > 1) sortints(HITCELL,hitcells); longcode = MASH(longcode,hitcells); /* divide cells according to which vertices are hit */ for (i = 0; i < hitcells; ++i) { j = v1 = v2 = HITCELL[i]; longcode = MASH(longcode,v2); k = 0; do { lj = lab[j]; if (ISMARKED2(lj)) HITS[k++] = lj; else lab[v2++] = lj; } while (ptn[j++] > level); longcode = MASH(longcode,k); v3 = v2; while (--k >= 0) { j = HITS[k]; CELLSTART[j] = v2; lab[v3++] = j; } if (v2 != v3 && v2 != v1) { ++*numcells; if (v2 == v1+1) CELLSTART[lab[v1]] = n; if (v3 == v2+1) CELLSTART[lab[v2]] = n; ptn[v2-1] = level; longcode = MASH(longcode,v2); if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1)) { ADDELEMENT(active,v1); ACTIVE[nactive++] = v1; } else { ADDELEMENT(active,v2); ACTIVE[nactive++] = v2; } } } } else /* non-trivial splitting */ { /* isplit is the start of the splitting cell. Set HITS[i] = hits of i for i in non-trivial cells, HITCELL[0..hitcells-1] = starts of hit non-trivial cells */ RESETMARKS1; hitcells = 0; do { vi = v[lab[isplit]]; di = d[lab[isplit]]; for (ii = 0; ii < di; ++ii) { j = e[vi+ii]; k = CELLSTART[j]; if (k != n) { if (ISNOTMARKED1(k)) { MARK1(k); HITCELL[hitcells++] = k; do { HITS[lab[k]] = 0; } while (ptn[k++] > level); } ++HITS[j]; } } } while (ptn[isplit++] > level); if (hitcells > 1) sortints(HITCELL,hitcells); /* divide cells according to hit counts */ longcode = MASH(longcode,hitcells); for (i = 0; i < hitcells; ++i) { v1 = HITCELL[i]; w1 = HITS[lab[v1]]; longcode = MASH(longcode,v1); v2 = v1+1; while (ptn[v2-1] > level && HITS[lab[v2]] == w1) ++v2; if (ptn[v2-1] <= level) continue; w2 = NAUTY_INFINITY; v3 = j = v2; do { lj = lab[j]; w3 = HITS[lj]; if (w3 == w1) { lab[j] = lab[v3]; lab[v3] = lab[v2]; lab[v2] = lj; ++v2; ++v3; } else if (w3 == w2) { lab[j] = lab[v3]; lab[v3] = lj; ++v3; } else if (w3 < w1) { lab[j] = lab[v2]; lab[v2] = lab[v1]; lab[v1] = lj; v3 = v2 + 1; v2 = v1 + 1; w2 = w1; w1 = w3; } else if (w3 < w2) { lab[j] = lab[v2]; lab[v2] = lj; v3 = v2 + 1; w2 = w3; } } while (ptn[j++] > level); longcode = MASH(longcode,w1); longcode = MASH(longcode,v2); if (j != v2) /* At least two fragments * v1..v2-1 = w1; v2..v3-1 = w2 */ { if (v2 == v1+1) CELLSTART[lab[v1]] = n; if (v3 == v2+1) CELLSTART[lab[v2]] = n; else for (k = v2; k < v3; ++k) CELLSTART[lab[k]] = v2; ++*numcells; ptn[v2-1] = level; if (j == v3) { /* Two fragments only */ if (v2-v1 <= v3-v2 && !ISELEMENT(active,v1)) { ADDELEMENT(active,v1); ACTIVE[nactive++] = v1; } else { ADDELEMENT(active,v2); ACTIVE[nactive++] = v2; } } else { /* Extra fragments: v3..j-1 > w2 */ longcode = MASH(longcode,v3); sortindirect(lab+v3,HITS,j-v3); ACTIVE[nactive++] = v2; ADDELEMENT(active,v2); if (v2-v1 >= v3-v2) { bigpos = -1; bigsize = v2-v1; } else { bigpos = nactive-1; bigsize = v3-v2; longcode = MASH(longcode,bigsize); } for (k = v3-1; k < j-1;) { ptn[k] = level; ++*numcells; l = k+1; ADDELEMENT(active,l); ACTIVE[nactive++] = l; w3 = HITS[lab[l]]; longcode = MASH(longcode,w3); for (k = l; k < j-1 && HITS[lab[k+1]] == w3; ++k) CELLSTART[lab[k+1]] = l; size = k-l+1; if (size == 1) CELLSTART[lab[l]] = n; else { CELLSTART[lab[l]] = l; if (size > bigsize) { bigsize = size; bigpos = nactive-1; } } } if (bigpos >= 0 && !ISELEMENT(active,v1)) { DELELEMENT(active,ACTIVE[bigpos]); ADDELEMENT(active,v1); ACTIVE[bigpos] = v1; } } } } } } longcode = MASH(longcode,*numcells); *code = CLEANUP(longcode); } /***************************************************************************** * * * cheapautom_sg(ptn,level,digraph,n) returns TRUE if the partition at the * * specified level in the partition nest (lab,ptn) {lab is not needed here} * * satisfies a simple sufficient condition for its cells to be the orbits of * * some subgroup of the automorphism group. Otherwise it returns FALSE. * * It always returns FALSE if digraph!=FALSE. * * * * nauty assumes that this function will always return TRUE for any * * partition finer than one for which it returns TRUE. * * * *****************************************************************************/ boolean cheapautom_sg(int *ptn, int level, boolean digraph, int n) { int i,k,nnt; if (digraph) return FALSE; k = n; nnt = 0; for (i = 0; i < n; ++i) { --k; if (ptn[i] > level) { ++nnt; while (ptn[++i] > level) {} } } return (k <= nnt + 1 || k <= 4); } /***************************************************************************** * * * bestcell_sg(g,lab,ptn,level,tc_level,m,n) returns the index in lab of * * the start of the "best non-singleton cell" for fixing. If there is no * * non-singleton cell it returns n. * * This implementation finds the first cell which is non-trivially joined * * to the greatest number of other cells, assuming equitability. * * This is not good for digraphs! * * * *****************************************************************************/ static int bestcell_sg(graph *g, int *lab, int *ptn, int level, int tc_level, int m, int n) { int nnt; int *d,*e; int i,k,di; int *work1b; int maxcnt; size_t *v,vi,j; SG_VDE(g,v,d,e); #if !MAXN DYNALLOC1(int,work1,work1_sz,n,"bestcell_sg"); DYNALLOC1(int,work2,work2_sz,n,"bestcell_sg"); DYNALLOC1(int,work3,work3_sz,n,"bestcell_sg"); DYNALLOC1(int,work4,work4_sz,n,"bestcell_sg"); #endif work1b = work1 + (n/2); #define START work1 #define SIZE work1b #define NNTCELL work2 #define HITS work3 #define COUNT work4 /* find non-singleton cells: put starts in START[0..nnt-1], sizes in SIZE[0..nnt-1]. Also NNTCELL[i] = n if {i} is a singelton, else index of nontriv cell containing i. */ i = nnt = 0; while (i < n) { if (ptn[i] > level) { START[nnt] = i; j = i; do NNTCELL[lab[j]] = nnt; while (ptn[j++] > level); SIZE[nnt] = j-i; ++nnt; i = j; } else { NNTCELL[lab[i]] = n; ++i; } } if (nnt == 0) return n; /* set COUNT[i] to # non-trivial neighbours of n.s. cell i */ for (i = 0; i < nnt; ++i) HITS[i] = COUNT[i] = 0; for (i = 0; i < nnt; ++i) { vi = v[lab[START[i]]]; di = d[lab[START[i]]]; for (j = 0; j < di; ++j) { k = NNTCELL[e[vi+j]]; if (k != n) ++HITS[k]; } for (j = 0; j < di; ++j) { k = NNTCELL[e[vi+j]]; if (k != n) { if (HITS[k] > 0 && HITS[k] < SIZE[k]) ++COUNT[i]; HITS[k] = 0; } } } /* find first greatest bucket value */ j = 0; maxcnt = COUNT[0]; for (i = 1; i < nnt; ++i) if (COUNT[i] > maxcnt) { j = i; maxcnt = COUNT[i]; } return (int)START[j]; } /***************************************************************************** * * * targetcell_sg(g,lab,ptn,level,tc_level,digraph,hint,m,n) returns the * * index in lab of the next cell to split. * * hint is a suggestion for the answer, which is obeyed if it is valid. * * Otherwise we use bestcell() up to tc_level and the first non-trivial * * cell after that. * * * *****************************************************************************/ int targetcell_sg(graph *g, int *lab, int *ptn, int level, int tc_level, boolean digraph, int hint, int m, int n) { int i; if (hint >= 0 && ptn[hint] > level && (hint == 0 || ptn[hint-1] <= level)) return hint; else if (level <= tc_level) return bestcell_sg(g,lab,ptn,level,tc_level,m,n); else { for (i = 0; i < n && ptn[i] <= level; ++i) {} return (i == n ? 0 : i); } } /***************************************************************************** * * * sortlists_sg(g) sorts the adjacency lists into numerical order * * * *****************************************************************************/ void sortlists_sg(sparsegraph *g) { int *d,*e; int n,i; size_t *v; sg_weight *wt; SWG_VDE(g,v,d,e,wt); n = g->nv; if (wt) { for (i = 0; i < n; ++i) if (d[i] > 1) sortweights(e+v[i],wt+v[i],d[i]); } else { for (i = 0; i < n; ++i) if (d[i] > 1) sortints(e+v[i],d[i]); } } /***************************************************************************** * * * put_sg(f,sg,digraph,linelength) writes the sparse graph to file f using * * at most linelength characters per line. If digraph then all directed * * edges are written; else one v-w for w>=v is written. * * * *****************************************************************************/ void put_sg(FILE *f, sparsegraph *sg, boolean digraph, int linelength) { int *d,*e; int n,di; int i,curlen,slen; size_t *v,vi,j; char s[12]; SG_VDE(sg,v,d,e); n = sg->nv; for (i = 0; i < n; ++i) { vi = v[i]; di = d[i]; if (di == 0) continue; slen = itos(i+labelorg,s); putstring(f,s); putstring(f," :"); curlen = slen + 2; for (j = 0; j < di; ++j) { if (!digraph && e[vi+j] < i) continue; slen = itos(e[vi+j]+labelorg,s); if (linelength && curlen + slen + 1 >= linelength) { putstring(f,"\n "); curlen = 2; } PUTC(' ',f); putstring(f,s); curlen += slen + 1; } PUTC('\n',f); } } /***************************************************************************** * * * sg_to_nauty(sg,g,reqm,&m) creates a nauty-format graph from a sparse * * graph. reqm is the required m value (computed from n if reqm=0), and * * m is the actual value used. g is dynamically generated if NULL is given. * * A pointer to g is returned. * * * *****************************************************************************/ graph* sg_to_nauty(sparsegraph *sg, graph *g, int reqm, int *pm) { int *d,*e; int m,n,i,di; size_t *v,vi,j; set *gi; SG_VDE(sg,v,d,e); n = sg->nv; if (reqm != 0 && reqm*WORDSIZE < n) { fprintf(ERRFILE,"sg_to_nauty: reqm is impossible\n"); exit(1); } if (reqm != 0) m = reqm; else m = (n+WORDSIZE-1)/WORDSIZE; *pm = m; if (g == NULL) { if ((g = (graph*)ALLOCS(n,m*sizeof(graph))) == NULL) { fprintf(ERRFILE,"sg_to_nauty: malloc failed\n"); exit(1); } } for (i = 0, gi = g; i < n; ++i, gi += m) { vi = v[i]; di = d[i]; EMPTYSET(gi,m); for (j = 0; j < di; ++j) ADDELEMENT(gi,e[vi+j]); } return g; } /***************************************************************************** * * * copy_sg(sg1,sg2) makes a copy of sg1 into sg2. * * If sg2 is not NULL, it is assumed that the vlen,dlen,elen fields are * * correct and v,d,e are dynamically allocated (or NULL); they are * * reallocated if necessary. If sg2==NULL, a new structure is allocated. * * A pointer to the copy is returned. * * The new graph e component is the same, no compression is done. * * * *****************************************************************************/ sparsegraph* copy_sg(sparsegraph *sg1, sparsegraph *sg2) { int *d1,*e1,*d2,*e2; int i,n; size_t *v1,*v2,k; sg_weight *wt1,*wt2; if (!sg2) { if ((sg2 = (sparsegraph*)ALLOCS(1,sizeof(sparsegraph))) == NULL) { fprintf(ERRFILE,"copy_sg: malloc failed\n"); exit(1); } SG_INIT(*sg2); } SWG_VDE(sg1,v1,d1,e1,wt1); n = sg1->nv; k = 0; for (i = 0; i < n; ++i) if (v1[i]+d1[i]>k) k = v1[i] + d1[i]; if (wt1) SWG_ALLOC(*sg2,n,k,"copy_sg malloc"); else { SG_ALLOC(*sg2,n,k,"copy_sg malloc"); DYNFREE(sg2->w,sg2->wlen); } SWG_VDE(sg2,v2,d2,e2,wt2); sg2->nv = n; sg2->nde = sg1->nde; memcpy(v2,v1,n*sizeof(size_t)); memcpy(d2,d1,n*sizeof(int)); memcpy(e2,e1,k*sizeof(int)); if (wt1) memcpy(wt2,wt1,k*sizeof(sg_weight)); return sg2; } /***************************************************************************** * * * nauty_to_sg(g,sg,m,n) creates a sparse graph from a nauty format graph * * If sg is not NULL, it is assumed that the vlen,dlen,elen fields are * * correct and v,d,e are dynamically allocated (or NULL); they are * * reallocated if necessary. If sg==NULL, a new structure is allocated. * * A pointer to the sparse graph is returned. * * * *****************************************************************************/ sparsegraph* nauty_to_sg(graph *g, sparsegraph *sg, int m, int n) { int *d,*e; int i,k; set *gi; size_t j,*v,nde; if (!sg) { if ((sg = (sparsegraph*)ALLOCS(1,sizeof(sparsegraph))) == NULL) { fprintf(ERRFILE,"nauty_to_sg: malloc failed\n"); exit(1); } SG_INIT(*sg); } nde = 0; for (gi = g + (size_t)m*(size_t)n; --gi >= g; ) if (*gi != 0) nde += POPCOUNT(*gi); sg->nv = n; sg->nde = nde; SG_ALLOC(*sg,n,nde,"nauty_to_sg"); SG_VDE(sg,v,d,e); j = 0; for (i = 0, gi = g; i < n; ++i, gi += m) { v[i] = j; for (k = -1; (k = nextelement(gi,m,k)) >= 0; ) e[j++] = k; d[i] = j - v[i]; } return sg; } /***************************************************************************** * * * distances_sg() assigns to each vertex v a value depending on the number * * of vertices at each distance from v, and what cells they lie in. * * If we find any cell which is split in this manner, we don't try any * * further cells. * * * *****************************************************************************/ 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) { int *d,*e; int i,k,dlim,wt; int di; int cell1,cell2,iv,liv,kcode; int head,tail; long longcode; size_t *v,vi,j; boolean success; SG_VDE(g,v,d,e); #if !MAXN DYNALLOC1(int,work1,work1_sz,n,"distances_sg"); DYNALLOC1(int,work4,work4_sz,n,"distances_sg"); DYNALLOC1(int,work3,work3_sz,n,"distances_sg"); #endif #define CELLCODE work1 #define QUEUE work4 #define DIST work3 for (i = n; --i >= 0;) invar[i] = 0; wt = 1; for (i = 0; i < n; ++i) { CELLCODE[lab[i]] = FUZZ1(wt); if (ptn[i] <= level) ++wt; } if (invararg > n || invararg == 0) dlim = n; else dlim = invararg+1; success = FALSE; for (cell1 = 0; cell1 < n; cell1 = cell2 + 1) { for (cell2 = cell1; ptn[cell2] > level; ++cell2) {} if (cell2 == cell1) continue; for (iv = cell1; iv <= cell2; ++iv) { liv = lab[iv]; QUEUE[0] = liv; DIST[liv] = 0; RESETMARKS1; MARK1(liv); longcode = 0; head = 0; tail = 1; while (tail < n && head < tail) { i = QUEUE[head++]; if (DIST[i] >= dlim) break; vi = v[i]; di = d[i]; for (j = 0; j < di; ++j) { k = e[vi+j]; if (ISNOTMARKED1(k)) { MARK1(k); DIST[k] = DIST[i] + 1; kcode = DIST[k]+CELLCODE[k]; ACCUM(longcode,FUZZ1(kcode)); QUEUE[tail++] = k; } } } invar[liv] = CLEANUP(longcode); if (invar[liv] != invar[lab[cell1]]) success = TRUE; } if (success) break; } } /***************************************************************************** * * * adjacencies_sg() assigns to each vertex v a code depending on which cells * * it is joined to and from, and how many times. It is intended to provide * * better partitioning that the normal refinement routine for digraphs. * * It will not help with undirected graphs in nauty at all. * * * *****************************************************************************/ 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) { int *d,*e; int vwt,wwt; int *ei,di,i; size_t *v,j; SG_VDE(g,v,d,e); #if !MAXN DYNALLOC1(int,work2,work2_sz,n,"adjacencies_sg"); #endif vwt = 1; for (i = 0; i < n; ++i) { work2[lab[i]] = vwt; if (ptn[i] <= level) ++vwt; invar[i] = 0; } for (i = 0; i < n; ++i) { vwt = FUZZ1(work2[i]); wwt = 0; di = d[i]; ei = e + v[i]; for (j = 0; j < di; ++j) { ACCUM(wwt,FUZZ2(work2[ei[j]])); ACCUM(invar[ei[j]],vwt); } ACCUM(invar[i],wwt); } } /***************************************************************************** * * * sparsenauty(g,lab,ptn,orbits,&options,&stats,h) * * is a slightly simplified interface to nauty(). It allocates enough * * workspace for 60 automorphisms and checks that the sparsegraph dispatch * * vector is in use. * * * *****************************************************************************/ void sparsenauty(sparsegraph *g, int *lab, int *ptn, int *orbits, optionblk *options, statsblk *stats, sparsegraph *h) { int m,n; if (options->dispatch != &dispatch_sparse) { fprintf(ERRFILE,"Error: sparsenauty() needs standard options block\n"); exit(1); } n = g->nv; m = SETWORDSNEEDED(n); #if !MAXN /* Don't increase 2*60*m in the following without also increasing the static decalaration of snwork[] above. */ DYNALLOC1(set,snwork,snwork_sz,2*60*m,"densenauty malloc"); #endif nauty((graph*)g,lab,ptn,NULL,orbits,options,stats, snwork,2*60*m,m,n,(graph*)h); } /***************************************************************************** * * * nausparse_check() checks that this file is compiled compatibly with the * * given parameters. If not, call exit(1). * * * *****************************************************************************/ void nausparse_check(int wordsize, int m, int n, int version) { if (wordsize != WORDSIZE) { fprintf(ERRFILE,"Error: WORDSIZE mismatch in nausparse.c\n"); exit(1); } #if MAXN if (m > MAXM) { fprintf(ERRFILE,"Error: MAXM inadequate in nausparse.c\n"); exit(1); } if (n > MAXN) { fprintf(ERRFILE,"Error: MAXN inadequate in nausparse.c\n"); exit(1); } #endif if (version < NAUTYREQUIRED) { fprintf(ERRFILE,"Error: nausparse.c version mismatch\n"); exit(1); } } /***************************************************************************** * * * nausparse_freedyn() - free the dynamic memory in this module * * * *****************************************************************************/ void nausparse_freedyn(void) { #if !MAXN DYNFREE(vmark1,vmark1_sz); DYNFREE(vmark2,vmark2_sz); DYNFREE(work1,work1_sz); DYNFREE(work2,work2_sz); DYNFREE(work3,work3_sz); DYNFREE(work4,work4_sz); DYNFREE(snwork,snwork_sz); #endif }