From 58acff54b1cd64cb23b9d0b1a304eb9db768e3eb Mon Sep 17 00:00:00 2001 From: Andrew Guschin Date: Sun, 13 Aug 2023 01:27:00 +0400 Subject: Initial commit --- nauty/nausparse.c | 1757 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1757 insertions(+) create mode 100644 nauty/nausparse.c (limited to 'nauty/nausparse.c') diff --git a/nauty/nausparse.c b/nauty/nausparse.c new file mode 100644 index 0000000..d3b30cc --- /dev/null +++ b/nauty/nausparse.c @@ -0,0 +1,1757 @@ +/***************************************************************************** +* * +* Sparse-graph-specific auxiliary source file for version 2.8 of nauty. * +* * +* Copyright (2004-2022) 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[] * +* 6-Apr-21 : increase work space in sparsenauty() * +* 16-Nov-22 : fix an error in the Traces utility comparelab_tr() * +* * +*****************************************************************************/ + +#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*500*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; + +#if !MAXN + DYNALLOC1(int,work1,work1_sz,n,"comparelab_tr"); +#endif +#define NGHCOUNTS work1 + + memset(NGHCOUNTS,0,n*sizeof(int)); + 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; + + mina = n; + for (j = 0; j < d1; ++j) { + (NGHCOUNTS[col[invlab1[e1[j]]]])++; + } + for (j = 0; j < d1; ++j) + { + k = col[invlab2[e2[j]]]; + if (NGHCOUNTS[k]) { + (NGHCOUNTS[k])--; + } + else { + if (k < mina) mina = k; + } + } + if (mina != n) + { + for (j = 0; j < d1; ++j) + { + k = col[invlab1[e1[j]]]; + if (NGHCOUNTS[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 500 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*500*m in the following without also increasing + the static decalaration of snwork[] above. */ + DYNALLOC1(set,snwork,snwork_sz,2*500*m,"densenauty malloc"); +#endif + + nauty((graph*)g,lab,ptn,NULL,orbits,options,stats, + snwork,2*500*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 +} -- cgit v1.2.3