From 08527faf160d498c5d540100b9bb2dbb0a0224d0 Mon Sep 17 00:00:00 2001 From: kaneko Date: Thu, 17 Jul 2025 13:42:08 +0900 Subject: [PATCH] feat(grcc): add QGRAF options This patch implements the following QGRAF options: onepi/onepr onshell/offshell nosigma/sigma nosnail/snail notadpole/tadpole simple/notsimple bipart/nonbipart cycli/cyclr floop The following options are not implemented: topol onevi/onevr onshellx/offshellx The default options are updated to generate the set of all connected graphs. Specifying additional options can generate a subset of the default output. Example: #define LOOPS "2" Vectors Q1,Q2,p1,...,p8; Set QQ:Q1,Q2; Set PP:p1,...,p8; Local F = topologies_(`LOOPS',2,{3,4},QQ,PP); Print +ss; .end This generates 24 terms, corresponding to 24 topologies, due to the updated default options. When applied to the 3-loop case, the same example yields 240 terms. Note: The above example crashes when called with "`LOOPS',-2,{3,4},QQ,PP". This behaviour was also seen in previous versions. --- sources/diawrap.cc | 407 +--------- sources/grcc.cc | 1884 ++++++++++++++++++++++++++++++------------- sources/grcc.h | 99 ++- sources/grccparam.h | 46 +- 4 files changed, 1461 insertions(+), 975 deletions(-) diff --git a/sources/diawrap.cc b/sources/diawrap.cc index c6d4ddb4..723acc9c 100644 --- a/sources/diawrap.cc +++ b/sources/diawrap.cc @@ -14,7 +14,7 @@ typedef struct ToPoTyPe { WORD *vertmax; Options *opt; int cldeg[MAXPOINTS], clnum[MAXPOINTS], clext[MAXPOINTS]; - int cmind[MAXLEGS+1],cmaxd[MAXLEGS+1]; + int cmind[MAXPOINTS], cmaxd[MAXPOINTS]; int ncl, nvert; } TOPOTYPE; @@ -199,21 +199,6 @@ int ReConvertParticle(Model *model,int grccnum) } // #] ReConvertParticle : -// #[ numParticle : - -int numParticle(MODEL *m,WORD n) -{ - int i; - for ( i = 0; i < m->nparticles; i++ ) { - if ( m->vertices[i]->particles[0].number == n ) return(i); - if ( m->vertices[i]->particles[1].number == n ) return(i); - } - MesPrint("numParticle: particle %d not found in model",n); - Terminate(-1); - return(-1); -} - -// #] numParticle : // #[ ProcessDiagram : void ProcessDiagram(EGraph *eg, void *ti) @@ -245,8 +230,7 @@ void ProcessDiagram(EGraph *eg, void *ti) // // Now get the nodes // - if ( ( info->flags & NONODES ) == 0 ) { - for ( i = 0; i < eg->nNodes; i++ ) { + for ( i = 0; i < eg->nNodes; i++ ) { // // node_(number,coupling,particle_1(momentum_1),...,particle_n(momentum_n)) // @@ -261,7 +245,7 @@ void ProcessDiagram(EGraph *eg, void *ti) // function for when we want to work with counterterms. // if ( !eg->isExternal(i) ) { - afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill) + afill = fill; *fill++ = 0; *fill++ = 0; FILLARG(fill) cfill = fill; *fill++ = 0; intr = eg->nodes[i]->intrct; for ( j = 0; j < model->interacts[intr]->nclist; j++ ) { @@ -289,10 +273,8 @@ void ProcessDiagram(EGraph *eg, void *ti) edge = eg->nodes[i]->edges[j]; vect = ABS(edge)-1; *fill++ = 6+FUNHEAD; - int a; - if ( edge < 0 ) { a = ReConvertParticle(model,-eg->edges[vect]->ptcl); } - else { a = ReConvertParticle(model,eg->edges[vect]->ptcl); } - *fill++ = a; + + *fill++ = ReConvertParticle(model,eg->edges[vect]->ptcl); // code of particle_j *fill++ = FUNHEAD+2; FILLFUN(fill) *fill++ = edge < 0 ? -MINVECTOR: -VECTOR; @@ -305,131 +287,26 @@ void ProcessDiagram(EGraph *eg, void *ti) *fill++ = 1; *fill++ = 1; *fill++ = 3; } startfill[1] = fill-startfill; - } - } - if ( ( info->flags & WITHEDGES ) == WITHEDGES ) { - for ( i = 0; i < eg->nEdges; i++ ) { - int n1 = eg->edges[i]->nodes[0]; - int n2 = eg->edges[i]->nodes[1]; -// int l1 = eg->edges[i]->nlegs[0]; -// int l2 = eg->edges[i]->nlegs[1]; - - startfill = fill; - *fill++ = EDGE; - *fill++ = 0; - FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge -// - *fill++ = ARGHEAD+FUNHEAD+6; - *fill++ = 0; - FILLARG(fill) - *fill++ = 6+FUNHEAD; - int a = ReConvertParticle(model,eg->edges[i]->ptcl); - *fill++ = a; - *fill++ = FUNHEAD+2; - FILLFUN(fill) - *fill++ = -VECTOR; - // Look up in set of internal momenta set - if ( i < info->numextern ) { // Look up in set of external momenta - *fill++ = SetElements[Sets[info->externalset].first+i]; - } - else { // Look up in set of internal momenta set - *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; - } - *fill++ = 1; *fill++ = 1; *fill++ = 3; -// - *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from - *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to - startfill[1] = fill - startfill; - } - } - if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) { - for ( i = 0; i < eg->econn->nblocks; i++ ) { - startfill = fill; - *fill++ = BLOCK; - *fill++ = 0; - FILLFUN(fill); - *fill++ = -SNUMBER; - *fill++ = i+1; - *fill++ = -SNUMBER; - *fill++ = eg->econn->blocks[i].loop; -// -// Now we have to make a list of all nodes inside this block -// - int bnodes[GRCC_MAXNODES], k; - WORD *argfill = fill, *funfill; - *fill++ = 0; *fill++ = 0; FILLARG(fill) - *fill++ = 0; - for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0; - for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) { - bnodes[eg->econn->blocks[i].edges[k][0]] = 1; - bnodes[eg->econn->blocks[i].edges[k][1]] = 1; - } - for ( k = 0; k < GRCC_MAXNODES; k++ ) { - if ( bnodes[k] == 0 ) continue; -// -// Now we put the node inside this argument. -// - funfill = fill; - *fill++ = NODEFUNCTION; - *fill++ = 0; - FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = k+1; - numlegs = eg->nodes[k]->deg; - for ( j = 0; j < numlegs; j++ ) { - edge = eg->nodes[k]->edges[j]; - vect = ABS(edge)-1; - *fill++ = -VECTOR; - if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta - *fill++ = SetElements[Sets[info->externalset].first+vect]; - } - else { // Look up in set of internal momenta set - *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; - } - } - funfill[1] = fill-funfill; - } - *fill++ = 1; *fill++ = 1; *fill++ = 3; - *argfill = fill - argfill; - argfill[ARGHEAD] = argfill[0] - ARGHEAD; - startfill[1] = fill-startfill; - } - } - if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { - for ( i = 0; i < eg->econn->nopic; i++ ) { - startfill = fill; - *fill++ = ONEPI; - *fill++ = 0; - FILLFUN(fill); - *fill++ = -SNUMBER; - *fill++ = i+1; - *fill++ = -ONEPI; - for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) { - *fill++ = -SNUMBER; - *fill++ = eg->econn->opics[i].nodes[j]+1; - } - startfill[1] = fill-startfill; - } } // // Topology counter. We have exaggerated a bit with the eye on the far future. // - if ( info->numtopo-1 < MAXPOSITIVE ) { + if ( info->numtopo < MAXPOSITIVE ) { *fill++ = TOPO; *fill++ = FUNHEAD+2; FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo-1); + *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo); } - else if ( info->numtopo-1 < FULLMAX-1 ) { + else if ( info->numtopo < FULLMAX-1 ) { *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+4; FILLFUN(fill) *fill++ = ARGHEAD+4; *fill++ = 0; FILLARG(fill) *fill++ = 4; - *fill++ = (WORD)((info->numtopo-1) & WORDMASK); + *fill++ = (WORD)(info->numtopo & WORDMASK); *fill++ = 1; *fill++ = 3; } else { // for now: science fiction *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+6; FILLFUN(fill) *fill++ = ARGHEAD+6; *fill++ = 0; FILLARG(fill) - *fill++ = 6; *fill++ = (WORD)((info->numtopo-1) >> BITSINWORD); - *fill++ = (WORD)((info->numtopo-1) & WORDMASK); + *fill++ = 6; *fill++ = (WORD)(info->numtopo >> BITSINWORD); + *fill++ = (WORD)(info->numtopo & WORDMASK); *fill++ = 0; *fill++ = 1; *fill++ = 5; } // @@ -452,54 +329,19 @@ void ProcessDiagram(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; -// MesPrint("<> %a",newterm[0],newterm); - Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; } // #] ProcessDiagram : -// #[ fendMG : - -Bool fendMG(EGraph *eg, void *ti) -{ - DUMMYUSE(eg); - DUMMYUSE(ti); - return True; -} - -// #] fendMG : // #[ ProcessTopology : Bool ProcessTopology(EGraph *eg, void *ti) { // // This routine is called for each new topology. -// New convention: return True; generate the corresponding diagrams if needed -// return False; skip diagram generation (when asked for). // TERMINFO *info = (TERMINFO *)ti; -#ifdef WITHEARLYVETO - if ( ( ( info->flags & CHECKEXTERN ) == CHECKEXTERN ) && info->currentMODEL != NULL ) { - int i, j; - int numlegs, vect, edge; - for ( i = 0; i < eg->nNodes; i++ ) { - if ( eg->isExternal(i) ) continue; - numlegs = eg->nodes[i]->deg; - for ( j = 0; j < numlegs; j++ ) { - edge = eg->nodes[i]->edges[j]; - vect = ABS(edge)-1; - if ( vect < info->numextern && info->legcouple[vect][numlegs] == 0 ) { - -// This cannot be. - - info->numtopo++; - return False; - } - } - } - } -#endif if ( ( info->flags & TOPOLOGIESONLY ) == 0 ) { info->numtopo++; return True; @@ -513,8 +355,6 @@ Bool ProcessTopology(EGraph *eg, void *ti) WORD *tail = tdia + tdia[1]; WORD *tend = term + *term; WORD *fill, *startfill; - Model *model = (Model *)info->currentModel; - MODEL *m = (MODEL *)info->currentMODEL; int i, j; int numlegs, vect, edge; @@ -534,26 +374,6 @@ Bool ProcessTopology(EGraph *eg, void *ti) *fill++ = 0; FILLFUN(fill) *fill++ = -SNUMBER; *fill++ = i+1; - if ( model != NULL && m != NULL ) { - if ( !eg->isExternal(i) ) { - WORD *afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill) - WORD *cfill = fill; *fill++ = 0; - int cpl = 2*eg->nodes[i]->extloop+eg->nodes[i]->deg-2; - *fill++ = SYMBOL; *fill++ = 4; - *fill++ = m->couplings[0]; - *fill++ = cpl; - *fill++ = 1; *fill++ = 1; *fill++ = 3; - *cfill = fill - cfill; - *afill = fill - afill; - if ( *afill == ARGHEAD+8 && afill[ARGHEAD+4] == 1 ) { - fill = afill; *fill++ = -SYMBOL; *fill++ = afill[ARGHEAD+3]; - } - } - else { - *fill++ = -SNUMBER; - *fill++ = 1; - } - } // // Now the momenta. // @@ -570,116 +390,6 @@ Bool ProcessTopology(EGraph *eg, void *ti) } startfill[1] = fill-startfill; } - if ( ( info->flags & WITHEDGES ) == WITHEDGES ) { - for ( i = 0; i < eg->nEdges; i++ ) { - int n1 = eg->edges[i]->nodes[0]; - int n2 = eg->edges[i]->nodes[1]; -// int l1 = eg->edges[i]->nlegs[0]; -// int l2 = eg->edges[i]->nlegs[1]; - startfill = fill; - *fill++ = EDGE; - *fill++ = 0; - FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge -// - *fill++ = -VECTOR; - if ( i < info->numextern ) { // Look up in set of external momenta - *fill++ = SetElements[Sets[info->externalset].first+i]; - } - else { // Look up in set of internal momenta set - *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; - } -// - *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from - *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to - startfill[1] = fill - startfill; - } - } -// -// Block information -// - if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) { - for ( i = 0; i < eg->econn->nblocks; i++ ) { - startfill = fill; - *fill++ = BLOCK; - *fill++ = 0; - FILLFUN(fill); - *fill++ = -SNUMBER; - *fill++ = i+1; - *fill++ = -SNUMBER; - *fill++ = eg->econn->blocks[i].loop; -// -// Now we have to make a list of all nodes inside this block -// - int bnodes[GRCC_MAXNODES], k; - WORD *argfill = fill, *funfill; - *fill++ = 0; *fill++ = 0; FILLARG(fill) - *fill++ = 0; - for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0; - for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) { - bnodes[eg->econn->blocks[i].edges[k][0]] = 1; - bnodes[eg->econn->blocks[i].edges[k][1]] = 1; - } - for ( k = 0; k < GRCC_MAXNODES; k++ ) { - if ( bnodes[k] == 0 ) continue; -// -// Now we put the node inside this argument. -// - funfill = fill; - *fill++ = NODEFUNCTION; - *fill++ = 0; - FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = k+1; - numlegs = eg->nodes[k]->deg; - for ( j = 0; j < numlegs; j++ ) { - edge = eg->nodes[k]->edges[j]; - vect = ABS(edge)-1; - *fill++ = -VECTOR; - if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta - *fill++ = SetElements[Sets[info->externalset].first+vect]; - } - else { // Look up in set of internal momenta set - *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; - } - } - funfill[1] = fill-funfill; - } - *fill++ = 1; *fill++ = 1; *fill++ = 3; - *argfill = fill - argfill; - argfill[ARGHEAD] = argfill[0] - ARGHEAD; - startfill[1] = fill-startfill; - } - -// if ( eg->econn->narticuls > 0 ) { -// startfill = fill; -// *fill++ = BLOCK; -// *fill++ = 0; -// FILLFUN(fill); -// for ( i = 0; i < eg->econn->snodes; i++ ) { -// if ( eg->econn->articuls[i] != 0 ) { -// *fill++ = -SNUMBER; -// *fill++ = i+1; -// } -// } -// startfill[1] = fill-startfill; -// } - } - if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { - for ( i = 0; i < eg->econn->nopic; i++ ) { - startfill = fill; - *fill++ = ONEPI; - *fill++ = 0; - FILLFUN(fill); - *fill++ = -SNUMBER; - *fill++ = i+1; - *fill++ = -ONEPI; - for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) { - *fill++ = -SNUMBER; - *fill++ = eg->econn->opics[i].nodes[j]+1; - } - startfill[1] = fill-startfill; - } - } // // Topology counter. We have exaggerated a bit with the eye on the far future. // @@ -718,15 +428,13 @@ Bool ProcessTopology(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; -//MesPrint("<> %a",*newterm,newterm); - Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; info->numtopo++; - return False; + return True; } -// #] ProcessTopology : +// #] ProcessTopology : // #[ GenDiagrams : int GenDiagrams(PHEAD WORD *term, WORD level) @@ -739,7 +447,7 @@ int GenDiagrams(PHEAD WORD *term, WORD level) int babble = 0; // Later we may set this at the FORM code level TERMINFO info; WORD inset,outset,*coupl,setnum,optionnumber = 0; - int i, j, cpl[GRCC_MAXNCPLG]; + int i, cpl[GRCC_MAXNCPLG]; int ninitl, initlPart[GRCC_MAXLEGS], nfinal, finalPart[GRCC_MAXLEGS]; for ( i = 0; i < GRCC_MAXNCPLG; i++ ) cpl[i] = 0; // @@ -751,7 +459,6 @@ int GenDiagrams(PHEAD WORD *term, WORD level) info.diaoffset = AR.funoffset; info.externalset = term[info.diaoffset+FUNHEAD+7]; info.internalset = term[info.diaoffset+FUNHEAD+9]; - info.flags = 0; inset = term[info.diaoffset+FUNHEAD+3]; outset = term[info.diaoffset+FUNHEAD+5]; coupl = term + info.diaoffset + FUNHEAD + 10; @@ -764,6 +471,7 @@ int GenDiagrams(PHEAD WORD *term, WORD level) if ( term[info.diaoffset+1] > *coupl+FUNHEAD+10 ) optionnumber = term[info.diaoffset+*coupl+FUNHEAD+11]; } + setnum = term[info.diaoffset+FUNHEAD+1]; m = AC.models[SetElements[Sets[setnum].first]]; @@ -778,59 +486,33 @@ int GenDiagrams(PHEAD WORD *term, WORD level) opt = new Options(); - opt->setOutAG(ProcessDiagram, &info); - opt->setOutMG(ProcessTopology, &info); -// opt->setEndMG(fendMG, &info); + opt->setOutAG(ProcessDiagram, (void*) &info); + opt->setOutMG(ProcessTopology, (void*) &info); opt->values[GRCC_OPT_1PI] = ( optionnumber & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; opt->values[GRCC_OPT_NoTadpole] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; -// -// Next are snails: -// opt->values[GRCC_OPT_No1PtBlock] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; -// - if ( ( optionnumber & WITHINSERTIONS ) == WITHINSERTIONS ) { - opt->values[GRCC_OPT_No2PtL1PI] = True; - opt->values[GRCC_OPT_NoAdj2PtV] = True; - opt->values[GRCC_OPT_No2PtL1PI] = True; - } - else { - opt->values[GRCC_OPT_NoAdj2PtV] = True; - } + opt->values[GRCC_OPT_No2PtL1PI] = ( optionnumber & WITHINSERTIONS ) == WITHINSERTIONS; opt->values[GRCC_OPT_SymmInitial] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; - opt->values[GRCC_OPT_SymmFinal] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; - -// opt->values[GRCC_OPT_Block] = ( optionnumber & WITHBLOCKS ) == WITHBLOCKS; + if ( ( optionnumber & TOPOLOGIESONLY ) == TOPOLOGIESONLY ) + info.flags |= TOPOLOGIESONLY; opt->setOutputF(False,""); - opt->setOutputP(False,""); opt->printLevel(babble); -// opt->values[GRCC_OPT_Step] = GRCC_AGraph; - // Load the various arrays. ninitl = Sets[inset].last - Sets[inset].first; for ( i = 0; i < ninitl; i++ ) { x = SetElements[Sets[inset].first+i]; initlPart[i] = ConvertParticle(model,x); - info.legcouple[i] = m->vertices[numParticle(m,x)]->couplings; } nfinal = Sets[outset].last - Sets[outset].first; for ( i = 0; i < nfinal; i++ ) { x = SetElements[Sets[outset].first+i]; finalPart[i] = ConvertParticle(model,x); - info.legcouple[i+ninitl] = m->vertices[numParticle(m,x)]->couplings; } info.numextern = ninitl + nfinal; - for ( i = 2; i <= MAXLEGS; i++ ) { - if ( m->legcouple[i] == 1 ) { - for ( j = 0; j < info.numextern; j++ ) { - if ( info.legcouple[j][i] == 0 ) { info.flags |= CHECKEXTERN; goto Go_on; } - } - } - } -Go_on:; // // Now we have to sort out the coupling constants. // The argument at coupl can be of type -SNUMBER, -SYMBOL or generic @@ -847,20 +529,11 @@ Go_on:; int nc = coupl[1]*2 + ninitl + nfinal - 2; int *scratch = (int *)Malloc1(nc*sizeof(int),"DistrN"); scratch[0] = -2; // indicating startup cq first call. - - if ( ( info.flags & TOPOLOGIESONLY ) == 0 ) { - while ( DistrN(nc,cpl,m->ncouplings,scratch) ) { - proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); - delete proc; - info.numtopo = 1; - } - } - else { - cpl[0] = nc; + while ( DistrN(nc,cpl,m->ncouplings,scratch) ) { proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); + ninitl, initlPart, nfinal, finalPart, cpl); delete proc; + info.numtopo = 1; } M_free(scratch,"DistrN"); opt->end(); @@ -889,9 +562,6 @@ Go_on:; t += 2; } } -/* - And now the generation: -*/ proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl); opt->end(); @@ -926,11 +596,15 @@ int processVertex(TOPOTYPE *TopoInf, int pointsremaining, int level) TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level]; TopoInf->clnum[TopoInf->ncl] = j; TopoInf->clext[TopoInf->ncl] = 0; + TopoInf->cmind[TopoInf->ncl] = 0; + TopoInf->cmaxd[TopoInf->ncl] = 0; TopoInf->ncl++; - MGraph *mgraph = new MGraph(1, TopoInf->ncl, TopoInf->cldeg, - TopoInf->clnum, TopoInf->clext, - TopoInf->cmind, TopoInf->cmaxd, TopoInf->opt); + MGraph *mgraph = new MGraph(1, + TopoInf->ncl, TopoInf->cldeg, + TopoInf->clnum, TopoInf->clext, + TopoInf->cmind, TopoInf->cmaxd, + TopoInf->opt); mgraph->generate(); @@ -989,8 +663,6 @@ int GenTopologies(PHEAD WORD *term, WORD level) info.numextern = nlegs; - for ( i = 0; i <= MAXLEGS; i++ ) { TopoInf.cmind[i] = TopoInf.cmaxd[i] = 0; } - t1 += 10; if ( t1 < tstop && t1[0] == -SETSET ) { TopoInf.vertmax = &(SetElements[Sets[t1[1]].first]); @@ -1002,27 +674,19 @@ int GenTopologies(PHEAD WORD *term, WORD level) if ( t1 < tstop && t1[0] == -SNUMBER ) { if ( ( t1[1] & NONODES ) == NONODES ) info.flags |= NONODES; if ( ( t1[1] & WITHEDGES ) == WITHEDGES ) info.flags |= WITHEDGES; - if ( ( t1[1] & WITHBLOCKS ) == WITHBLOCKS ) info.flags |= WITHBLOCKS; - if ( ( t1[1] & WITHONEPI ) == WITHONEPI ) info.flags |= WITHONEPI; opt->values[GRCC_OPT_1PI] = ( t1[1] & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; -// opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; - opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOSNAILS ) == NOSNAILS; + opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; opt->values[GRCC_OPT_No1PtBlock] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; - opt->values[GRCC_OPT_NoExtSelf] = ( t1[1] & NOEXTSELF ) == NOEXTSELF; - - if ( ( t1[1] & WITHINSERTIONS ) == WITHINSERTIONS ) { - opt->values[GRCC_OPT_No2PtL1PI] = True; - opt->values[GRCC_OPT_NoAdj2PtV] = True; - opt->values[GRCC_OPT_No2PtL1PI] = True; - } + opt->values[GRCC_OPT_No2PtL1PI] = ( t1[1] & WITHINSERTIONS ) == WITHINSERTIONS; opt->values[GRCC_OPT_SymmInitial] = ( t1[1] & WITHSYMMETRIZE ) == WITHSYMMETRIZE; } info.numdia = 0; info.numtopo = 1; - opt->setOutAG(ProcessDiagram, &info); - opt->setOutMG(ProcessTopology, &info); + opt->setOutAG(ProcessDiagram, (void*) &info); + opt->setOutMG(ProcessTopology, (void*) &info); + // // Now we should sum over all possible vertices and run MGraph for // each combination. This is done by recursion in the processVertex routine @@ -1037,6 +701,7 @@ int GenTopologies(PHEAD WORD *term, WORD level) } for ( i = 0; i < nlegs; i++ ) { TopoInf.cldeg[i] = 1; TopoInf.clnum[i] = 1; TopoInf.clext[i] = -1; + TopoInf.cmind[i] = 0; TopoInf.cmaxd[i] = 0; } int points = 2*nloops-2+nlegs; diff --git a/sources/grcc.cc b/sources/grcc.cc index f69d9a5f..8479f55f 100644 --- a/sources/grcc.cc +++ b/sources/grcc.cc @@ -16,8 +16,20 @@ extern "C" { #include #include -// #define DEBUG9 +//============================================================== +// new default values of options +// #define OLDDEFAULT + +// new function for the output in python format +// #define NEWOUTPY + +//============================================================== +// For debugging +// +// #define MONITOR // in gaph elimination +// #define CHECK // check consistency +//============================================================== // optimization : best combination depends on process by process #define SIMPSEL @@ -45,19 +57,26 @@ extern "C" { using namespace Grcc; #endif +#define MAXSTR 1024 + //-------------------------------------------------------------- // Macro functions #define CLWIGHTD(x) (5*(x)) #define CLWIGHTO(x) (3*(x)-2) +#define MAXSTRLEN 81 + +#define MASK(n) ((1ul) << (n)) + +#define BOOLSTR(x) ((x) ? "True" : "False") //-------------------------------------------------------------- // Static variables -static OptDef optDef[] = { +static OptDef optDef0[] = { {"Step", "Generate particle assigned graphs", GRCC_AGraph, 0}, {"Outgrf", "Output to file (out.grf)", False, 0}, {"Outgrp", "Output to file (out.grp)", False, 0}, - {"1PI", "Only 1PI graphs", True , 0}, + {"OPI", "Only 1PI graphs", True , 0}, {"NoSelfLoop", "Exclude graphs with loops consist of 1 edge", True , 0}, {"NoTadpole", "Exclude graphs with tadpoles (2 edge connected)", True , 0}, {"No1PtBlock", "Exclude graphs with tadpole blocks (2 node connected)", False, 0}, @@ -65,11 +84,51 @@ static OptDef optDef[] = { {"NoExtSelf", "Exclude graphs with 2-pt subgraphs at ext. particles", False, 0}, {"NoAdj2PtV", "Exclude graphs with an edge connecting 2-pt vertices", False, 0}, {"Block", "Exclude graphs with more than one block", False, 0}, + {"NoMultiEdge","Exclude graphs with multi-edges", False, 0}, {"SymmInitial","Symmetrize initial particles", False, 0}, {"SymmFinal", "Symmetrize final particles", False, 0}, }; +static OptDef optDef1[] = { + {"Step", "Generate particle assigned graphs", GRCC_AGraph, 0}, + {"Outgrf", "Ouput to file (out.grf)", False, 0}, + {"Outgrp", "Ouput to file (out.grp)", False, 0}, + {"OPI", "Only 1PI graphs", False, 0}, + {"NoSelfLoop", "Exclude graphs with loops consist of 1 edge", False, 0}, + {"NoTadpole", "Exclude graphs with tadpoles (2 edge connected)", False, 0}, + {"No1PtBlock", "Exclude graphs with tadpole blocks (2 node connected)", False, 0}, + {"No2PtL1PI", "Exclude graphs with 2-point subgraphs", False, 0}, + {"NoExtSelf", "Exclude graphs with 2-pt subgraphs at ext. particles", False, 0}, + {"NoAdj2PtV", "Exclude graphs with an edge connecting 2-pt vertices", False, 0}, + {"Block", "Exclude graphs with more than one block", False, 0}, + {"NoMultiEdge","Exclude graphs with multi-edges", False, 0}, + {"SymmInitial","Symmetrize initial particls", False, 0}, + {"SymmFinal", "Symmetrize final particls", False, 0}, +}; +#ifdef OLDDEFAULT +static OptDef *optDef = &(optDef0[0]); +static int nOptDef = sizeof(optDef0)/sizeof(OptDef); +#else +static OptDef *optDef = &(optDef1[0]); +static int nOptDef = sizeof(optDef1)/sizeof(OptDef); +#endif + + +static OptQGDef optQGDef[] = { + {"onepi", "onepr", "one-particle illreducible"}, + {"onshell", "offshell", "without self-energy part at external particles"}, + {"nosigma", "sigma", "no 2-point functions"}, + {"nosnail", "snail", "without snail"}, + {"notadpole", "tadpole", "without tadpole"}, + {"simple", "notsimple", "without self-loop nor multi-edge"}, + {"bipart", "nonbipart", "only bipartite graph"}, + {"cycli", "cyclr", "???"}, + {"floop", "", "without fermion loops of odd length"}, +#ifdef GRCC_QGRAF_OPT_TOPOL + {"topol", "?", "topology"}, +#endif +}; -static int nOptDef = sizeof(optDef)/sizeof(OptDef); +static int nOptQGDef = sizeof(optQGDef)/sizeof(OptQGDef); static int prlevel = 2; static ErExit *erExit = NULL; @@ -118,8 +177,11 @@ Options::Options(void) { if (nOptDef != GRCC_OPT_Size) { fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); + fprintf(GRCC_Stderr, "nOptDef=%d, GRCC_OPT_Size=%d\n", + nOptDef, GRCC_OPT_Size); exit(1); } + model = NULL; proc = NULL; sproc = NULL; @@ -132,7 +194,33 @@ Options::Options(void) argemg = NULL; argag = NULL; - setDefaultValue(); + setDefaultValues(); + + // QGraf options + if (nOptQGDef != GRCC_QGRAF_OPT_Size) { + fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); + fprintf(GRCC_Stderr, "nOptQGDef=%d, GRCC_QGRAF_OPT_Size=%d\n", + nOptQGDef, GRCC_QGRAF_OPT_Size); + exit(1); + } + nqgopt = 0; + for (int j = 0; j < nOptQGDef; j++) { + qgref[nqgopt].name = optQGDef[j].name; + qgref[nqgopt].index = j; + qgref[nqgopt].sign = +1; + nqgopt++; + if (strlen(optQGDef[j].cname) > 0) { + qgref[nqgopt].name = optQGDef[j].cname; + qgref[nqgopt].index = j; + qgref[nqgopt].sign = -1; + nqgopt++; + } + } + + // defulat values + for (int j = 0; j < GRCC_QGRAF_OPT_Size; j++) { + qgopt[j] = 0; + } // for output out = new Output(this); @@ -158,7 +246,7 @@ Options::~Options(void) } //-------------------------------------------------------------- -void Options::setDefaultValue(void) +void Options::setDefaultValues(void) { int j; @@ -173,6 +261,21 @@ void Options::setDefaultValue(void) prlevel = 1; } +//-------------------------------------------------------------- +void Options::setOldDefaultValues(void) +{ + int j; + + // defualt values + for (j = 0; j < GRCC_OPT_Size; j++) { + values[j] = optDef0[j].defaultv; + } + + values[GRCC_OPT_Step] = GRCC_AGraph; + + // print level + prlevel = 1; +} //-------------------------------------------------------------- void Options::setOutMG(OutEGB *omg, void *pt) @@ -237,6 +340,81 @@ int Options::getValue(int ind) return -1; } +//-------------------------------------------------------------- +void Options::setQGrafOpt(int *qg) +{ + // reset options : default is to generate all connected graphs + +#ifdef OLDDEFAULT + values[GRCC_OPT_1PI] = False; + values[GRCC_OPT_NoSelfLoop] = False; + values[GRCC_OPT_NoTadpole] = False; + values[GRCC_OPT_No1PtBlock] = False; + values[GRCC_OPT_No2PtL1PI] = False; + values[GRCC_OPT_NoExtSelf] = False; + values[GRCC_OPT_NoAdj2PtV] = False; + values[GRCC_OPT_Block] = False; + values[GRCC_OPT_SymmInitial] = False; + values[GRCC_OPT_SymmFinal] = False; +#endif + + for (int j = 0; j < GRCC_QGRAF_OPT_Size; j++) { + qgopt[j] = qg[j]; + } + + // {"onepi", "onepr"}, + if (qgopt[GRCC_QGRAF_OPT_ONEPI] > 0) { + values[GRCC_OPT_1PI] = True; + } else if (qgopt[GRCC_QGRAF_OPT_ONEPI] < 0) { + values[GRCC_OPT_1PI] = False; + } + + // {"onshell", "offshell"} + if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) { + values[GRCC_OPT_NoExtSelf] = 1; + } else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) { + values[GRCC_OPT_NoExtSelf] = 0; + } + + // {"nosigma", "sigma"} + + // {"nosnail", "snail"}, + if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] > 0) { + values[GRCC_OPT_NoTadpole] = True; + values[GRCC_OPT_No1PtBlock] = True; + } else if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] < 0) { + values[GRCC_OPT_NoTadpole] = False; + values[GRCC_OPT_No1PtBlock] = False; + } + + // {"notadpole", "tadpole"} + if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) { + values[GRCC_OPT_NoTadpole] = True; + } else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) { + values[GRCC_OPT_NoTadpole] = False; + } + + // {"simple", "notsimple"}, + if (qgopt[GRCC_QGRAF_OPT_SIMPLE] > 0) { + values[GRCC_OPT_NoSelfLoop] = True; + values[GRCC_OPT_NoMultiEdge] = True; + } else if (qgopt[GRCC_QGRAF_OPT_SIMPLE] < 0) { + values[GRCC_OPT_NoSelfLoop] = False; + values[GRCC_OPT_NoMultiEdge] = False; + } + + // {"bipart", "nonbipart"}, + // {"cycli", "cyclr"}, + +#ifdef GRCC_QGRAF_OPT_TOPOL + // {"topol", ""}, + if (qgopt[GRCC_QGRAF_OPT_TOPOL] > 0) { + values[GRCC_OPT_SymmInitial] = True; + values[GRCC_OPT_SymmFinal] = True; + } +#endif +} + //-------------------------------------------------------------- void Options::print(void) { @@ -251,6 +429,10 @@ void Options::print(void) j, optDef[j].name, values[j], optDef[j].defaultv); } printf(" outgrf=%s, outgrp=%s\n", out->outgrf, out->outgrp); + printf(" GRCC_QGRAF_OPT_Size=%d:\n", GRCC_QGRAF_OPT_Size); + for (j=0; j < GRCC_QGRAF_OPT_Size; j++) { + printf(" %4d %-10s = %2d\n", j, optQGDef[j].name, qgopt[j]); + } } //-------------------------------------------------------------- @@ -259,6 +441,17 @@ const OptDef *Options::getDef(void) return optDef; } +//-------------------------------------------------------------- +const OptDef *Options::getOldDef(void) +{ + return optDef0; +} + +//-------------------------------------------------------------- +const OptQGDef *Options::getQGDef(void) +{ + return optQGDef; +} //-------------------------------------------------------------- void Options::setOutputF(Bool outf, const char *fname) @@ -318,7 +511,7 @@ void Options::begin(Model *mdl) //-------------------------------------------------------------- void Options::end(void) { - if (prlevel > 0) { + if (prlevel > 1) { printf("Optimization: "); #ifdef SIMPSEL printf("SIMPSEL=1 "); @@ -527,7 +720,6 @@ void Options::endSubProc(void) printf("* discarded for refinement: %ld\n", mgraph->discardRefine); printf("* discarded for disconnected: %ld\n", mgraph->discardDisc); printf("* discarded for duplication: %ld\n", mgraph->discardIso); - printf("* discarded by outMG: %ld\n", mgraph->discardMG); #endif } @@ -547,13 +739,6 @@ void Options::newMGraph(MGraph *mgr) { MGraph *mgraph = mgr; -#ifdef DEBUG - if (proc != NULL) { - printf("+++ New EGraph (MG): %ld\n", proc->mgrcount); - } else if (sproc != NULL) { - printf("+++ New EGraph (MG): %ld\n", sproc->mgrcount); - } -#endif if (proc != NULL) { proc->mgrcount++; mgr->egraph->mId = proc->mgrcount; @@ -584,15 +769,6 @@ void Options::newAGraph(EGraph *egraph) { Fraction sf, zr; -#ifdef DEBUG - if (proc != NULL) { - printf("+++ New EGraph (AG): %ld (%ld)\n", - proc->agrcount, proc->mgrcount); - } else if (sproc != NULL) { - printf("+++ New EGraph (AG): %ld (%ld)\n", - sproc->agrcount, sproc->mgrcount); - } -#endif if (proc != NULL) { proc->agrcount++; egraph->sId = proc->agrcount; @@ -622,11 +798,6 @@ void Options::newAGraph(EGraph *egraph) out->outEGraphP(egraph); } if (outag != NULL) { -#ifdef DEBUG1 - printf("call outag\n"); - egraph->model->prModel(); - egraph->print(); -#endif (*outag)(egraph, argag); } @@ -688,6 +859,7 @@ void Output::setOutgrf(const char *fname) } if (outgrf != NULL) { + // delete outgrf; free(outgrf); } if (fname == NULL || strlen(fname) < 1) { @@ -717,6 +889,10 @@ void Output::setOutgrp(const char *fname) //-------------------------------------------------------------- Bool Output::outBeginF(Model *mdl, Bool pr) { + time_t tp; + struct tm *tm; + char sdate[MAXSTRLEN]; + model = mdl; if (outgrf == NULL || strlen(outgrf) < 1) { @@ -739,8 +915,12 @@ Bool Output::outBeginF(Model *mdl, Bool pr) } return False; } + tp = time(NULL); + tm = localtime(&tp); + strftime(sdate, MAXSTRLEN, "%Y/%m/%d %H:%M:%S", tm); + fprintf(outgrfp, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"); - fprintf(outgrfp, "%% Generated by 'grcc'\n"); + fprintf(outgrfp, "%% Generated by grcc at \"%s\"\n", sdate); fprintf(outgrfp, "Version={2,2,0,0};\n"); if (model != NULL) { fprintf(outgrfp, "Model=\"./%s.mdl\";\n", model->name); @@ -755,6 +935,10 @@ Bool Output::outBeginF(Model *mdl, Bool pr) //-------------------------------------------------------------- Bool Output::outBeginP(Model *mdl, Bool pr) { + time_t tp; + struct tm *tm; + char sdate[MAXSTRLEN]; + model = mdl; if (outgrp == NULL || strlen(outgrp) < 1) { @@ -777,8 +961,12 @@ Bool Output::outBeginP(Model *mdl, Bool pr) } return False; } + tp = time(NULL); + tm = localtime(&tp); + strftime(sdate, MAXSTRLEN, "%Y/%m/%d %H:%M:%S", tm); + fprintf(outgrpp, "################################\n"); - fprintf(outgrpp, "# Generated by 'grcc'\n"); + fprintf(outgrpp, "# Generated by 'grcc' at \"%s\"\n", sdate); return True; } @@ -879,7 +1067,7 @@ void Output::outProcBeginP(Process *prc) } outproc = True; - fprintf(outgrfp, "# Process=%d;\n", proc->id); + fprintf(outgrpp, "# Process=%d;\n", proc->id); } //-------------------------------------------------------------- @@ -900,11 +1088,16 @@ void Output::outProcBegin0(int next, int couple, int loop) fprintf(outgrfp, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"); fprintf(outgrfp, "Process=%d;\n", procId); fprintf(outgrfp, "External=%d;\n", next); - for (k = 0, ex = 0; k < next; k++, ex++) { + int ninit = (next + 1) / 2; + for (k = 0, ex = 0; k < ninit; k++, ex++) { fprintf(outgrfp, "%4d= initial undef;\n", ex); } + for (k = 0; k < next - ninit; k++, ex++) { + fprintf(outgrfp, "%4d= final undef;\n", ex); + } fprintf(outgrfp, "Eend;\n"); fprintf(outgrfp, "GRCC_PHI=%d; ", couple); + // fprintf(outgrfp, "PHI=%d; ", couple); fprintf(outgrfp, "Loop=%d;\n", loop); fprintf(outgrfp, "OPI=%s;\n", (opt->values[GRCC_OPT_1PI] > 0? "Yes" : "No")); fprintf(outgrfp, "Assign=No;\n"); @@ -1014,9 +1207,6 @@ void Output::outProcEndP(void) if (outgrpp == NULL) { return; } - if (proc != NULL) { - fprintf(outgrpp, "# Pend=%d;\n", proc->id); - } fflush(outgrpp); } @@ -1028,13 +1218,6 @@ void Output::outEGraphF(EGraph *egraph) Bool popt; EFLine *fl; -#ifdef DEBUG9 - if (egraph->mgraph != NULL) { - printf("outEGraph:sId=%ld\n", egraph->mId); - // egraph->mgraph->print(); - egraph->mgraph->mconn->print(); - } -#endif if (outgrfp == NULL) { return; } @@ -1105,8 +1288,9 @@ void Output::outEGraphF(EGraph *egraph) if (lg != 0) { fprintf(outgrfp, ", "); } - ed = V2Iedge(egraph->nodes[nd]->edges[lg]); - fprintf(outgrfp, "%4d", egraph->nodes[nd]->edges[lg]); + ed = V2Iedge(egraph->nodes[nd]->edges[lg]); + // fprintf(outgrfp, "%4d", egraph->nodes[nd]->edges[lg]); + fprintf(outgrfp, "%4d", abs(egraph->nodes[nd]->edges[lg])); ptcl = egraph->edges[ed]->ptcl; if (mdl != NULL && ptcl != 0 && egraph->assigned) { if (egraph->nodes[nd]->edges[lg] < 0) { @@ -1122,17 +1306,17 @@ void Output::outEGraphF(EGraph *egraph) fprintf(outgrfp, "};\n"); } // print Fermion line as comment line - if (egraph->nflines >= 0) { + if (egraph->nFlines >= 0) { fprintf(outgrfp, "%% FLines=%d; FSign=%d; sId=%ld;\n", - egraph->nflines, egraph->fsign, egraph->sId); + egraph->nFlines, egraph->fsign, egraph->sId); } - for (j = 0; j < egraph->nflines; j++) { + for (j = 0; j < egraph->nFlines; j++) { fl = egraph->flines[j]; fprintf(outgrfp, "%% %4d", j); if (fl->ftype == FL_Open) { - fprintf(outgrfp, "[Open]=["); + fprintf(outgrfp, "[FOpen]=["); } else if (fl->ftype == FL_Closed) { - fprintf(outgrfp, "[Loop]=["); + fprintf(outgrfp, "[FLoop]=["); } else { fprintf(outgrfp, " ?%d", fl->ftype); } @@ -1140,7 +1324,7 @@ void Output::outEGraphF(EGraph *egraph) if (k != 0) { fprintf(outgrfp, ", "); } - fprintf(outgrfp, "%d", Abs(fl->elist[k])); + fprintf(outgrfp, "%d", fl->elist[k]); } fprintf(outgrfp, "];\n"); } @@ -1149,72 +1333,72 @@ void Output::outEGraphF(EGraph *egraph) } //-------------------------------------------------------------- +#ifdef NEWOUTPY +void Output::outEGraphP(EGraph *eg) +{ + if (outgrpp == NULL) { + return; + } + + if (eg->assigned) { + printf("outEGraphP:egraph:egraph[%ld]\n", eg->mId-1); + eg->printPy(outgrpp, eg->mId); + return; + } else { + printf("outEGraphP:mgraph:egraph[%ld]\n", eg->mId-1); + eg->mgraph->printPy(outgrpp, eg->mId); + return; + } +} +#else void Output::outEGraphP(EGraph *eg) { - int j, nd, lg, ed, nin, nfi; - char nl = '\n'; + int j, nd, lg, ed, k; ENode *node; static char undef[] = "Undef"; char *s; + EFLine *fl; - // model - fprintf(outgrpp, "egraph[%ld] = {%c", eg->sId-1, nl); - if (model==NULL) { - fprintf(outgrpp, " \"Model\": [\"Undef\", 1],%c", nl); - } else { - fprintf(outgrpp, " \"Model\": [\"%s\", %d],%c", - model->name, model->ncouple, nl); + if (outgrpp == NULL) { + return; } - // process - nin = nfi = 0; - for (j = 0; j < eg->nNodes; j++) { - if (eg->isExternal(j)) { - if (eg->nodes[j]->extloop == GRCC_AT_Final) { - nfi++; - } else { - nin++; - } - } - } - fprintf(outgrpp, " \"Process\": [[%d, %d], %d, [", - nin, nfi, eg->nLoops); - if (proc != NULL && model != NULL) { - for (j = 0; j < model->ncouple; j++) { - if (j != 0) { - fprintf(outgrpp, ", "); - } - fprintf(outgrpp, "%d", proc->clist[j]); - } + if (eg->assigned) { + fprintf(outgrpp, "egraph[%ld] = {\n", eg->sId-1); + } else { + fprintf(outgrpp, "egraph[%ld] = {\n", eg->mId-1); } - fprintf(outgrpp, "]],%c", nl); - // option of the process - fprintf(outgrpp, " \"Opt\": ["); - if (eg->opt != NULL) { - for (j = 0; j < GRCC_OPT_Size; j++) { - if (j != 0) { - fprintf(outgrpp, ", "); - } - fprintf(outgrpp, "%d", eg->opt->values[j]); - } + // process + if (proc != NULL) { + fprintf(outgrpp, " \"Process\": {\n"); + proc->outProcP(outgrpp); + fprintf(outgrpp, " },\n"); } - fprintf(outgrpp, "],%c", nl); // graph - fprintf(outgrpp, " \"Id\": [%ld, %ld, %ld],%c", - eg->mId, eg->aId, eg->sId, nl); - fprintf(outgrpp, " \"GStat\": [%d, %d, %d],%c", - eg->assigned, eg->nNodes, eg->nEdges, nl); - fprintf(outgrpp, " \"Sym\": [%d, %ld, %ld, %ld, %ld],%c", - eg->fsign, eg->nsym, eg->esym, eg->nsym1, eg->multp, nl); + fprintf(outgrpp, " \"GId\": [%ld, %ld, %ld],\n", + eg->mId, eg->aId, eg->sId); + fprintf(outgrpp, " \"GParam\": {\"Assigned\":%d, " + "\"NNodes\":%d, \"NEdges\":%d, \"NFLines\":%d},\n", + eg->assigned, eg->nNodes, eg->nEdges, eg->nFlines); + fprintf(outgrpp, " \"Sym\": {\"FSign\":%d, \"NSym\":%ld, " + "\"ESym\":%ld, \"NSym1\":%ld, \"Multp\":%ld},\n", + eg->fsign, eg->nsym, eg->esym, eg->nsym1, eg->multp); // nodes - fprintf(outgrpp, " \"Nodes\": [%c", nl); + fprintf(outgrpp, " \"Nodes\": [\n"); for (nd = 0; nd < eg->nNodes; nd++) { node = eg->nodes[nd]; - if (model==NULL || !eg->assigned) { + if (model==NULL) { s = undef; + } else if(!eg->assigned) { + if (eg->isExternal(nd)) { + ed = Abs(eg->nodes[nd]->edges[0])-1; + s = model->particleName(eg->edges[ed]->ptcl); + } else { + s = undef; + } } else if (eg->isExternal(nd)) { s = model->particleName(node->intrct); } else { @@ -1228,27 +1412,56 @@ void Output::outEGraphP(EGraph *eg) } fprintf(outgrpp, "%d", eg->nodes[nd]->edges[lg]); } - fprintf(outgrpp, "]],%c", nl); + fprintf(outgrpp, "]],\n"); } - fprintf(outgrpp, " ],%c", nl); + fprintf(outgrpp, " ],\n"); // edges - fprintf(outgrpp, " \"Edges\": [%c", nl); + fprintf(outgrpp, " \"Edges\": [\n"); for (ed = 0; ed < eg->nEdges; ed++) { if (model != NULL) { s = model->particleName(eg->edges[ed]->ptcl); - fprintf(outgrpp, " [%d, \"%s\", [%d, %d]],%c", - ed+1, s, eg->edges[ed]->nodes[0], eg->edges[ed]->nodes[1], nl); + fprintf(outgrpp, " [%d, \"%s\", [%d, %d]],\n", + ed+1, s, eg->edges[ed]->nodes[0], + eg->edges[ed]->nodes[1]); } else { - fprintf(outgrpp, " [%d, \"Undef\", [%d, %d]],%c", - ed+1, eg->edges[ed]->nodes[0], eg->edges[ed]->nodes[1], nl); + fprintf(outgrpp, " [%d, \"Undef\", [%d, %d]],\n", + ed+1, eg->edges[ed]->nodes[0], + eg->edges[ed]->nodes[1]); } } - fprintf(outgrpp, " ],%c", nl); + fprintf(outgrpp, " ],\n"); + if (eg->nFlines >= 0) { + // print Fermion line as comment line + fprintf(outgrpp, " \"FLines\": [\n"); + for (j = 0; j < eg->nFlines; j++) { + fl = eg->flines[j]; + fprintf(outgrpp, " [%d, ", j); + if (fl->ftype == FL_Open) { + fprintf(outgrpp, "\"FOpen\", "); + } else if (fl->ftype == FL_Closed) { + fprintf(outgrpp, "\"FLoop\", "); + } else { + fprintf(outgrpp, "\"Undef%d\", ", fl->ftype); + } + fprintf(outgrpp, "["); + for (k = 0; k < fl->nlist; k++) { + if (k != 0) { + fprintf(outgrpp, ", "); + } + fprintf(outgrpp, "%d", fl->elist[k]); + } + fprintf(outgrpp, "]],\n"); + } + fprintf(outgrpp, " ]\n"); + } + + // end of a graph fprintf(outgrpp, "};\n"); } +#endif //-------------------------------------------------------------- void Output::outModelF(void) @@ -1367,31 +1580,44 @@ void Output::outModelP(void) // copupling constants fprintf(mdlfp, "Model={\n"); fprintf(mdlfp, " \"Model\": ["); - fprintf(mdlfp, "\"%s\", %d, [", model->name, model->ncouple); + fprintf(mdlfp, "\"%s\", %d, [", + model->name, model->ncouple); for (j = 0; j < model->ncouple; j++) { if (j != 0) { fprintf(mdlfp, ", "); } - fprintf(mdlfp, "%s", model->cnlist[j]); + fprintf(mdlfp, "\"%s\"", model->cnlist[j]); + } + fprintf(mdlfp, "], "); + if (model->defpart == GRCC_DEFBYCODE) { + fprintf(mdlfp, "\"ByCode\""); + } else { + fprintf(mdlfp, "\"ByName\""); } fprintf(mdlfp, "],\n"); fprintf(mdlfp, " \"Particles\": [\n"); for (j = 1; j < model->nParticles; j++) { pt = model->particles[j]; - ptn = pt->typeGName(); - fprintf(mdlfp, " [\"%s\", \"%s\", \"%s\"],\n", - pt->name, pt->aname, ptn); + ptn = pt->typeName(); + fprintf(mdlfp, " [\"%s\", %d, \"%s\", %d, \"%s\", %d],\n", + pt->name, pt->pcode, pt->aname, pt->acode, + ptn, pt->extonly); } + fprintf(mdlfp, " ],\n"); fprintf(mdlfp, " \"Interactions\": [\n"); for (j = 0; j < model->nInteracts; j++) { vt = model->interacts[j]; - fprintf(mdlfp, " [\"%s\", %d, [", vt->name, vt->nplist); + fprintf(mdlfp, " [\"%s\", %d, %d, [", vt->name, vt->id, vt->nplist); for (k = 0; k < vt->nplist; k++) { if (k != 0) { fprintf(mdlfp, ", "); } - fprintf(mdlfp, "%s", model->particleName(vt->plist[k])); + if (model->defpart == GRCC_DEFBYCODE) { + fprintf(mdlfp, "%d", vt->plist[k]); + } else { + fprintf(mdlfp, "\"%s\"", model->particleName(vt->plist[k])); + } } fprintf(mdlfp, "], ["); for (k = 0; k < model->ncouple; k++) { @@ -1403,7 +1629,7 @@ void Output::outModelP(void) fprintf(mdlfp, "]],\n"); } fprintf(mdlfp, " ],\n"); - fprintf(mdlfp, "],\n"); + fprintf(mdlfp, "};\n"); fclose(mdlfp); delete[] mdlfn; } @@ -1571,6 +1797,7 @@ Interaction::Interaction(Model *modl, int iid, const char *nam, int icd, int *cp { static char buff[100]; static int nbuff=100; + static Bool prmsg = True; Particle *p; int j, ndir, sdir, jdir, nmaj, jmaj, ngho, sgho, jgho, ptcl; Bool ok; @@ -1672,14 +1899,16 @@ Interaction::Interaction(Model *modl, int iid, const char *nam, int icd, int *cp fprintf(GRCC_Stderr, "*** Interaction: odd number of Majorana particles\n"); ok = False; } - if (ndir + nmaj + ngho > 2) { - fprintf(GRCC_Stderr, "+++ Interaction: more than 2 Dirac/Majorana/Ghost " + if (ndir + nmaj + ngho > 2 && prmsg) { + fprintf(GRCC_Stderr, "+++ Interaction: more than 2 " + "Dirac/Majorana/Ghost " "particles are interacting.\n"); - fprintf(GRCC_Stderr, " Interaction has %d Dirac, %d Ghost and %d Majorana " - "particles.\n", ndir, nmaj, ngho); - fprintf(GRCC_Stderr, " Sign factors related to fermions will not " - "be calculated correctly.\n"); - mdl->skipFLine = True; + fprintf(GRCC_Stderr, " Interaction has %d Dirac, %d Ghost " + "and %d Majorana particles.\n", + ndir, nmaj, ngho); + fprintf(GRCC_Stderr, " Sign factors related to fermions may " + "be inconsistent with your calculation method.\n"); + prmsg = False; } if (!ok) { prInteraction(); @@ -1778,7 +2007,7 @@ Model::Model(MInput *minp) maxcpl = 0; maxloop = 0; ncplgcp = 0; - skipFLine = False; + cplgvl = NULL; addParticle(&pundef); } @@ -1788,13 +2017,15 @@ Model::~Model(void) { int j; - for (j = ncplgcp-1; j >= 0; j--) { - cplgvl[j] = delintdup(cplgvl[j]); + if (ncplgcp > 0) { + for (j = ncplgcp-1; j >= 0; j--) { + cplgvl[j] = delintdup(cplgvl[j]); + } + delete[] cplgvl; + cplgnvl = delintdup(cplgnvl); + cplglg = delintdup(cplglg); + cplgcp = delintdup(cplgcp); } - delete[] cplgvl; - cplgnvl = delintdup(cplgnvl); - cplglg = delintdup(cplglg); - cplgcp = delintdup(cplgcp); for (j = GRCC_MAXMINTERACT-1; j >= 0; j--) { if (interacts[j] != NULL) { @@ -2148,9 +2379,6 @@ void Model::addInteractionEnd(void) } } } -#ifdef DEBUG - prModel(); -#endif } //-------------------------------------------------------------- @@ -2316,6 +2544,56 @@ int Model::findInteractionCode(int icd) return -1; } +//-------------------------------------------------------------- +void Model::printMInput(MInput *min) +{ + printf(" \"Model\": [\"%s\", %d, [", + min->name, min->ncouple); + for (int j = 0; j < min->ncouple; j++) { + if (j != 0) { + printf(", "); + } + printf("\"%s\"", min->cnamlist[j]); + } + printf("], "); + if (min->defpart == GRCC_DEFBYCODE) { + printf("\"ByCode\""); + } else { + printf("\"ByName\""); + } + printf("],\n"); +} + +//-------------------------------------------------------------- +void Model::printPInput(PInput *pin) +{ + printf(" [\"%s\"(%d), \"%s\"(%d), \"%s\"(%d), %d],\n", + pin->name, pin->pcode, pin->aname, pin->acode, + pin->ptypen, pin->ptypec, pin->extonly); +} + +//-------------------------------------------------------------- +void Model::printIInput(IInput *iin) +{ + int j; + + printf(" [\"%s\"(%d), %d, [", iin->name, iin->icode, iin->nplistn); + for (j = 0; j < iin->nplistn; j++) { + if (j != 0) { + printf(", "); + } + printf("\"%s\"(%d)", iin->plistn[j], iin->plistc[j]); + } + printf("], ["); + for (j = 0; j < GRCC_MAXNCPLG; j++) { + if (j != 0) { + printf(", "); + } + printf("%d", iin->cvallist[j]); + } + printf("],\n"); +} + //************************************************************** // proc.cc @@ -2859,6 +3137,7 @@ SProcess::~SProcess(void) delete pnclass; delete astack; delete mgraph; + model = NULL; } //-------------------------------------------------------------- @@ -3037,13 +3316,14 @@ Process::Process(int pid, Model *modl, Options *optn, int nini, int *intlPrt, in { // Define a process and construct a set of sprocesses // model : Model object - // opt : Option object + // optn : Option object // initlPart : list of particle-id of initlPart particles // finalPart : list of particle-id of final particles // coupling : list of coupling constants int j, lp2; Bool ok; + char buff[MAXSTR]; id = pid; model = modl; @@ -3121,10 +3401,10 @@ Process::Process(int pid, Model *modl, Options *optn, int nini, int *intlPrt, in erEnd("Process: illegal input"); } -#ifdef DEBUG - // print message - prProcess(); -#endif + if (opt->out->outgrp != NULL) { + snprintf(buff, MAXSTR, "out%d.prp", pid); + prProcessP(buff); + } // construct sprocesses mkSProcess(); @@ -3135,7 +3415,10 @@ Process::~Process(void) { initlPart = delintdup(initlPart); finalPart = delintdup(finalPart); - delete sproc; + for (int j = 0; j < nSubproc; j++) { + delete sptbl[j]; + } + // delete sproc; } //-------------------------------------------------------------- @@ -3146,6 +3429,76 @@ void Process::prProcess(void) prIntArray(model->ncouple, clist, "\n"); } +//-------------------------------------------------------------- +void Process::prProcessP(const char *fname) +{ + FILE *fp; + + if ((fp = fopen(fname, "w")) == NULL) { + fprintf(GRCC_Stderr, "*** cannot open \"%s\"\n", fname); + return; + } + fprintf(fp, "Process = {\n"); + outProcP(fp); + fprintf(fp, "};\n"); + fclose(fp); +} + +//-------------------------------------------------------------- +void Process::outProcP(FILE *fp) +{ + int j; + + if (model == NULL) { + return; + } + fprintf(fp, " \"Model\": \"%s\",\n", model->name); + fprintf(fp, " \"Initial\": ["); + for (j = 0; j < ninitl; j++) { + if (j != 0) { + fprintf(fp, ", "); + } + fprintf(fp, "\"%s\"", model->particleName(initlPart[j])); + } + fprintf(fp, "],\n"); + fprintf(fp, " \"Final\": ["); + for (j = 0; j < nfinal; j++) { + if (j != 0) { + fprintf(fp, ", "); + } + fprintf(fp, "\"%s\"", model->particleName(finalPart[j])); + } + fprintf(fp, "],\n"); + fprintf(fp, " \"Couple\": ["); + for (j = 0; j < model->ncouple; j++) { + if (j != 0) { + fprintf(fp, ", "); + } + fprintf(fp, "%d", clist[j]); + } + fprintf(fp, "],\n"); + fprintf(fp, " \"Options\": {\n"); + for (j = 0; j < GRCC_OPT_Size; j++) { + if (j == GRCC_OPT_Outgrf || j == GRCC_OPT_Outgrp) { + continue; + } + fprintf(fp, " \"%s\":%d,\n", optDef[j].name, opt->values[j]); + } + fprintf(fp, " \"%s\":", optDef[GRCC_OPT_Outgrf].name); + if (opt->out->outgrf != NULL) { + fprintf(fp, "\"%s\",\n", opt->out->outgrf); + } else { + fprintf(fp, "0,\n"); + } + fprintf(fp, " \"%s\":", optDef[GRCC_OPT_Outgrp].name); + if (opt->out->outgrp != NULL) { + fprintf(fp, "\"%s\",\n", opt->out->outgrp); + } else { + fprintf(fp, "0,\n"); + } + fprintf(fp, " }\n"); +} + //-------------------------------------------------------------- void Process::mkSProcess(void) { @@ -3160,6 +3513,9 @@ void Process::mkSProcess(void) int nl[GRCC_MAXMINTERACT]; double proct0 = 0, proct1 = 0; + if (model->ncplgcp < 1) { + erEnd("function 'addInteractionEnd' has not been called"); + } // output file to start process opt->beginProc(this); @@ -3179,12 +3535,6 @@ void Process::mkSProcess(void) // for partitions of (leg, order) while (nextPart(ctotal, model->ncplgcp, model->cplgcp, nl, &r)) { -#ifdef DEBUG - printf("nextPart:ctotal=%d, nc=%d, c=", ctotal, model->ncplgcp); - prIntArray(model->ncplgcp, model->cplgcp, ", nl="); - prIntArray(model->ncplgcp, nl, ""); - printf(", r=%d\n", r); -#endif // count the total number of vertices and legs. nvtx = 0; nleg = 0; @@ -3299,13 +3649,6 @@ void Process::mkSProcess(void) nc++; } } -#ifdef DEBUG1 - for (j = 0; j < nc; j++) { - printf("%d/%d: deg=%d, typ=%d, ptcl=%d, cple=%d, num=%d\n", - j, nc, cls[j].cldeg, cls[j].cltyp, cls[j].ptcl, - cls[j].cple, cls[j].clnum); - } -#endif // create a sprocess ??? to be rewritten sproc = new SProcess(model, this, opt, nSubproc, clist, nc, cls); @@ -3315,11 +3658,6 @@ void Process::mkSProcess(void) } sptbl[nSubproc++] = sproc; -#ifdef DEBUG - // print sprocesses - sproc->prSProcess(); -#endif - // generate M-graphs sproc->generate(); @@ -3339,8 +3677,8 @@ void Process::mkSProcess(void) } else { ngraphs = nAGraphs; } - delete sproc; - sproc = NULL; + // delete sproc; + // sproc = NULL; } // ending time @@ -3617,11 +3955,8 @@ MNode::MNode(int vid, int vclss, NCInput *mgi) extloop = mgi->cltyp; // external node or not cmindeg = mgi->cmind; // min(deg of connectable vertex) cmaxdeg = mgi->cmaxd; // max(deg of connectable vertex) -#ifdef DEBUG3 - printf("MNode:id=%d, freelg=%d, cmind=%d, cmaxd=%d\n", - id, freelg, cmindeg, cmaxdeg); -#endif } + //-------------------------------------------------------------- MNode::MNode(int vid, int vdeg, int vextlp, int vclss, int cmin, int cmax) { @@ -3638,9 +3973,6 @@ MNode::MNode(int vid, int vdeg, int vextlp, int vclss, int cmin, int cmax) extloop = vextlp; // external node or not cmindeg = cmin; // min(deg of connectable vertex) cmaxdeg = cmax; // max(deg of connectable vertex) -#ifdef DEBUG3 - printf("MNode:id=%d, freelg=%d, cmind=%d, cmaxd=%d\n", id, freelg, cmindeg, cmaxdeg); -#endif } //=============================================================== @@ -3650,15 +3982,6 @@ MGraph::MGraph(int pid, int ncl, int *cldeg, int *clnum, int *cltyp, int *cmind, { int nn, ne, j, k; -#ifdef DEBUGM - printf("MGraph::MGraph(pid=%d, ncl=%d)\n", pid, ncl); - for (j = 0; j < ncl; j++){ - printf("%d: deg=%d, num=%d, typ=%d, mind=%d, maxd=%d\n", - j, cldeg[j], clnum[j], cltyp[j], cmind[j], cmaxd[j]); - } - opts->print(); -#endif - // initial conditions nClasses = ncl; clist = new int[ncl]; @@ -3708,13 +4031,7 @@ MGraph::MGraph(int pid, int ncl, int *cldeg, int *clnum, int *cltyp, int *cmind, egraph = new EGraph(nNodes, nEdges, maxdeg); nn = 0; nExtern = 0; -#ifdef DEBUG3 - printf("MGraph::MGraph: nClasses=%d\n", nClasses); -#endif for (j = 0; j < nClasses; j++) { -#ifdef DEBUG3 - printf("cmind[%d]=%d, cmaxd[%d]=%d\n", j, cmind[j], j, cmaxd[j]); -#endif for (k = 0; k < clist[j]; k++, nn++) { nodes[nn] = new MNode(nn, cldeg[j], cltyp[j], j, cmind[j], cmaxd[j]); egraph->setExtLoop(nn, cltyp[j]); @@ -3733,16 +4050,6 @@ MGraph::MGraph(int pid, int ncl, NCInput *mgi, Options *opts) { int nn, ne, j, k; -#ifdef DEBUGM - printf("MGraph::MGraph(pid=%d, ncl=%d)\n", pid, ncl); - for (j = 0; j < ncl; j++){ - printf("%d: deg=%d, num=%d, typ=%d, mind=%d, maxd=%d\n", - j, mgi[j].cldeg, mgi[j].clnum, mgi[j].cltyp, - mgi[j].cmind, mgi[j].cmaxd); - } - opts->print(); -#endif - // initial conditions nClasses = ncl; clist = new int[ncl]; @@ -3780,6 +4087,10 @@ MGraph::MGraph(int pid, int ncl, NCInput *mgi, Options *opts) } pId = pid; nNodes = nn; + if (nNodes < 0) { + printf("*** nNodes = %d\n", nNodes); + erEnd("MGraph::MGrap : nNodes < 0"); + } nodes = new MNode*[nNodes]; group = new SGroup(); #ifdef ORBITS @@ -3792,14 +4103,7 @@ MGraph::MGraph(int pid, int ncl, NCInput *mgi, Options *opts) egraph = new EGraph(nNodes, nEdges, maxdeg); nn = 0; nExtern = 0; -#ifdef DEBUG3 - printf("MGraph::MGraph: nClasses=%d\n", nClasses); -#endif for (j = 0; j < nClasses; j++) { -#ifdef DEBUG3 - printf("cmind[%d]=%d, cmaxd[%d]=%d\n", - j, mgi[j].cmind, j, mgi[j].cmaxd); -#endif for (k = 0; k < clist[j]; k++, nn++) { nodes[nn] = new MNode(nn, j, mgi+j); egraph->setExtLoop(nn, mgi[j].cltyp); @@ -3850,9 +4154,10 @@ void MGraph::init(void) // work space for isIsomorphic modmat = newMat(nNodes, nNodes, 0); - // work space for bisearchM + // work space for bisearchME bidef = newArray(nNodes, 0); bilow = newArray(nNodes, 0); + bicol = newArray(nNodes, 0); bicount = 0; } @@ -3863,6 +4168,7 @@ MGraph::~MGraph(void) // group->delGroup(); + bicol = deleteArray(bicol); bilow = deleteArray(bilow); bidef = deleteArray(bidef); @@ -3928,6 +4234,37 @@ void MGraph::print(void) printAdjMat(curcl); } +//--------------------------------------------------------------- +void MGraph::printPy(FILE *fp, long mId) +{ + fprintf(fp, "#TGraph : (gseq, gid, nodes, amat)\n"); + fprintf(fp, "(%ld, (%d, %d),\n", mId, pId, -1); + fprintf(fp, " [ # nodes: (cid, deg0, loop, part)\n"); + for (int j = 0; j < nNodes; j++) { + fprintf(fp, " (%2d,%2d,%2d,%2d), # %2d\n", + nodes[j]->clss, nodes[j]->deg, nodes[j]->extloop, 0, j); + } + fprintf(fp, " ],\n"); + fprintf(fp, " [\n"); + fprintf(fp, " #"); + for (int j2 = 0; j2 < nNodes; j2++) { + fprintf(fp, " %4d", j2); + } + fprintf(fp, "\n"); + for (int j1 = 0; j1 < nNodes; j1++) { + fprintf(fp, " ["); + for (int j2 = 0; j2 < nNodes; j2++) { + fprintf(fp, " %2d, ", adjMat[j1][j2]); + } + fprintf(fp, "], # %2d\n", j1); + } + fprintf(fp, " ], # nodes\n"); + fprintf(fp, " [ # group of 2 elements\n"); + fprintf(fp, " [],\n"); + fprintf(fp, " ], # group\n"); + fprintf(fp, ")\n"); +} + //--------------------------------------------------------------- Bool MGraph::isConnected(void) { @@ -4154,19 +4491,22 @@ void MGraph::biconnME(void) { // Count the number of 1PI components. - int j, root, vr, next, nart; + int j, j1, root, vr, next, nart; MCOpi mopi; MCBlock mblk; + ULong momset; // initialization bicount = 0; - mconn->init(); + mconn->initCEdges(this); for (j = 0; j < nNodes; j++) { bidef[j] = -1; bilow[j] = -1; + bicol[j] = 0; } + bipart = True; // find a root for the root of bisearch // root must be an external node @@ -4189,7 +4529,7 @@ void MGraph::biconnME(void) } if (root >= 0) { - bisearchME(root, -1, 0, &mopi, &mblk, &next, &nart); + bisearchME(root, -1, 0, 1, &mopi, &mblk, &momset, &next, &nart); mconn->nlpopic = 0; mconn->nctopic = 0; @@ -4217,6 +4557,15 @@ void MGraph::biconnME(void) mconn->nselfloops += adjMat[j][j]/2; } + mconn->nmultiedges = 0; + for (j = 0; j < nNodes; j++) { + for (j1 = j+1; j1 < nNodes; j1++) { + if (adjMat[j][j1] > 1) { + mconn->nmultiedges++; + } + } + } + mconn->na1blocks = 0; for (j = 0; j < mconn->nblocks; j++) { if (mconn->blocks[j].nartps == 1) { @@ -4225,6 +4574,33 @@ void MGraph::biconnME(void) } mconn->neblocks = mconn->nblocks + mconn->nctopic; +#ifdef CHECK + if (True) { + if (isExternal(root)) { + momset |= MASK(root); + } + ULong m = 0; + for (j = 0; j < nExtern; j++) { + m |= MASK(j); + } + bool ok = True; + if (momset != m) { + printf("*** momset = %ld != %ld\n", momset, m); + ok = False; + } + if (mconn->nbacked != nLoops) { + printf("*** nbacked = %d != %c\n", mconn->nbacked, nLoops); + ok = False; + } + if (! ok) { + print(); + egraph->print(); + mconn->print(); + erEnd("biconnME: illegal connection"); + } + } +#endif + // no vertex. } else { // no vertices @@ -4236,7 +4612,8 @@ void MGraph::biconnME(void) } //--------------------------------------------------------------- -void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int *next, int *nart) +void MGraph::bisearchME(int nd, int pd, int ned, int col, + MCOpi *mopi, MCBlock *mblk, ULong *momset, int *next, int *nart) { // Search biconnected component // visit : pd --> nd --> td @@ -4266,23 +4643,30 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int // Bool newv; int td, td0, lp, conn, next1, nart1, nart2, nvisit, ndart; - int opit, opit0, blkt; + int opit, opit0, blkt, l; MCOpi mopi1; MCBlock mblk1; + ULong momset1, momset2, momset3; mopi->init(); mblk->init(); - *next = 0; // # external below 'nd' inclusive - *nart = 0; // # articulation points 'nd' inclusive - nart1 = 0; // # articulation points below 'nd'. - nart2 = 0; // # 'nd' is articulation point then 1 - ndart = 0; // # the number of blocks attached to 'nd' + *momset = 0; // # the set of momenta + *next = 0; // # external below 'nd' inclusive + *nart = 0; // # articulation points 'nd' inclusive + nart1 = 0; // # articlulation points below 'nd'. + nart2 = 0; // # 'nd' is articlulation point then 1 + ndart = 0; // # the number of blocks attached to 'nd' newv = (bidef[nd] < 0); + if (! newv) { + printf("*** nd=%d, pd=%d, bidef[nd]=%d\n", nd, pd, bidef[nd]); + erEnd("bisearchME: illegal connection"); + } bidef[nd] = bicount; bilow[nd] = bicount; + bicol[nd] = col; bicount++; opit0 = mconn->opistkptr; @@ -4294,11 +4678,15 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int // external node and not the root if (isExternal(nd)) { if (pd >= 0) { + // not the root + *momset = MASK(nd); mopi->next = 1; mopi->nlegs = 1; + mopi->mom0lg = 0; mblk->nartps = 1; *next = 1; *nart = 1; + mconn->addCEdge(nd, pd, *momset); return; } } @@ -4314,16 +4702,22 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int } nvisit++; td0 = td; + + momset1 = 0; // external node if (isExternal(td) && bidef[td] < 0) { - mopi->nlegs++; - mopi->next++; + bicol[td] = - col; (*next)++; ndart++; mconn->addArtic(nd, 1); mblk->nartps++; nart2 = 1; + momset1 = MASK(td); + mopi->nlegs++; + mopi->next++; + mconn->addCEdge(td, nd, momset1); + *momset |= momset1; if (newv) { mconn->pushNode(td); } @@ -4334,14 +4728,19 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int mopi->loop += lp; mopi->nedges += lp; ndart += lp; + bipart = False; mconn->addArtic(nd, lp); mconn->addBlockSelf(nd, lp); mblk->nartps++; nart2 = 1; - + for (l = 0; l < lp; l++) { + int m = nExtern + mconn->nbacked; + mconn->nbacked++; + momset1 = MASK(m); + mconn->addCEdge(td, nd, momset1); + } // back to the parent : pd --> nd --> pd, pd is a vertex - // Back edges when the connection is multi-edges (ned > 1) } else if (td == pd) { if (ned > 1) { bilow[nd] = Min(bilow[nd], bidef[pd]); @@ -4352,14 +4751,56 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int nvisit--; } - // back edge :td is already visited + // td is already visited } else if (bidef[td] >= 0) { + bipart &= (bicol[td] == - col); bilow[nd] = Min(bilow[nd], bidef[td]); if (bidef[nd] >= bidef[td]) { + // new back-edge is found mopi->loop += conn; mopi->nedges += conn; mblk->loop += conn; mconn->pushEdge(td, nd); + for (l = 0; l < conn; l++) { + int m = nExtern + mconn->nbacked; + mconn->nbacked++; + momset1 = MASK(m); + *momset |= momset1; + mconn->addCEdge(td, nd, momset1); + } + } else { + // reverse direction of back-edge + int cn = 0; + + for (int ed = 0; ed < mconn->sedges; ed++) { + if (mconn->cedges[ed].nodes[0] == nd && + mconn->cedges[ed].nodes[1] == td && + mconn->cedges[ed].momdir > 0) { + cn++; + *momset &= ~ (mconn->cedges[ed].momset); + } else if (mconn->cedges[ed].nodes[1] == nd && + mconn->cedges[ed].nodes[0] == td && + mconn->cedges[ed].momdir < 0) { + cn++; + *momset &= ~ (mconn->cedges[ed].momset); + } + } + if (cn != conn) { + printf("*** revisit back-ed (%d --> %d) cn=%d != conn=%d\n", + nd, td, cn, conn); + printf(" sedges = %d\n", mconn->sedges); + for (int ed = 0; ed < mconn->sedges; ed++) { + printf(" %d : (%d --> %d) : [%2d*%2ld]\n", ed, + mconn->cedges[ed].nodes[0], + mconn->cedges[ed].nodes[1], + mconn->cedges[ed].momdir, + mconn->cedges[ed].momset); + } + print(); + egraph->print(); + mconn->print(); + erEnd("bisearchME: illegal connection"); + } } // new node @@ -4377,17 +4818,33 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int mconn->pushEdge(nd, td); // visit child - bisearchME(td, nd, adjMat[td][nd], &mopi1, &mblk1, &next1, &nart1); + bisearchME(td, nd, adjMat[td][nd], - col, + &mopi1, &mblk1, &momset1, &next1, &nart1); + // momset + momset2 = 0; + for (l = 0; l < conn - 1; l++) { + int m = nExtern + mconn->nbacked; + mconn->nbacked++; + momset3 = MASK(m); + mconn->addCEdge(nd, td, momset3); + momset2 |= momset3; + } + mconn->addCEdge(td, nd, momset1 | momset2); + *momset |= momset1; // articulation point of bridge if (bilow[td] > bidef[nd]) { // new OPI component (not including 'td') + int mom0lg = (momset1 == 0) ? 1 : 0; mopi1.nlegs++; + mopi1.mom0lg += mom0lg; mconn->addOPIc(&mopi1, opit); + mopi1.init(); - mopi1.nlegs = 1; + mopi1.nlegs = 1; + mopi1.mom0lg = mom0lg; if (pd >= 0 || !isExternal(nd)) { // opi mconn->addBridge(nd, td, next1, nExtern); @@ -4429,14 +4886,14 @@ void MGraph::bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int mopi->loop += mopi1.loop; mopi->nedges += mopi1.nedges; mopi->ctloop += mopi1.ctloop; + mopi->mom0lg += mopi1.mom0lg; mblk->nartps += mblk1.nartps; mblk->loop += mblk1.loop; - *next += next1; - *nart += nart1; - } - + *next += next1; + *nart += nart1; + } // of "if (bidef[td] < 0)" } // end of for // no connection except for 'pd'==>'nd'. Then 'nd' is an art. point @@ -4483,10 +4940,6 @@ BigInt MGraph::generate(void) MNodeClass *cl; -#ifdef DEBUGM - printf("MGraph::generate:\n"); - // opt->printLevel(2); -#endif // Initial classification of nodes. cl = new MNodeClass(nNodes, nClasses); cl->init(clist, maxdeg, adjMat); @@ -4495,9 +4948,6 @@ BigInt MGraph::generate(void) delete cl; -#ifdef DEBUGM - printf("MGraph::generate:end:%ld\n", cDiag); -#endif return cDiag; } @@ -4511,10 +4961,6 @@ void MGraph::connectClass(MNodeClass *cl) xcl = refineClass(cl); -#ifdef DEBUG3 - printf("connectClass\n"); - print(); -#endif if (xcl == NULL) { #ifdef MONITOR discardRefine++; @@ -4544,10 +4990,6 @@ void MGraph::connectNode(int so, int ss, MNodeClass *cl) return; } -#ifdef DEBUG3 - \printf("connectNode\n"); - print(); -#endif for (sn = ss; sn < cl->flist[sc+1]; sn++) { connectLeg(so, sn, so, sn, cl); return; @@ -4571,10 +5013,6 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) return; } #endif -#ifdef DEBUG3 - printf("connectLeg:0\n"); - print(); -#endif // There remains no free legs in the node 'sn' : move to next node. if (nodes[sn]->freelg < 1) { @@ -4584,19 +5022,12 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) discardDisc++; #endif } else { -#ifdef DEBUG3 - printf("call connectNode:1\n"); -#endif // next node in the current class. connectNode(so, sn+1, cl); } return; } -#ifdef DEBUG3 - printf("connectLeg:2\n"); - print(); -#endif // connect a free leg of the current node 'sn'. for (to1 = to; to1 < cl->nClasses; to1++) { tc = cl->clord[to1]; @@ -4609,27 +5040,12 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) if (ts1 >= nNodes) { continue; } -# ifdef DEBUG3 - printf("connectLeg:0: sn=%d, ts1=%d, ?(%d <= %d <= %d)\n", - sn, ts1, nodes[sn]->cmindeg, nodes[ts1]->deg, nodes[sn]->cmaxdeg); - printf("connectLeg:0: ts1=%d, sn=%d, ?(%d <= %d <= %d)\n", - ts1, sn, nodes[ts1]->cmindeg, nodes[sn]->deg, nodes[ts1]->cmaxdeg); -# endif if ((nodes[sn]->cmindeg > 0 && nodes[ts1]->deg < nodes[sn]->cmindeg) ||(nodes[sn]->cmaxdeg > 0 && nodes[ts1]->deg > nodes[sn]->cmaxdeg) | (nodes[ts1]->cmindeg > 0 && nodes[sn]->deg < nodes[ts1]->cmindeg) ||(nodes[ts1]->cmaxdeg > 0 && nodes[sn]->deg > nodes[ts1]->cmaxdeg)) { -# ifdef DEBUG3 - printf("connectLeg:1: sn=%d, ts1=%d, !(%d <= %d <= %d)\n", - sn, ts1, nodes[sn]->cmindeg, nodes[ts1]->deg, nodes[sn]->cmaxdeg); - printf("connectLeg:2: ts1=%d, sn=%d, !(%d <= %d <= %d)\n", - ts1, sn, nodes[ts1]->cmindeg, nodes[sn]->deg, nodes[ts1]->cmaxdeg); -# endif continue; } -#ifdef DEBUG3 - printf("connectLeg:3: sn=%d, ts1=%d\n", sn, ts1); -#endif #endif for (tn = ts1; tn < cl->flist[tc+1]; tn++) { if (sc == tc && sn > tn) { @@ -4652,9 +5068,11 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) maxself = nodes[sn]->freelg/2; } - // If we can assume no tadpole, the following line can be used. - if ((opt->values[GRCC_OPT_1PI] > 0 - || opt->values[GRCC_OPT_NoTadpole] > 0) && nExtern > 2) { + if ((nExtern > 1) && (opt->values[GRCC_OPT_1PI] > 0 + || opt->values[GRCC_OPT_NoTadpole] > 0) + && nNodes > 1) { + // there are two or more nodes in the graph : + // avoid to generate tadpole maxself = Min((nodes[sn]->deg-2)/2, maxself); } @@ -4667,9 +5085,6 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) cl->incMat(sn, tn, ncm); // next connection -#ifdef DEBUG3 - printf("call connectLeg:1: %d--%d\n", sn, sn); -#endif connectLeg(so, sn, to1, tn+1, cl); // restore the configuration @@ -4687,6 +5102,10 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) if (nNodes > 2 && nodes[sn]->deg == nodes[tn]->deg) { maxcon = Min(maxcon, nodes[sn]->deg-1); } + + if (opt->values[GRCC_OPT_NoMultiEdge] > 0) { + maxcon = Min(maxcon, 1); + } #ifdef CHECK if ((adjMat[sn][tn] != 0) || @@ -4709,9 +5128,6 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) cl->incMat(tn, sn, ncm); // next connection -#ifdef DEBUG3 - printf("call connectLeg:2: %d--%d\n", sn, tn); -#endif connectLeg(so, sn, to1, tn+1, cl); // restore configuration @@ -4732,82 +5148,64 @@ Bool MGraph::isOptM(void) { Bool ok = True; - opi = (mconn->nopic == 1); - opiloop = (mconn->nlpopic <= 1); - tadpole = (mconn->ne0bridges >= ((nExtern == 1) ? 0 : 1)); - selfloop = (mconn->nselfloops > 0); - tadblock = (mconn->na1blocks > 0); - block = (mconn->neblocks == 1); - extself = (mconn->ne1bridges > 0); + opi = (mconn->nopic == 1); + opiloop = (mconn->nlpopic <= 1); + tadpole = (mconn->ne0bridges >= ((nExtern == 1) ? 0 : 1)); + selfloop = (mconn->nselfloops > 0); + multiedge = (mconn->nmultiedges > 0); + tadblock = (mconn->na1blocks > 0); + block = (mconn->neblocks == 1); + extself = (mconn->ne1bridges > 0); if (opt->values[GRCC_OPT_1PI] > 0) { ok = ok && opi; -#ifdef DEBUGM - printf("isOptM: 1PI:%d nopic=%d\n", ok, mconn->nopic); -#endif } else if (opt->values[GRCC_OPT_1PI] < 0) { ok = ok && !opi; -#ifdef DEBUGM - printf("isOptM: -1PI:%d nopic=%d\n", ok, mconn->nopic); -#endif } if (opt->values[GRCC_OPT_NoExtSelf] > 0) { ok = ok && !extself; -#ifdef DEBUGM - printf("isOptM: NoExtSelf:%d, extself=%d\n", ok, extself); -#endif } else if (opt->values[GRCC_OPT_NoExtSelf] < 0) { ok = ok && extself; -#ifdef DEBUGM - printf("isOptM:-NoExtSelf:%d, extself=%d\n", ok, extself); -#endif } if (opt->values[GRCC_OPT_NoTadpole] > 0) { +#ifdef OLDOPT ok = ok && !tadpole; -#ifdef DEBUGM - printf("isOptM: NoTadPole:%d, tadpole=%d\n", ok, tadpole); +#else + if (nExtern > 1) { + ok = ok && !tadpole; + } #endif } else if (opt->values[GRCC_OPT_NoTadpole] < 0) { - ok = ok && tadpole; -#ifdef DEBUGM - printf("isOptM: -NoTadPole:%d\n", ok); +#ifdef OLDOPT + ok = ok && !tadpole; +#else + if (nExtern > 1) { + ok = ok && tadpole; + } #endif } if (opt->values[GRCC_OPT_NoSelfLoop] > 0) { -#ifdef DEBUGM - printf("isOptM: NoSelfLoop:%d\n", ok); -#endif } else if (opt->values[GRCC_OPT_NoSelfLoop] < 0) { ok = ok && selfloop; -#ifdef DEBUGM - printf("isOptM: -NoSelfLoop:%d\n", ok); -#endif + } + if (opt->values[GRCC_OPT_NoMultiEdge] > 0) { + } else if (opt->values[GRCC_OPT_NoMultiEdge] < 0) { + ok = ok && multiedge; } if (opt->values[GRCC_OPT_No1PtBlock] > 0) { - ok = ok && !tadblock; -#ifdef DEBUGM - printf("isOptM: No1PtBlock:%d\n", ok); -#endif + if (nExtern > 1) { + ok = ok && !tadblock; + } } else if (opt->values[GRCC_OPT_No1PtBlock] < 0) { - ok = ok && tadblock; -#ifdef DEBUGM - printf("isOptM: -NoSelfLoop:%d\n", ok); -#endif + if (nExtern > 1) { + ok = ok && tadblock; + } } if (opt->values[GRCC_OPT_Block] > 0) { ok = ok && block; -#ifdef DEBUGM - printf("isOptM: Block:%d\n", ok); -#endif } else if (opt->values[GRCC_OPT_Block] < 0) { ok = ok && !block; -#ifdef DEBUGM - printf("isOptM: -Block:%d\n", ok); -#endif } -#ifdef DEBUGM - printf("isOptM:%d\n", ok); -#endif return ok; } @@ -4874,16 +5272,7 @@ void MGraph::newGraph(MNodeClass *cl) } else { egraph->mId = cDiag; } -#ifdef DEBUG1 - printf("call outmg\n"); - egraph->model->prModel(); - opt->print(); - egraph->print(); -#endif ok = (*(opt->outmg))(egraph, opt->argmg); -#ifdef DEBUG - printf("ok=%d\n", ok); -#endif } if (ok) { @@ -4911,44 +5300,22 @@ void MGraph::newGraph(MNodeClass *cl) #ifdef MONITOR printf("\n"); - printf("Graph : %ld (%ld) 1PIComp=%d 1PILoop=%d", - cDiag, ngen, n1PIComps, n1PILoop); + printf("Graph : %ld (%ld)", cDiag, ngen); printf(" sym. factor = (%ld*%ld)\n", nsym, esym); printAdjMat(cl); // cl->printMat(); # ifdef ORBITS orbits->print(); # endif -# ifdef DEBUGM - printf("refine: %ld\n", - nCallRefine); - printf("discarded for refinement: %ld\n", - discardRefine); - printf("discarded for disconnected: %ld\n", - discardDisc); - printf("discarded for duplication: %ld\n", - discardIso); -# endif #endif // go to next step -#ifdef DEBUGM - printf("newGraph:%ld:accepted\n", ngen); -#endif opt->newMGraph(this); } else { -#ifdef DEBUGM - printf("deleted by setOutMG-function\n"); -#endif } } else { -#ifdef DEBUGM - printf("newGraph:%ld:discardOptE\n", ngen); -#endif } } else { -#ifdef DEBUGM - printf("newGraph:%ld:discardOptM\n", ngen); -#endif + ; } } } @@ -5061,6 +5428,20 @@ void MOrbits::toOrbits(void) //************************************************************** // n-edge connected components +//============================================================ +// class MCEdge +//------------------------------------------------------------ +MCEdge::MCEdge(void) +{ + ; +} + +//------------------------------------------------------------ +MCEdge::~MCEdge(void) +{ + ; +} + //============================================================ // class MCOpi: one n-edge connected component //------------------------------------------------------------ @@ -5078,14 +5459,14 @@ MCOpi::~MCOpi(void) //------------------------------------------------------------ void MCOpi::init(void) { - nlegs = 0; - nedges = 0; + nodes = NULL; nnodes = 0; + nlegs = 0; next = 0; + nedges = 0; loop = 0; ctloop = 0; - nnodes = 0; - nodes = NULL; + mom0lg = 0; } //============================================================ @@ -5127,9 +5508,13 @@ void MCBlock::init(void) //------------------------------------------------------------ MConn::MConn(int nnod, int nedg) { + // nnod : the number of nodes + // nedg : the number of edges + snodes = nnod; sedges = nedg; + cedges = new MCEdge[sedges]; opics = new MCOpi[snodes]; bridges = new MCBridge[sedges]; blocks = new MCBlock[sedges]; @@ -5155,6 +5540,7 @@ MConn::~MConn(void) delete[] blocks; delete[] bridges; delete[] opics; + delete[] cedges; } //------------------------------------------------------------ @@ -5162,52 +5548,106 @@ void MConn::init(void) { int j; - nopic = 0; - nlpopic = 0; - nctopic = 0; - nbridges = 0; - ne0bridges = 0; - ne1bridges = 0; - nselfloops = 0; - - nblocks = 0; - na1blocks = 0; - narticuls = 0; - neblocks = 0; - - opistkptr = 0; - nopisp = 0; - blkstkptr = 0; - nblksp = 0; + nopic = 0; + nlpopic = 0; + nbacked = 0; + nctopic = 0; + nbridges = 0; + ne0bridges = 0; + ne1bridges = 0; + nselfloops = 0; + nmultiedges = 0; + nbacked = 0; + + nblocks = 0; + na1blocks = 0; + narticuls = 0; + neblocks = 0; + + opistkptr = 0; + nopisp = 0; + blkstkptr = 0; + nblksp = 0; for (j = 0; j < snodes; j++) { articuls[j] = 0; } -} - -//------------------------------------------------------------ -void MConn::pushNode(int nd) -{ - if (opistkptr >= snodes) { - erEnd("MConn::pushNode: opistkptr >= snodes"); +} + +//------------------------------------------------------------ +void MConn::pushNode(int nd) +{ + if (opistkptr >= snodes) { + erEnd("MConn::pushNode: opistkptr >= snodes"); + } + opistk[opistkptr++] = nd; +} + +//------------------------------------------------------------ +void MConn::pushEdge(int n0, int n1) +{ + if (blkstkptr >= sedges) { + erEnd("MConn::pushEdge: blkstkptr >= sedges"); + } + if (n0 <= n1) { + blkstk[blkstkptr][0] = n0; + blkstk[blkstkptr][1] = n1; + } else { + blkstk[blkstkptr][0] = n1; + blkstk[blkstkptr][1] = n0; + } + blkstkptr++; +} + + +//------------------------------------------------------------ +void MConn::initCEdges(MGraph *mg) +{ + int ed, n0, n1, e; + + ed = 0; + for (n0 = 0; n0 < mg->nNodes; n0++) { + for (e = 0; e < mg->adjMat[n0][n0]/2; e++, ed++) { + cedges[ed].nodes[0] = n0; + cedges[ed].nodes[1] = n0; + cedges[ed].momdir = 0; + } + for (n1 = n0+1; n1 < mg->nNodes; n1++) { + for (e = 0; e < mg->adjMat[n0][n1]; e++, ed++) { + cedges[ed].nodes[0] = n0; + cedges[ed].nodes[1] = n1; + cedges[ed].momdir = 0; + } + } + } + if (ed != sedges) { + printf("*** ed=%d != sedges=%d\n", ed, sedges); + erEnd("MConn::initCEdge: table overflow"); } - opistk[opistkptr++] = nd; } //------------------------------------------------------------ -void MConn::pushEdge(int n0, int n1) +void MConn::addCEdge(int n0, int n1, ULong momset) { - if (blkstkptr >= sedges) { - erEnd("MConn::pushEdge: blkstkptr >= sedges"); - } + int m0, m1, dir, ed; + if (n0 <= n1) { - blkstk[blkstkptr][0] = n0; - blkstk[blkstkptr][1] = n1; + m0 = n0; + m1 = n1; + dir = 1; } else { - blkstk[blkstkptr][0] = n1; - blkstk[blkstkptr][1] = n0; + m0 = n1; + m1 = n0; + dir = -1; + } + for (ed = 0; ed < sedges; ed++) { + if (cedges[ed].nodes[0] == m0 && cedges[ed].nodes[1] == m1 && + cedges[ed].momdir == 0) { + cedges[ed].momdir = dir; + cedges[ed].momset = momset; + return; + } } - blkstkptr++; } //------------------------------------------------------------ @@ -5229,6 +5669,7 @@ void MConn::addOPIc(MCOpi *mopi, int stp) opics[nopic].next = mopi->next; opics[nopic].loop = mopi->loop; opics[nopic].ctloop = mopi->ctloop; + opics[nopic].mom0lg = mopi->mom0lg; opics[nopic].nodes = opisp + nopisp; nopisp += nn; opistkptr = stp; @@ -5312,18 +5753,31 @@ void MConn::print(void) int j, k; printf("+++ MConn object: snodes=%d, sedges=%d\n", snodes, sedges); - printf(" nopic=%d, nlpopic=%d, nctopic=%d, nbridges=%d, " - "ne0bridges=%d, ne1bridges=%d\n", - nopic, nlpopic, nctopic, nbridges, ne0bridges, ne1bridges); + printf(" nopic=%d, nlpopic=%d, nbacked=%d, nctopic=%d, " + "nbridges=%d, ne0bridges=%d, ne1bridges=%d\n", + nopic, nlpopic, nbacked, nctopic, + nbridges, ne0bridges, ne1bridges); printf(" nblocks=%d, neblocks=%d, na1blocks=%d, narticuls=%d, nselfloop=%d\n", nblocks, neblocks, na1blocks, narticuls, nselfloops); + printf(" nmultiedges=%d\n", nmultiedges); + + printf(" cEdges (%d)\n", sedges); + if (sedges > 0) { + for (j = 0; j < sedges; j++) { + printf(" %2d: (%d,%d)[%2d*%2ld] \n", j, + cedges[j].nodes[0], cedges[j].nodes[1], + cedges[j].momdir, cedges[j].momset); + } + } + printf("\n"); printf(" 1PI components (%d)\n", nopic); for (j = 0; j < nopic; j++) { printf(" %d: nleg=%d, nnodes=%d, nedge=%d, ", j, opics[j].nlegs, opics[j].nnodes, opics[j].nedges); - printf("next=%d, loop=%d, ctlp=%d: [", - opics[j].next, opics[j].loop, opics[j].ctloop); + printf("next=%d, loop=%d, ctlp=%d: m0lg=%d [", + opics[j].next, opics[j].loop, opics[j].ctloop, + opics[j].mom0lg); for (k = 0; k < opics[j].nnodes; k++) { printf(" %d", opics[j].nodes[k]); } @@ -5334,7 +5788,7 @@ void MConn::print(void) if (nbridges > 0) { printf(" "); for (j = 0; j < nbridges; j++) { - printf("(%d,%d)[ex=%d] ", + printf("(%d,%d)[mom=%d] ", bridges[j].nodes[0], bridges[j].nodes[1], bridges[j].next); } printf("\n"); @@ -5366,6 +5820,25 @@ void MConn::print(void) } } +//------------------------------------------------------------ +void MConn::prEdges(void) +{ + int j; + + printf(" cEdges (%d)", sedges); + if (sedges > 0) { + for (j = 0; j < sedges; j++) { + if (j % 5 == 0) { + printf("\n "); + } + printf("(%d,%d)[%2d*%3ld] ", + cedges[j].nodes[0], cedges[j].nodes[1], + cedges[j].momdir, cedges[j].momset); + } + printf("\n"); + } +} + //************************************************************** // sgroup.cc //============================================================== @@ -5600,8 +6073,14 @@ void ENode::setId(EGraph *egrph, const int nid) //-------------------------------------------------------------- void ENode::setExtern(int typ, int pt) { - ndtype = typ; - intrct = pt; + intrct = pt; + if (typ == GRCC_AT_Initial || typ == GRCC_AT_Final) { + extloop = typ; + } else { + fprintf(GRCC_Stderr, "** ENode::setExternal:: illegal typ=%d\n", + typ); + erEnd("ENode::setExternal:: illegal typ"); + } } //-------------------------------------------------------------- @@ -5653,6 +6132,9 @@ EEdge::EEdge(void) emom = NULL; lmom = NULL; + + momset = 0; + momdir = 0; } //-------------------------------------------------------------- @@ -5742,7 +6224,10 @@ void EEdge::print(void) if (egraph == NULL || !egraph->assigned) { printf("[%d, %d]", nodes[0], nodes[1]); } else { - printf("[(%d,%d), (%d,%d)]", nodes[0], nlegs[0], nodes[1], nlegs[1]); + printf("[(%d,%d), (%d,%d))", nodes[0], nlegs[0], nodes[1], nlegs[1]); + } + if (momdir != 0) { + printf("[%2d*%2lu]", momdir, momset); } if (egraph != NULL && egraph->bicount > 0) { @@ -5769,6 +6254,12 @@ void EEdge::print(void) if (zero) { printf(" 0"); } + } else { + if (egraph == NULL) { + printf(" egraph=NULL"); + } else { + printf(" egraph->bicount=%d", egraph->bicount); + } } printf("\n"); } @@ -5862,11 +6353,10 @@ EGraph::EGraph(int nnodes, int nedges, int mxdeg) extMom = new int[sEdges+1]; fsign = 1; - nflines = -1; + nFlines = -1; for (j = 0; j < GRCC_MAXFLINES; j++) { flines[j] = NULL; } - } //-------------------------------------------------------------- @@ -5946,7 +6436,6 @@ void EGraph::copy(EGraph *eg) void EGraph::setExtLoop(int nd, int val) { // set the node 'nd' being an external node (-1) or a looped vertex - nodes[nd]->extloop = val; } @@ -5973,11 +6462,11 @@ void EGraph::fromDGraph(DGraph *dg) int n0, n1, e; if (dg->nnodes > GRCC_MAXNODES) { - fprintf(GRCC_Stderr, "*** too many nodes\n"); + fprintf(GRCC_Stderr, "*** too many nodes (GRCC_MAXNODES)\n"); exit(1); } if (dg->nedges > GRCC_MAXEDGES) { - fprintf(GRCC_Stderr, "*** too many edges\n"); + fprintf(GRCC_Stderr, "*** too many edges (GRCC_MAXEDGES)\n"); exit(1); } @@ -6028,7 +6517,7 @@ void EGraph::fromDGraph(DGraph *dg) esym = 1; extperm = 1; multp = 1; -// maxdeg = maxdeg; + maxdeg = maxdeg; for (n0 = 0; n0 < dg->nnodes; n0++) { nodes[n0]->deg = 0; @@ -6040,6 +6529,7 @@ void EGraph::fromDGraph(DGraph *dg) edges[e]->nodes[0] = n0; edges[e]->nodes[1] = n1; edges[e]->ext = (isExternal(n0) || isExternal(n1)); + edges[e]->momdir = 0; nodes[n0]->edges[nodes[n0]->deg++] = I2Vedge(e, -1); nodes[n1]->edges[nodes[n1]->deg++] = I2Vedge(e, +1); @@ -6055,7 +6545,7 @@ void EGraph::fromDGraph(DGraph *dg) tc = 0; for (n0 = 0; n0 < nNodes; n0++) { if (isExternal(n0)) { - setExtern(n0, 1, GRCC_ND_Initial); + setExtern(n0, 1, GRCC_AT_Initial); } else { tc += 2*nodes[n0]->extloop + nodes[n0]->deg - 2; } @@ -6077,7 +6567,7 @@ void EGraph::fromMGraph(MGraph *mg) // This function should be called after // EGraph(), setExtLoop() and endSetExtLoop(). - int n0, n1, ed, e; + int n0, n1, m0, m1, ed, e; int j, ni, nf; PNodeClass *pnc; int k; @@ -6092,7 +6582,9 @@ void EGraph::fromMGraph(MGraph *mg) } #endif mgraph = mg; - sproc = mgraph->opt->sproc; + opt = mgraph->opt; + econn = mgraph->mconn; + sproc = opt->sproc; if (sproc == NULL) { proc = NULL; model = NULL; @@ -6100,8 +6592,6 @@ void EGraph::fromMGraph(MGraph *mg) proc = sproc->proc; model = sproc->model; } - opt = mgraph->opt; - econn = mgraph->mconn; nNodes = mg->nNodes; nEdges = mg->nEdges; @@ -6182,25 +6672,26 @@ void EGraph::fromMGraph(MGraph *mg) if (proc != NULL) { ni = proc->ninitl; for (j = 0; j < ni; j++) { - setExtern(j, proc->initlPart[j], GRCC_ND_Initial); + setExtern(j, proc->initlPart[j], GRCC_AT_Initial); } nf = proc->nfinal; for (j = 0; j < nf; j++) { - setExtern(j+ni, proc->finalPart[j], GRCC_ND_Final); + setExtern(j+ni, proc->finalPart[j], GRCC_AT_Final); } nExtern = ni + nf; } else if (sproc != NULL) { pnc = sproc->pnclass; for (j = 0; j < sproc->pnclass->nclass; j++) { - if (pnc->type[j] == GRCC_AT_Initial || pnc->type[j] == GRCC_AT_External) { + if (pnc->type[j] == GRCC_AT_Initial + || pnc->type[j] == GRCC_AT_External) { for (k = pnc->cl2nd[j]; k < pnc->cl2nd[j+1]; k++) { - setExtern(j, k, GRCC_ND_Initial); + setExtern(j, k, GRCC_AT_Initial); ni++; } } else if (sproc->pnclass->type[j] == GRCC_AT_Final) { for (k = pnc->cl2nd[j]; k < pnc->cl2nd[j+1]; k++) { - setExtern(j, k, GRCC_ND_Final); + setExtern(j, k, GRCC_AT_Final); nf++; } } @@ -6215,6 +6706,39 @@ void EGraph::fromMGraph(MGraph *mg) } } + for (ed = 0; ed < nEdges; ed++) { + edges[ed]->momdir = 0; + } + for (e = 0; e < econn->sedges; e++) { + m0 = econn->cedges[e].nodes[0]; + m1 = econn->cedges[e].nodes[1]; + + bool found = False; + for (ed = 0; ed < nEdges; ed++) { + if (edges[ed]->momdir != 0) { + continue; + } + n0 = edges[ed]->nodes[0]; + n1 = edges[ed]->nodes[1]; + if (m0 == n0 && m1 == n1) { + edges[ed]->momdir = econn->cedges[e].momdir; + edges[ed]->momset = econn->cedges[e].momset; + found = True; + break; + } else if (m0 == n1 && m1 == n0) { + edges[ed]->momdir = - econn->cedges[e].momdir; + edges[ed]->momset = econn->cedges[e].momset; + found = True; + break; + } + } + + if (! found) { + printf("*** EGraph::fromMGraph:edge (%d->%d) is not found\n", + m0, m1); + } + } + bicount = -1; } @@ -6234,9 +6758,9 @@ ENode *EGraph::setExtern(int n0, int pt, int ndtyp) } // particle comes into the node - if (ndtyp == GRCC_ND_Initial) { + if (ndtyp == GRCC_AT_Initial) { npt = pt; - } else if (ndtyp == GRCC_ND_Final) { + } else if (ndtyp == GRCC_AT_Final) { npt = -pt; } else { npt = 0; @@ -6253,7 +6777,7 @@ ENode *EGraph::setExtern(int n0, int pt, int ndtyp) // particle flows on e from leg=0 to 1. if (model != NULL){ npt = model->normalParticle(npt); - ept = model->normalParticle(-npt); + ept = model->normalParticle(npt); } else { npt = 0; ept = 0; @@ -6326,18 +6850,74 @@ void EGraph::print(void) } printf("\n"); - if (nflines > 0) { + if (nFlines > 0) { prFLines(); } } + +//-------------------------------------------------------------- +void EGraph::printPy(FILE *fp, long mId) +{ + if (!assigned) { + mgraph->printPy(fp, mId); + return; + } + + fprintf(fp, "#AGraph : (gseq, {spid, gid, ...}, nodes, node-group)\n"); + fprintf(fp, "(%ld,\n", mId); + + // {'proc':0, 'subproc':0, 'tgraph':1, ... } + fprintf(fp, " {'proc':%d, 'subproc':%d, ", proc->id, sproc->id); + fprintf(fp, "'tgraph':%ld, ", mId); + fprintf(fp, "'OPI':%s, ", BOOLSTR(opt->values[GRCC_OPT_1PI])); + fprintf(fp, "'agraph':%ld, \n", aId); + fprintf(fp, " 'nnodes':%d, 'nnsym':%ld, 'nesym':%ld, },\n", + nNodes, nsym, esym); + // #nodes + fprintf(fp, " [ #nodes\n"); + for (int nd = 0; nd < nNodes; nd++) { + ENode *end = nodes[nd]; + + // #node 0: external=phi + // #node 4: phi3=0 + fprintf(fp, " [ #node %d: ", nd); + fprintf(fp, "\n"); + + // {'deg':1, 'cid':-3, ..., } + fprintf(fp, " {'deg':%d, 'cid':%d, 'extern':%s, ", + end->deg, 0, BOOLSTR(isExternal(nd))); + fprintf(fp, "'intr':%d, 'loop':%d, },\n", + end->intrct, end->extloop); + + // legs : (edge, leg, particle) + // [ (0, 1, 0), ... ], + fprintf(fp, " [ "); + for (int lg = 0; lg < end->deg; lg++) { + int ed = V2Iedge(end->edges[lg]); + int ptcl = edges[ed]->ptcl; + int edlg = (edges[ed]->nodes[0] == nd) ? 1 : 0; + int nnd = edges[ed]->nodes[edlg]; + int nlg = edges[ed]->nlegs[edlg]; + fprintf(fp, "(%d, %d, %d), ", nnd, ptcl, nlg); + } + fprintf(fp, "],\n"); + fprintf(fp, " ],\n"); // end #node + } + fprintf(fp, " ], #node\n"); + fprintf(fp, " [ # group of 1 elements\n"); + fprintf(fp, " [], #\n"); + fprintf(fp, " ], # group\n"); + fprintf(fp, ")\n"); +} + //-------------------------------------------------------------- void EGraph::prFLines(void) { int j; printf(" Fermion lines %d, sign=%d (mId=%ld, aId=%ld)\n", - nflines, fsign, mId, aId); - for (j = 0; j < nflines; j++) { + nFlines, fsign, mId, aId); + for (j = 0; j < nFlines; j++) { printf("%4d ", j); flines[j]->print("\n"); } @@ -6436,6 +7016,244 @@ int EGraph::groupLMom(int *grp, int *ed2gr) return ngrp; } +//-------------------------------------------------------------- +Bool EGraph::optQGrafM(Options *opt) +{ + int *qgopt = opt->qgopt; + int nopis[GRCC_MAXEDGES]; + ULong mext = 0; + mext = (~ mext) << nExtern; + +#ifdef PRINT + printf("optQGrafM:"); + print(); +#endif + +#ifdef PRINT + printf("optQGrafM: %8ld\n", mId); + econn->print(); +#endif + + // count the number of self-energy 1PI components. + int maxlegs = 0; + for (int j = 0; j < econn->nopic; j++) { + maxlegs = Max(maxlegs, econn->opics[j].nlegs); + } + if (maxlegs >= GRCC_MAXEDGES) { + printf("*** table overflow\n"); + abort(); + } + for (int k = 0; k < GRCC_MAXEDGES; k++) { + nopis[k] = 0; + } + for (int j = 0; j < econn->nopic; j++) { + nopis[econn->opics[j].nlegs]++; + } + + // + if (qgopt[GRCC_QGRAF_OPT_ONEPI] > 0) { + if (econn->nopic != 1) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_ONEPI] < 0) { + if (econn->nopic < 2) { + return False; + } + } + + if (qgopt[GRCC_QGRAF_OPT_ONSHELL] != 0) { + if (nExtern == 1) { + if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) { + if (econn->nopic != 1) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) { + if (econn->nopic == 1) { + return False; + } + } + } else { + if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) { + if (econn->ne1bridges > 0) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) { + if (econn->ne1bridges <= 0) { + return False; + } + } + } + } + + if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] != 0) { + if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] > 0) { + if (nExtern == 1) { + if (econn->nblocks != 1) { + return False; + } + } else { + if (mgraph->mconn->ne0bridges >= 1 || + mgraph->mconn->na1blocks >= 1) { + return False; + } + } + } else if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] < 0) { + if (nExtern == 1) { + if (econn->nblocks == 1) { + return False; + } + } else { + if (mgraph->mconn->ne0bridges < 1 && + mgraph->mconn->na1blocks < 1) { + return False; + } + } + } + } + + if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] != 0) { + if (nExtern == 1) { + if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) { + if (econn->nopic != 1) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) { + if (econn->nopic == 1) { + return False; + } + } + } else { + if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) { + if (mgraph->mconn->ne0bridges != 0) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) { + if (mgraph->mconn->ne0bridges == 0) { + return False; + } + } + } + } + + if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] != 0) { + bool ok = True; + if (nExtern != 2) { + if (nopis[2] > 1) { + ok = False; + } + } + for (int j = 0; j < econn->nopic; j++) { + if (econn->opics[j].nlegs >= 2 && + econn->opics[j].nlegs == econn->opics[j].mom0lg) { + ok = False; + } + } + + // loop momenta + for (int j = 0; j < econn->sedges - 1; j++) { + ULong momj = econn->cedges[j].momset; + if (momj == 0) { + continue; + } + int extj = 0; + if (econn->cedges[j].nodes[0] < nExtern) { + extj = 1; + } else if (econn->cedges[j].nodes[1] < nExtern) { + extj = 1; + } + for (int k = j+1; k < econn->sedges; k++) { + ULong momk = econn->cedges[k].momset; + if (momk == 0) { + continue; + } + int extk = 0; + if (econn->cedges[k].nodes[0] < nExtern) { + extk = 1; + } else if (econn->cedges[k].nodes[1] < nExtern) { + extk = 1; + } + if ( momj == momk) { + if (nExtern == 2) { + if (extj + extk != 2) { + ok = False; + } + } else { + ok = False; + } + } + } + } + if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] > 0) { + if (! ok) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] < 0) { + if (ok) { + return False; + } + } + } + + if (qgopt[GRCC_QGRAF_OPT_SIMPLE] > 0) { + if (mgraph->selfloop || mgraph->multiedge) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_SIMPLE] < 0) { + if (!mgraph->selfloop && !mgraph->multiedge) { + return False; + } + } + + if (qgopt[GRCC_QGRAF_OPT_BIPART] > 0) { + if (! mgraph->bipart) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_BIPART] < 0) { + if (mgraph->bipart) { + return False; + } + } + // GRCC_QGRAF_OPT_CYCLI + if (qgopt[GRCC_QGRAF_OPT_CYCLI] != 0) { + int nb = 0; + for (int k = 0; k < econn->nblocks; k++) { + if (econn->blocks[k].loop > 0) { + nb++; + } + } + if (qgopt[GRCC_QGRAF_OPT_CYCLI] > 0) { + if (nb > 1) { + return False; + } + } else if (qgopt[GRCC_QGRAF_OPT_CYCLI] < 0) { + if (nb <= 1) { + return False; + } + } + } + + return True; +} + +//-------------------------------------------------------------- +Bool EGraph::optQGrafA(Options *opt) +{ +#ifdef PRINT + printf("optQGrafA: %8ld\n", mId); + econn->print(); +#endif + if (opt->qgopt[GRCC_QGRAF_OPT_FLOOP] != 0) { + for (int fl=0; fl < nFlines; fl++) { + if (flines[fl]->ftype == FL_Closed) { + if (flines[fl]->nlist % 2 != 0) { + return False; + } + } + } + } + return True; +} + //-------------------------------------------------------------- Bool EGraph::isOptE(void) { @@ -6446,10 +7264,11 @@ Bool EGraph::isOptE(void) Bool ok; int minopi2p; - if (opt->values[GRCC_OPT_No2PtL1PI] == 0 - && opt->values[GRCC_OPT_NoAdj2PtV] == 0) { - return True; - } + + // if (opt->values[GRCC_OPT_No2PtL1PI] == 0 + // && opt->values[GRCC_OPT_NoAdj2PtV] == 0) { + // return True; + // } for (ed = 0; ed < nEdges; ed++) { edges[ed]->cut = False; @@ -6457,6 +7276,11 @@ Bool EGraph::isOptE(void) biconnE(); + // QGraf options + if (! optQGrafM(opt)) { + return False; + } + if (opt->values[GRCC_OPT_NoAdj2PtV] > 0) { if (nadj2ptv > 0) { return False; @@ -6696,11 +7520,6 @@ void EGraph::biconnE(void) #ifdef CHECK chkMomConsv(); #endif -#ifdef DEBUG - printf("biconnE:opiext=%d, opiloop=%d, opi2plp=%d, nopi2p=%d, nopicomp=%d, bconn=%d\n", - opiext, opiloop, opi2plp, nopi2p, nopicomp, bconn); -#endif - return; } @@ -7062,20 +7881,32 @@ int EGraph::isFermion(int ed) ptcl = edges[ed]->ptcl; ptype = model->particles[Abs(ptcl)]->ptype; - return (ptype == GRCC_PT_Dirac || ptype == GRCC_PT_Majorana || ptype == GRCC_PT_Ghost); + return (ptype == GRCC_PT_Dirac || ptype == GRCC_PT_Majorana + || ptype == GRCC_PT_Ghost); } //-------------------------------------------------------------- int EGraph::fltrace(int fk, int nd0, int *fl) { + // fk : kind of fermion + // = (GRCC_PT_Dirac, GRCC_PT_Majorana or GRCC_PT_Ghost) + // nd0 : the last node visited + // fl : list of signed edge on the fermion line. + // fl[j] : + // (V2Iedge(fl[j]), V2Ileg(fl[j])) is the next node + // fl[0] should be already defined + int nfl, k, i, nd, nl, e, ed, el, fgcnt, fkind, lk; nfl = 1; + // maximal possible length of a fline is nEdges for (k = 0; k < nEdges; k++) { if (k >= nfl) { printf("*** fltrace:illegal contorl: k=%d, nEdges=%d\n", k, nEdges); break; } + + // get next node : nd0 ---- nd (nl) e = fl[k]; ed = V2Iedge(e); el = V2Ileg(e); @@ -7095,6 +7926,8 @@ int EGraph::fltrace(int fk, int nd0, int *fl) fgcnt = 1; break; } + + // find leg of nd for going next fgcnt = 0; ed = 0; lk = 0; @@ -7102,6 +7935,8 @@ int EGraph::fltrace(int fk, int nd0, int *fl) e = nodes[nd]->edges[i]; ed = V2Iedge(e); fkind = model->particles[Abs(edges[ed]->ptcl)]->ptype; + // lk = 1 : (# of fkind particle) is even + // lk = -1 : (# of fkind particle) is odd if (fkind == fk) { if (lk == 1) { lk = -1; @@ -7118,18 +7953,21 @@ int EGraph::fltrace(int fk, int nd0, int *fl) if (fkind == fk) { // i should be the neighbor of nl // (nl, i) or (i, nl) should be a pair of fkind + // (nl, i) : lk = -1 + // (i, nl) : lk = 1 if ((lk == 1 && nl == i+1) || (lk == -1 && nl == i-1)) { edges[ed]->visited = True; fgcnt++; if (fgcnt == 1) { + // -e = (signed edge points the next node) fl[nfl++] = - e; break; } } } } - } - } + } // not visited + } // end of for k if (fgcnt == 0) { printf("*** fline: Fermion number is not conserved\n"); printf(" nd=%d, e=%d, ed=%d, fgcnt=%d\n", @@ -7141,9 +7979,9 @@ int EGraph::fltrace(int fk, int nd0, int *fl) nd, e, ed, fgcnt); } } - if (k >= nNodes || k >= nfl) { + if (k >= nEdges || k >= nfl) { printf("*** fline: illegal control\n"); - printf(" nfl=%d, ", nfl); + printf(" nEdges=%d, nfl=%d, k=%d, ", nEdges, nfl, k); prIntArray(nfl, fl, "\n"); erEnd("fline: illegal control"); } @@ -7157,7 +7995,7 @@ void EGraph::getFLines(void) int nextn, exto[GRCC_MAXNODES]; int e, ed, nd, floop, nswap, nfl, fkind, el, ptcl; - nflines = 0; + nFlines = 0; for (ed = 0; ed < nEdges; ed++) { edges[ed]->visited = False; } @@ -7187,8 +8025,12 @@ void EGraph::getFLines(void) continue; } - e = I2Vedge(ed, el); - fl[0] = - e; + // el ed el1 + // x----<----x e = (signed edge of (ed, el1)) + // nd ptcl(>= 0) + // external + + fl[0] = - I2Vedge(ed, el); nfl = 1; fkind = model->particles[Abs(edges[ed]->ptcl)]->ptype; edges[ed]->visited = True; @@ -7214,8 +8056,12 @@ void EGraph::getFLines(void) el = 1; } nd = edges[ed]->nodes[el]; - e = I2Vedge(ed, 1-el); - fl[0] = e; + + // el ed el1 + // x----<----x e = (signed edge of (ed, el1)) + // nd ptcl + + fl[0] = - I2Vedge(ed, el); nfl = 1; edges[ed]->visited = True; @@ -7241,21 +8087,22 @@ void EGraph::addFLine(const FLType ft, int fk, int nfl, int *fl) { int j; - if (nflines >= GRCC_MAXFLINES) { + if (nFlines >= GRCC_MAXFLINES) { erEnd("too many Fermion lines (GRCC_MAXEDGES)"); } - if (flines[nflines] == NULL) { - flines[nflines] = new EFLine(); + if (flines[nFlines] == NULL) { + flines[nFlines] = new EFLine(); } - flines[nflines]->ftype = ft; - flines[nflines]->fkind = fk; - flines[nflines]->nlist = nfl; + flines[nFlines]->ftype = ft; + flines[nFlines]->fkind = fk; + flines[nFlines]->nlist = nfl; for (j = 0; j < nfl; j++) { - flines[nflines]->elist[j] = fl[j]; + flines[nFlines]->elist[j] = fl[j]; } - nflines++; + nFlines++; } + //************************************************************** // assign.cc //============================================================== @@ -7264,22 +8111,6 @@ void EGraph::addFLine(const FLType ft, int fk, int nfl, int *fl) // method : selection of assignable node -#ifdef DEBUGM -static int nordleg = 0; -static int nopleg = 0; -static int niso = 0; -static int niso1 = 0; -static int niso11 = 0; -static int niso111 = 0; -static int niso112 = 0; -static int niso12 = 0; -static int niso13 = 0; -static int niso14 = 0; -static int niso2 = 0; -static int nivord = 0; -static int nextonly = 0; -#endif - //=============================================================== // class NCand NCand::NCand(const NCandSt sta, const int dega, const int nilst, int *ilst) @@ -7456,15 +8287,6 @@ Assign::Assign(SProcess *sprc, MGraph *mgr, PNodeClass *pnc) erEnd("Assign: astack == NULL"); } -#ifdef DEBUG - if (pnclass == NULL) { - printf("+++ Assign::Assign : pnclass = NULL\n"); - } else { - printf("+++ Assign::Assign : pnclass :\n"); - pnclass->prPNodeClass(); - } -#endif - nNodes = mgraph->nNodes; nEdges = mgraph->nEdges; nExtern = sproc->nExtern; @@ -7500,23 +8322,8 @@ Assign::Assign(SProcess *sprc, MGraph *mgr, PNodeClass *pnc) #endif if (ok) { - // for debugging -#ifdef DEBUG0 - printf("+++ End of initialization : candidate:\n"); - prCand("init"); - printf("\n"); -#endif - // start assignment assignAllVertices(); -#ifdef DEBUGM - printf("mId=%ld, sId=%ld, ordleg=%d, ivord=%d, extonly=%d ", - egraph->mId, egraph->sId, nordleg, nivord, nextonly); - printf("niso=%d [%d(%d {%d, %d}, %d, %d, %d) %d]\n", - niso, niso1, niso11, niso111, niso112, niso12, niso13, niso14, - niso2); -#endif - } else { // cannot assign } @@ -7619,10 +8426,6 @@ Bool Assign::assignAllVertices(void) checkCand("assignAllVertices:1"); #endif -#ifdef DEBUG0 - printf("+++ particle assignment for '%ld'\n", mgraph->mId); -#endif - // start main part #ifdef SIMPSEL ok = selectVertexSimp(-1); @@ -7630,16 +8433,6 @@ Bool Assign::assignAllVertices(void) ok = selectVertex(); #endif -#ifdef DEBUG0 - printf("\n"); - printf("+++ Total %ld assigned graphs for '%ld'\n", - nAGraphs, mgraph->mId); - printf("result: %ld ", nAGraphs); - wAGraphs.print(" "); - printf("%ld ", nAOPI); - wAOPI.print("\n"); -#endif - #ifdef CHECK checkCand("assignAllVertices:2"); #endif @@ -7670,6 +8463,10 @@ Bool Assign::selectVertexSimp(int lastv) ok = allAssigned(); if (ok) { + if (!egraph->optQGrafA(opt)) { + return False; + } + // egraph->biconnE(); // necessary ??? opt->newAGraph(egraph); } @@ -7705,6 +8502,10 @@ Bool Assign::selectVertex(void) ok = allAssigned(); if (ok) { + if (!egraph->optQGrafA(opt)) { + return False; + } + egraph->biconnE(); opt->newAGraph(egraph); } @@ -7971,9 +8772,6 @@ Bool Assign::allAssigned(void) // check duplication by violating ordering condition if (!isOrdLegs()) { -#ifdef DEBUGM - nordleg++; -#endif return False; } @@ -7986,9 +8784,6 @@ Bool Assign::allAssigned(void) ok = isIsomorphic(cl, &nsym, &esym, &nsym1); if (!ok || nsym < 1 || esym < 1) { -#ifdef DEBUGM - niso++; -#endif return False; } @@ -8003,16 +8798,6 @@ Bool Assign::allAssigned(void) wAOPI.add(1,nsym*esym); } -#ifdef DEBUG1 - for (int j = 0; j < model->ncouple; j++) { - if (cplleft[j] != 0) { - printf("nAgraphs=%ld: 0 != cplleft =", nAGraphs); - prIntArray(model->ncouple, cplleft, "\n"); - break; - } - } -#endif - #ifdef CHECK checkAG("allAssigned"); #endif @@ -8027,20 +8812,12 @@ Bool Assign::allAssigned(void) ext = egraph->edges[e]->ext; if (!ext) { if (model->particles[p]->extonly) { -#ifdef DEBUGM - nextonly++; -#endif return False; } } } #endif -#ifdef DEBUG0 - printf("Assigned graph = %ld, sym = (%ld, %ld) ", nAGraphs, nsym, esym); - prCand("allAssigned "); -#endif - return True; } @@ -8170,7 +8947,7 @@ Bool Assign::fromMGraph(void) erEnd("fromMGraph: ptcl=0"); } #endif - ok = assignPLeg(n, 0, - ptcl); + ok = assignPLeg(n, 0, -ptcl); if (!ok) { // impossible config return False; @@ -8397,14 +9174,7 @@ Bool Assign::fillEGraph(int aid, BigInt nsym, BigInt esym, BigInt nsym1) // analyse fermion line and determine Fermi statistical sign factor - if (!model->skipFLine) { - egraph->getFLines(); - } else { - if (prlevel > 0) { - fprintf(GRCC_Stderr, "+++ Sign factors related to " - "Dirac/Majorana/Ghost particles are not calculated.\n"); - } - } + egraph->getFLines(); #ifdef CHECK checkAG("fillEGraph:0"); @@ -8742,9 +9512,6 @@ Bool Assign::assignPLeg(int n, int ln, int pt) } if (!isOrdPLeg(n, ln, pt)) { -#ifdef DEBUGM - nopleg++; -#endif return False; } @@ -9129,9 +9896,6 @@ Bool Assign::isIsomorphic(MNodeClass *cl, BigInt *nsym, BigInt *esym, BigInt *ns cmp = cmpPermGraph(p, cl); if (cmp < 0) { // duplicated graph -#ifdef DEBUGM - niso1++; -#endif return False; } else if(cmp == 0) { // not duplicated (*nsym)++; @@ -9163,9 +9927,6 @@ Bool Assign::isIsomorphic(MNodeClass *cl, BigInt *nsym, BigInt *esym, BigInt *ns // calculate permutations of edges *esym = edgeSym(); if (*esym < 1) { -#ifdef DEBUGM - niso2++; -#endif return False; } @@ -9187,18 +9948,11 @@ int Assign::cmpPermGraph(int *p, MNodeClass *cl) if (p == NULL) { erEnd("Assign::cmpPermGraph: p==NULL"); } -#endif -#ifdef DEBUG1 - printf("cmpPermGraph:0: p="); - prIntArray(nNodes, p, "\n"); #endif for (n = 0; n < nNodes; n++) { if (!isATExternal(pnclass->type[pnclass->nd2cl[n]])) { cmp = cmpNodes(n, p[n], cl); if (cmp != 0) { -#ifdef DEBUGM - if (cmp < 0) { niso11++; } -#endif return cmp; } } @@ -9221,9 +9975,6 @@ int Assign::cmpPermGraph(int *p, MNodeClass *cl) } cmp = mgraph->adjMat[n1][n2] - mgraph->adjMat[p1][p2]; if (cmp != 0) { -#ifdef DEBUGM - if (cmp < 0) { niso12++; } -#endif return cmp; } @@ -9243,9 +9994,6 @@ int Assign::cmpPermGraph(int *p, MNodeClass *cl) cmp = njn - njp; if (cmp != 0) { -#ifdef DEBUGM - if (cmp < 0) { niso13++; } -#endif return cmp; } @@ -9258,9 +10006,6 @@ int Assign::cmpPermGraph(int *p, MNodeClass *cl) for (j = 0; j < njn; j++) { cmp = jn[j] - jp[j]; if (cmp != 0) { -#ifdef DEBUGM - if (cmp < 0) { niso14++; } -#endif return cmp; } } @@ -9282,17 +10027,11 @@ int Assign::cmpNodes(int nd0, int nd1, MNodeClass *cn) // Wether two nodes are in a same class or not. cmp = cn->ndcl[nd0] - cn->ndcl[nd1]; if (cmp != 0) { -#ifdef DEBUGM - if (cmp < 0) { niso111++; } -#endif return cmp; } // interaction cmp = nodes[nd0]->cand->ilist[0] - nodes[nd1]->cand->ilist[0]; -#ifdef DEBUGM - if (cmp < 0) { niso112++; } -#endif return cmp; } @@ -9807,14 +10546,25 @@ static void erEnd(const char *msg) if (erExit != NULL) { (*erExit)(msg, erExitArg); } - fprintf(GRCC_Stderr, "*** Error : %s\n", msg); + fprintf(GRCC_Stderr, "*** Error : %s\n\n", msg); GRCC_ABORT(); } +#define MAXPART 100 //------------------------------------------------------------ static Bool nextPart(int nelem, int nclist, int *clist, int *nl, int *r) { - int rem, pn, j; + // Generate configuration nl[] sequentially such that + // sum_j^{nclist} nl[j]*clist[j] = nelem + // 0 <= nl[j] (j = 0, ..., nclis-1) + // For control + // for the first call : *r < 0 + // otherwise : *r >= 0, nl is the last configuration + // Returns + // True : succeeded + // False : no more configuration + + int rem, pn, j, c; if (*r < 0) { *r = 0; @@ -9829,7 +10579,7 @@ static Bool nextPart(int nelem, int nclist, int *clist, int *nl, int *r) } else { rem = 0; } - for (int c = 0; c < 100; c++) { + for (c = 0; c < MAXPART; c++) { rem += nl[nclist-1]*clist[nclist-1]; for (pn = nclist-2; pn >= 0 && nl[pn] == 0; pn--) { ; @@ -9847,7 +10597,7 @@ static Bool nextPart(int nelem, int nclist, int *clist, int *nl, int *r) return True; } } - printf("*** nextPart : too many repetition\n"); + erEnd("*** nextPart : illegal control : too many repetition"); return False; } diff --git a/sources/grcc.h b/sources/grcc.h index 2a781a0c..4efb9d8e 100644 --- a/sources/grcc.h +++ b/sources/grcc.h @@ -5,6 +5,7 @@ #include #include +#define CHECK //============================================================== extern "C" { #include "grccparam.h" @@ -17,24 +18,12 @@ extern "C" { namespace Grcc { #endif -//============================================================== -// For debugging -// -// #define MONITOR // in graph elimination -// #define CHECK // check consistency -#define DEBUGABORT -// #define DEBUG -// #define DEBUG0 -// #define DEBUG1 -// #define DEBUGM - //============================================================== // Common types // +#define True 1 +#define False 0 typedef int Bool; -const int True = 1; -const int False = 0; - //============================================================== // Big integer (may be replaced by suitable one) @@ -63,6 +52,7 @@ class ENode; class Interaction; class MGraph; class MNodeClass; +class MCEdge; class MCOpi; class MCBridge; class MCBlock; @@ -100,9 +90,13 @@ class Options { void *argag; // additional argument to outag //switches for graph generation - int values[GRCC_OPT_Size]; // array of options + OptQGRef qgref[2*GRCC_QGRAF_OPT_Size]; // array of reference to QG-options + int values[GRCC_OPT_Size]; // array of options + int qgopt[GRCC_QGRAF_OPT_Size]; // array of QGRAF options - int DUMMYPADDING; + int nqgopt; // effective length of qgref + + int DUMMYPADDING; // measuring time double time0; @@ -111,11 +105,14 @@ class Options { //------------------ Options(void); ~Options(void); - void setDefaultValue(void); + void setDefaultValues(void); + void setOldDefaultValues(void); void setValue(int ind, int val); int getValue(int ind); + void setQGrafOpt(int *qgopt); + void print(void); void setOutputF(Bool outgrf, const char *fname); @@ -128,6 +125,8 @@ class Options { void setEndMG(OutEGB *omg, void *pt); void setErExit(ErExit *ere, void *pt); const OptDef *getDef(void); + const OptDef *getOldDef(void); + const OptQGDef *getQGDef(void); void begin(Model *mdl); void end(void); @@ -293,9 +292,6 @@ class Model { int ncplgcp; // the number of classes int defpart; // GRCC_DEFBYNAME or GRCC_DEFBYCODE - int skipFLine; // skip analysis fermion line - - int DUMMYPADDING; // methods Model(MInput *minp); @@ -317,6 +313,10 @@ class Model { int *allParticles(int *len); int findMClass(const int cpl, const int dgr); void prParticleArray(int n, int *a, const char *msg); + + static void printMInput(MInput *min); + static void printPInput(PInput *pin); + static void printIInput(IInput *iin); }; //************************************************************** @@ -439,7 +439,6 @@ class Process { int maxnlegs; // the maximum possible degree of node int clist[GRCC_MAXNCPLG]; // coupling constans of the process - // table of sprocesses int nSubproc; // the number of sprocesses SProcess *sptbl[GRCC_MAXSUBPROCS]; // the table of sprocesses @@ -472,6 +471,8 @@ class Process { ~Process(void); void prProcess(void); + void outProcP(FILE *fp); + void prProcessP(const char *fname); void mkSProcess(void); }; @@ -510,7 +511,7 @@ class ENode { int id; // id of the enode int maxdeg; // maximum degree int deg; // degree - int extloop; // loop inside this node + int extloop; // loop inside this node or AT_Initial etc. int ndtype; // type of the node int intrct; // assigned interaction/particle (ext.) @@ -552,6 +553,10 @@ class EEdge { int *lmom; // loop momenta int *extMom; // set of external momenta. + // momentum obtain in searchME + ULong momset; // set of momenta in bit string (leaf --> root) + int momdir; // direction (leg=0 --> leg=1) + Bool cut; int visited; int conid; // connected component @@ -559,6 +564,7 @@ class EEdge { int opicomp; // id of 1PI component int dir; // direction of momentum + int DUMMYPADDING; //-------------------------------- // functions @@ -648,7 +654,7 @@ class EGraph { // Fermion lines EFLine *flines[GRCC_MAXFLINES]; - int nflines; + int nFlines; int DUMMYPADDING1; @@ -659,8 +665,12 @@ class EGraph { void copy(EGraph *eg); void print(void); + void printPy(FILE *fp, long mId); void fromDGraph(DGraph *dg); void fromMGraph(MGraph *mgraph); + + Bool optQGrafM(Options *opt); + Bool optQGrafA(Options *opt); Bool isOptE(void); ENode *setExtern(int n0, int pt, int ndtp); @@ -691,7 +701,6 @@ class EGraph { void addFLine(const FLType ft, int fk, int nfl, int *fl); int legParticle(int ed, int lg); - }; //************************************************************** @@ -811,9 +820,11 @@ class MGraph { Bool opiloop; Bool extself; Bool selfloop; + Bool multiedge; Bool tadpole; Bool tadblock; Bool block; + Bool bipart; int DUMMYPADDING1; @@ -826,6 +837,7 @@ class MGraph { // work space for biconnected component int *bidef; int *bilow; + int *bicol; int bicount; int DUMMYPADDING2; @@ -841,6 +853,7 @@ class MGraph { Bool isExternal(int nd) { return (nodes[nd]->extloop < 0); }; void printAdjMat(MNodeClass *cl); void print(void); + void printPy(FILE *fp, long mId); Bool isConnected(void); Bool visit(int nd); @@ -848,15 +861,14 @@ class MGraph { void permMat(int size, int *perm, int **mat0, int **mat1); int compMat(int size, int **mat0, int **mat1); MNodeClass *refineClass(MNodeClass *cl); - void bisearchME(int nd, int pd, int ned, MCOpi *mopi, MCBlock *mblk, int *next, int *nart); + void bisearchME(int nd, int pd, int ned, int col, MCOpi *mopi, MCBlock *mblk, ULong *momset, int *next, int *nart); void biconnME(void); - Bool isOptM(void); + Bool isOptM(void); void connectClass(MNodeClass *cl); void connectNode(int sc, int ss, MNodeClass *cl); void connectLeg(int sc, int sn, int tc, int ts, MNodeClass *cl); void newGraph(MNodeClass *cl); - }; //=============================================================== @@ -895,6 +907,20 @@ class MNodeClass { void reorder(MGraph *mg); }; +//=============================================================== +// class of an edge +//-------------------------------------------------------------- +class MCEdge { + public: + Edge2n nodes; // nodes at the both size of the edge (leaf --> root) + ULong momset; // set of momenta in bit string + int momdir; // set of momenta in bit string + + int DUMMYPADDING; + + MCEdge(void); + ~MCEdge(void); +}; //=============================================================== // class of 1PI component @@ -903,11 +929,14 @@ class MCOpi { public: int *nodes; // array of nodes in int nnodes; // # nodes in the 1PI component - int nlegs; // # leg of the 1PI component + int nlegs; // # leg (bridges) of the 1PI component int next; // # external particles of the 1PI comp. int nedges; // # edges in the 1PI comp. int loop; // # loops in the 1PI comp. int ctloop; // # loops in the counter terms in the OP comp. + int mom0lg; // # leg (bridges) with 0 momentum + + int DUMMYPADDING; MCOpi(void); ~MCOpi(void); @@ -950,6 +979,7 @@ class MCBlock { class MConn { public: // 2-edge connected components + MCEdge *cedges; // table of edges MCOpi *opics; // table of n-edge-connected components MCBridge *bridges; // table of bridges MCBlock *blocks; // table of blocks @@ -966,10 +996,16 @@ class MConn { int nopic; // # 1PI components (n1PIComps) int nlpopic; // # looped 1PI components (n1PIComps) int nctopic; // # 1PI components of one counter term. + + // edges + int nbacked; // # back edges + + // bridges int nbridges; // # bridges int ne0bridges; // # bridges whose next=0 int ne1bridges; // # bridges whose next=1 int nselfloops; + int nmultiedges; // blocks (node-connected) int nblocks; // # blocks @@ -989,14 +1025,17 @@ class MConn { ~MConn(void); void init(void); + void initCEdges(MGraph *); void pushNode(int nd); void pushEdge(int n0, int n1); + void addCEdge(int n0, int n1, ULong momset); void addOPIc(MCOpi *mopi, int stp); void addBridge(int n0, int n1, int nex, int nextot); void addArtic(int nd, int mul); void addBlock(MCBlock *eblk, int stp); void addBlockSelf(int nd, int mul); void print(void); + void prEdges(void); }; //=============================================================== @@ -1042,7 +1081,7 @@ typedef enum { //=============================================================== class NCand { - // List of candidats for the assignment of interactions to a vertex + // List of candidates for the assignment of interactions to a vertex public: int deg; // degree of the node @@ -1060,7 +1099,7 @@ class NCand { //=============================================================== class ECand { - // List of candidats for the assignment of particles to an edge + // List of candidates for the assignment of particles to an edge public: Bool det; // determined or not diff --git a/sources/grccparam.h b/sources/grccparam.h index cbe96ced..869c1093 100644 --- a/sources/grccparam.h +++ b/sources/grccparam.h @@ -12,7 +12,7 @@ #define GRCC_MAXNCPLG 4 #define GRCC_MAXLEGS 10 #define GRCC_MAXMPARTICLES 50 -#define GRCC_MAXMINTERACT 200 +#define GRCC_MAXMINTERACT 500 #define GRCC_MAXSUBPROCS 500 #define GRCC_MAXNODES 20 #define GRCC_MAXEDGES 100 @@ -60,9 +60,6 @@ #define GRCC_PTS_Vector "vector" #define GRCC_PTS_Ghost "ghost" -#define GRCC_PTO_General 0 -#define GRCC_PTO_ExtOnly 1 - #define GRCC_PT_Table { GRCC_PTS_Undef, GRCC_PTS_Scalar, GRCC_PTS_Dirac, GRCC_PTS_Majorana, GRCC_PTS_Vector, GRCC_PTS_Ghost} #define GRCC_PT_GTable { "Undef", "Scalar", "Fermion", "Majorana", "Vector", "Ghost"} @@ -94,9 +91,12 @@ #define GRCC_OPT_NoExtSelf 8 #define GRCC_OPT_NoAdj2PtV 9 #define GRCC_OPT_Block 10 -#define GRCC_OPT_SymmInitial 11 -#define GRCC_OPT_SymmFinal 12 -#define GRCC_OPT_Size 13 +#define GRCC_OPT_NoMultiEdge 11 +#define GRCC_OPT_SymmInitial 12 +#define GRCC_OPT_SymmFinal 13 +#define GRCC_OPT_Size 14 + +typedef unsigned long ULong; typedef struct { const char *name; @@ -105,6 +105,38 @@ typedef struct { int DUMMYPADDING; } OptDef; +typedef struct { + const char *name; + const char *cname; + const char *mean; +} OptQGDef; + +typedef struct { + const char *name; + int index; + int sign; +} OptQGRef; + +/*--------------------------------------------------------------- + * QGRAF options + */ +#define GRCC_QGRAF_OPT_ONEPI 0 +#define GRCC_QGRAF_OPT_ONSHELL 1 +#define GRCC_QGRAF_OPT_NOSIGMA 2 +#define GRCC_QGRAF_OPT_NOSNAIL 3 +#define GRCC_QGRAF_OPT_NOTADPOLE 4 +#define GRCC_QGRAF_OPT_SIMPLE 5 +#define GRCC_QGRAF_OPT_BIPART 6 +#define GRCC_QGRAF_OPT_CYCLI 7 +#define GRCC_QGRAF_OPT_FLOOP 8 +#define GRCC_QGRAF_OPT_TOPOL 9 + +#ifdef GRCC_QGRAF_OPT_TOPOL +#define GRCC_QGRAF_OPT_Size 10 +#else +#define GRCC_QGRAF_OPT_Size 9 +#endif + /*--------------------------------------------------------------- * Conversion of * eind : index (ed) of EGraph.edges[ed]