neuron 7.5-1 / usr / share / nrn / lib / hoc / binfo.hoc

// like a file in that scanvar returns next vector element
begintemplate VectorStream
public scanvar, i, vec
objref vec
proc init() {
	i = -1
	vec = $o1
}
func scanvar() {
	return vec.x[i += 1]
}
endtemplate VectorStream

// each cell has a gid and each piece has a special spgid which will be
// equal to gid on the piece containing the output port. The idea is that
// from a spgid, one can derive the gid.
// Also the set of sid for each cell is encoded to make one cells set distinct
// from any other cell.

{load_file("stdlib.hoc")}

begintemplate CellBalancePiece
public sid, idvec, cx, host, merge, root
objref idvec, root
proc init() {local i, n
	sid = $2
	host = -1
	cx = $o1.scanvar
	n = $o1.scanvar
	idvec = new Vector()
	for i=0, n-1 {
		idvec.append($o1.scanvar)
	}
}
// all subtrees of the idvec connected to 0 of the first idvec
// $o1 is the cell SectionRef list
proc merge() {local i  localobj sr
if(0) for i=0,idvec.size-1 {
sr = $o1.object(idvec.x[i])
sr.parent printf("original %s connect ", secname())
sr.sec printf("%s(%g), (%g) : %d\n", secname(), section_orientation(), parent_connection(), idvec.x[i])
}
	root = $o1.object(idvec.x[0])
	for i=1, idvec.size-1 {
		sr = $o1.object(idvec.x[i])
//sr.sec {printf("disconnect %d %s\n", i, secname())}
//$o1.object(idvec.x[0]).sec { printf("reconnect to %s(0)\n", secname())}
		sr.sec { disconnect() }
		$o1.object(idvec.x[0]).sec { connect sr.sec(0), 0 }
	}
}
endtemplate CellBalancePiece

begintemplate CellBalanceInfo
public gid,  cx, subtrees, pcx, host, distinct_hosts, nsid, spgid
public multisplit, tmphost, root
objref subtrees, root
proc init() {local isid, nsub, isub
	tmphost = -1 // prevent multiple pieces on one host for trial
	host = -1
	subtrees = new List()
	gid = $o1.scanvar
	spgid = gid
	cx = $o1.scanvar
	nsid = $o1.scanvar
	for isid = 0, nsid-1 {
		nsub = $o1.scanvar
		for isub = 0, nsub-1 {
			subtrees.append(new CellBalancePiece($o1, isid))
		}
	}
	if (subtrees.count == 1) {
//		printf("%d only has 1 piece. Do not split\n", gid)
		nsid = 0
		subtrees.remove_all
	}
}

func pcx() {local i, c
	c = 0
	if (subtrees.count == 0 || $1 == 0) {
		c = cx
	}else for i = 0, subtrees.count-1 {
		c += subtrees.object(i).cx
	}
	return c
}
func distinct_hosts() {local i, j, h
	for i=0, subtrees.count-2 {
		h = subtrees.object(i).host
		for j=i+1, subtrees.count-1 {
			if (h == subtrees.object(j).host) {
				return 0
			}
		}
	}
	return 1
}

// $o1 NetCon with cell output source
// $2 is the binfo msgid.
// $o3 is the ParallelContext
// $4 is the possibly fake pnm.myid
// on entry the entire cell exists, on exit only the subtrees
// for this host. spgid will be gid for the subtree left containing
// the output port and otherwise will be derived from $2 and the subtree index
// The spgid will be registered for this machine and associated with the
// root of the subtree (but without it being an output port). The exception
// is that the subtree with the output port will have gid registered and
// be associated with that output port.
func multisplit() {local i  localobj srout, cell
	$o1.preloc() srout = new SectionRef() pop_section()
	cell = $o1.precell
	// spgid = gid + $2*(subtree_index + 1) except the subtree
	// with the output port has spgid = gid
	msdiv(cell, srout, $2, $o3, $4)
	if (spgid == gid) {
		$o3.cell(gid, $o1, 1)
	}
	return spgid
}

proc msdiv() {local i, sid, msgid, px, cnt, x, max_sid \
  localobj cell, srout, pc, sl, allsr, srlist, sidvec, sr, vsav, vdispl, hosts, nil
	max_sid = 100
	cell = $o1
	srout = $o2
	msgid = $3
	pc = $o4
	allsr = new List()
	sl = new SectionList()
	// if top level cell use currently accessed section
	if (object_id(cell) == 0) {	
		sl.wholetree
	}else{
		if (section_exists("soma", cell)) {
			$o1.soma { sl.wholetree }
		}else{ // at least there should be an all SectionList
			i = 0
			forsec cell.all {
				if (i == 0) {
					sl.wholetree
				}
				i += 1
			}
		}
	}
	forsec sl { allsr.append(new SectionRef()) }
	root = allsr.object(0)
	for i=0, subtrees.count-1 {
		subtrees.object(i).merge(allsr)
	}

	// in case the voltages are encoding information, save them
	vdispl = new Vector(allsr.count)
	cnt = 0
	for i=0, allsr.count-1 allsr.o(i).sec {
		vdispl.x[i] = cnt
		cnt += nseg + 2
	}
	vsav = new Vector(cnt)
	cnt = 0
	for i=0, allsr.count-1 allsr.object(i).sec {
		for (x) {
			vsav.x[cnt] = v(x)
			cnt += 1
		}
	}

	for i=0, allsr.count-1 allsr.object(i).sec {
		v = -1
	}
	srlist = new List()
	sidvec = new Vector()
	hosts = new Vector()
	sid = 0
	for i=1, subtrees.count-1 {
		sr = allsr.object(subtrees.object(i).idvec.x[0])
		sr.sec px = parent_connection()
		sr.parent if (v(px) == -1) {v(px) = sid  sid += 1 }
	}
	if (sid >= max_sid) {
		printf("%d gid %d with %d pieces has %d >= max %d sids\n", pc.id, gid, subtrees.count, sid, max_sid)
		execerror("Cell has too many sids. Increase max_sid in binfo.hoc")
	}
	for (i = subtrees.count-1; i > 0; i -= 1) {
		sr = allsr.object(subtrees.object(i).idvec.x[0])
		srlist.append(sr)
		sr.sec px = parent_connection()
		sr.parent { sidvec.append(gid*max_sid + v(px)) }
		hosts.append(subtrees.object(i).host)
	}
	sidvec.append(-1)
	allsr.object(0).root srlist.append(new SectionRef())
	hosts.append(subtrees.object(0).host)
	
	if (0 && pc.id == 0) for i=0, sidvec.size-1 srlist.object(i).sec {
		printf("%d %ld %d %s\n", i, sidvec.x[i], hosts.x[i], secname())
	}
	multisplit_divide(sidvec, srlist, hosts, $5, pc)

	// restore saved voltages for still existing sections
	for i=0, allsr.count-1 {
		sr = allsr.object(i)
		if (sr.exists) sr.sec {
			cnt = vdispl.x[i]
			for (x) {
				v(x) = vsav.x[cnt]
				cnt += 1
			}
		}
	}

	spgid = -1
	for i=0, srlist.count-1 {
		sr = srlist.object(i)
		if (sr.exists()) {
			spgid = gid + msgid*(i + 1)
			if (srout.exists()) srout.root if (sr.is_cas()) {
				spgid = gid
			}
			pc.set_gid2node(spgid, pc.id)
			if (spgid != gid) sr.sec {
				pc.cell(spgid, new NetCon(&v(.5), nil), 0)
//printf("%d spgid=%d for %s\n", pc.id, spgid, secname())
			} // else pc.cell will be called on outport after return
		}
	}
}

// division information format is a leaf to root ordered list of
// sid sectionname hostid
// where each triple refers to the subtree for which the triple defines
// the root of the subtree.
// Note that root of the whole cell defines a subtree which is what is
// left over when all other subtrees have been disconnected. Some point
// on this subtree (typically for convenience, soma(0)) should be the
// last item in the list and may have a sid of -1. In fact the only information
// used in the last item is the host field which defines where it will
// continue to exists.
// Note that consistency requires that if two pairs have the same parent
// then they must also have the same sid.
// Also required is that there cannot be more than two distinct connection
// points into the root subtree from other subtrees. Furthermore, there
// can only be at  most one connection point into a non-root subtree from
// other (child) subtrees.

// args are sid Vector, SectionRefList, host vector
// Note that any subtree with a host != pc.id is deleted.

// This algorithm allows multiple subtrees from the same cell on
// one cpu. In fact all the subtrees can be on one cpu.

// Note: occasionally the root  parent piece (containing the soma)
// is an orphan at its sid == 0 (root)
// in the sense that
// no other piece connects to that location. This is naturally handled
// in that we do not call multisplit at that point and so there is
// no wasteful backbone.


proc multisplit_divide() {local i, x \
 localobj sids, subroots, subroot, hosts, sr, sl1, sl2, pc, parents, px
	pc = $o5
	sids = $o1
	subroots = $o2
	hosts = $o3
	// sl1 is the whole tree
	sl1 = new SectionList()
	$o2.object(0).sec { sl1.wholetree }
	// create a parallel list of true_parent SectionRef and parent
	// connection points. We use that to help construct the parents
	// parent to child multisplit connection near the end.
	parents = new List()
	px = new Vector()
	for i=0, subroots.count - 2 { // do not use the root subtree
		subroot = subroots.object(i)
		if (subroot.has_parent) { // if not then root(0)
			subroot.parent parents.append(new SectionRef())
			subroot.sec px.append(parent_connection())
		}else{
			subroot.root parents.append(new SectionRef())
			px.append(0)
		}
	}	
	// now disconnect leaving only subtrees
	for i=0, subroots.count - 2 {
		subroots.object(i).sec disconnect()
	}
	// delete any subtree not on this cpu.
	for i=0, subroots.count - 1 {
		if (hosts.x[i] != $4) {
			sl1 = new SectionList()
			subroots.object(i).sec sl1.wholetree()
			forsec sl1 delete_section()
		}
	}
	// now we can do the multisplit for each subroot
	for i=0, subroots.count - 2 {
		if (subroots.object(i).exists()) subroots.object(i).sec {
//printf("%d %s pc.multisplit(%g, %d, 2)\n", pc.id, secname(), 0, sids.x[i])
			pc.multisplit(0, sids.x[i], 2)
			// and it will be helpful to initialize v marks
			sl1 = new SectionList()
			sl1.wholetree()
			forsec sl1 v = -1
		}
	}
	// mark all the subtree parents
	for i=0, subroots.count - 2 {
		if (parents.object(i).exists()) parents.object(i).sec {
			v(px.x[i]) = -1
		}
	}
	// and we can do the multisplit for the parents
	// but we should only do any given parent location once
	// (that is why we find the v mark useful)
	for i=0, subroots.count - 2 {
		if (parents.object(i).exists()) parents.object(i).sec {
			if (v(px.x[i]) == -1) { // multisplit needs to be done
//printf("%d %s pc.multisplit(%g, %d, 2)\n", pc.id, secname(), px.x[i], sids.x[i])
				pc.multisplit(px.x[i], sids.x[i], 2)
				v(px.x[i]) = sids.x[i]
			}else if (v(px.x[i]) != sids.x[i]) { // sanity check
printf("i=%d px.x[i]=%d v=%g sids.x[i]=%g\n", i, px.x[i], v(px.x[i]), sids.x[i])
				execerror("Subtrees at same parent with different sid", "")
			}
		}
	}
}

endtemplate CellBalanceInfo

begintemplate BalanceInfo
public bilist, cx, npiece, locality_balance
public nhost, stat, basename, write_balhost, msgid, write_hostcontext
public items, gids, sindices, mymetis2, mybal, cbindex
public base_gid, thishost_gid, gid2cx, write_colgid, ihost
objref bilist
strdef basename
objref items, gids, sindices, cbindex // when only a few are read

proc init() {
	msgid = -1
	nhost = -1
	ihost = -1
	bilist = new List()
	if (numarg() == 0) {
		return
	}
	if (numarg() > 0) {
		basename = $s1
	}
	if (numarg() == 1) {
		read_all()
	}else if (numarg() == 3) {
		read_host($2, $3, 1)
	}else if (numarg() == 4) {
		// balance file must have nhost*nthread rank lines
		// read nthread contiguous lines starting at rank*nthread
		read_host($2, $3, $4) //rank, nhost, nthread
	}
}

proc read_all() {local i1, i2, n1, n2  localobj f, s
	f = new File()
	s = new String()
	sprint(s.s, "%s.dat", basename)
	f.ropen(s.s)
	n1 = f.scanvar
	for i1=0, n1 - 1 {
		n2 = f.scanvar
		for i2=0, n2 -1 {
			bilist.append(new CellBalanceInfo(f))
		}
	}
	f.close()
}

proc read_host() {local i, j, n, k, i1, i2, n1, n2, rank, np, nthread, md, spgid \
  localobj f, s, cb, pc
	pc = new ParallelContext()
	md = 2^31 - 1 // maximum possible gid
	rank = $1
	np = $2
	nthread = $3
	items = new Vector()
	gids = new Vector()
	sindices = new Vector()

	f = new File()
	s = new String()
	sprint(s.s, "%s.%d.dat", basename, np*nthread)
	f.ropen(s.s)
	msgid = f.scanvar
	nhost = f.scanvar
	if (nhost != np*nthread) { execerror("read_host() wrong number of lines", "") }
	ihost = rank
	for i=0, rank*nthread-1 {f.gets(s.s)}
    for j=0, nthread-1 {
	if (f.scanvar != (rank*nthread + j)) { execerror("read_host() format error", "") }
	n = f.scanvar // number of triples
	for i=0, n-1 {
		items.append(f.scanvar)
		gids.append(f.scanvar)
		sindices.append(f.scanvar)
	}
    }
	cbindex = new Vector(n)
	f.close
if (0) {
printf("msgid=%d nhost=%d ihost=%d\n", msgid, nhost, ihost)
print "items" items.printf
print "gids" gids.printf
print "sindices" sindices.printf
}
	if (items.size == 0) { return }
	i = 0
	k = 0
	sprint(s.s, "%s.dat", basename)
	f.ropen(s.s)
	n1 = f.scanvar
	for i1=0, n1 - 1 {
		n2 = f.scanvar
		for i2=0, n2 - 1 {
			if (k == items.x[i]) {
				cb = new CellBalanceInfo(f)
				cbindex.x[i] = bilist.count
				bilist.append(cb)
//printf("k=%d items.x[%d]=%d gid=%d %d\n", k, i, items.x[i], gids.x[i], \
//bilist.object(i).gid)
				spgid = cb.gid + msgid*(cb.subtrees.count + 1)
				if (spgid > md) {
printf("%d Cell %d with %d pieces has spgid=%ld > 2^31 - 1. Decrease msgid %d in %s.%d.dat file\n", pc.id, cb.gid, cb.subtrees.count, spgid, msgid, basename, np*nthread)
					execerror("msgid too large relative to number of pieces of cell")
				}

				while(1) {
					if (gids.x[i] != cb.gid) {
						execerror("gid inconsistency in read_host ", "")
					}
					cb.host = ihost
					if (cb.subtrees.count > 0) {
		cb.subtrees.object(sindices.x[i]).host = ihost
					}
					i += 1
					if ( i >= items.size) {
						f.close
						return
					}
					if (k != items.x[i]) { break }
				}
			}else{
				s = new CellBalanceInfo(f) // skip
			}
			k += 1
		}
	}
	f.close()
}

// legacy name
proc mymetis2() {
	nhost = $1
	balance($1, 0, bilist)
}

// $1 total number of cpus
// $2 cpu group size
// Return balance
// whole cells are used to balance the $1/$2 group partitions
// and then the cpus in the groups are balanced using the cell
// pieces.
func locality_balance() {local i, j, gsize, ngroup, pw, mw, cavg \
  localobj clist, cb, wt, ipart, wpart, indices
	nhost = $1
	gsize = $2
	ngroup = $1/$2
	wt = new Vector()
	for i=0, bilist.count-1 {
		cb = bilist.o(i)
		if (cb.subtrees.count == 0) {
			wt.append(cb.cx)
		}else{
			wt.append(0)
			for j=0, cb.subtrees.count-1 {
				wt.x[i] += cb.subtrees.object(j).cx
			}
		}
	}
	cavg = wt.sum/nhost
	lpt(wt, ngroup, ipart, wpart)
	indices = new Vector()
	clist = new List()
	mw = 0
	for i=0, ngroup-1 { // should parallelize this
		indices.indvwhere(ipart, "==", i)
		clist.remove_all()
		for j=0, indices.size-1 {
			clist.append(bilist.o(indices.x[j]))
		}
		pw = balance(gsize, i*gsize, clist)
		if (pw > mw) { mw = pw }
	}
	return mw/cavg
}

// $1 is the number of cpus in the group
// $2 is the cpu offset for the group
// $o3 is the list of CellBalance objects
// Returns the max partition weight
// and fills the CellBalance objects (or pieces if they
// exist, with the hostid (cpu index + offset)
func balance() {local i, j, k, ncpu, mw, coff \
  localobj wt, cb, cpu, wcpu
	ncpu = $1
	coff = $2
	wt = new Vector()
	for i=0, $o3.count-1 {
		cb = $o3.object(i)
		if (cb.subtrees.count == 0) {
			wt.append(cb.cx)
		} else for j=0, cb.subtrees.count-1 {
			wt.append(cb.subtrees.object(j).cx)
		}
	}
//	printf("%d pieces  size max=%g min=%g\n", wt.size, wt.max, wt.min)
	mw = lpt(wt, ncpu, cpu, wcpu)
//printf("npiece = %d   npart = %d\n", wt.size, $1)
//printf("max and min complexity %g %g   avg = %g\n", wcpu.max, wcpu.min, wcpu.mean)
//printf("load balance %g\n", b)
	k = -1
	for i=0, $o3.count-1 {
		cb = $o3.object(i)
		if (cb.subtrees.count == 0) {
			cb.host = cpu.x[k += 1] + coff
		} else for j=0, cb.subtrees.count-1 {
			cb.subtrees.object(j).host = cpu.x[k += 1] + coff
		}
	}
	return mw
}

// least processing time algorithm
// $o1 is vector of weights  $2 is number of partitions
// return value is the max partition weight, return in $o3 a
// partition indices parallel to weights($o1), return in $o4
// a vector of partition weights
func lpt() {local i, j  localobj wx, ix, pw
	wx = $o1.sortindex.reverse
	ix = new Vector($o1.size)
	pw = new Vector($2)
	for i=0, $o1.size-1 {
		j = wx.x[i]
		w = $o1.x[j]
		ip = pw.min_ind
		pw.x[ip] += w
		ix.x[j] = ip
	}
	$o3 = ix
	$o4 = pw
	return pw.max
}

func cx() {local i, c, j
	j = 0
	if (numarg() == 1) {
		j = $1
	}
	c = 0
	for i = 0, bilist.count-1 {
		c += bilist.object(i).pcx(j)
	}
	return c
}
func npiece() {local i, n
	n = 0
	for i=0, bilist.count-1 {
		if ( bilist.object(i).subtrees.count == 0) {
			n += 1
		}else{
			n += bilist.object(i).subtrees.count
		}
	}
	return n
}
proc stat() {local i, j, c, cmax, mcp, mxs, mcpi, mxsi  localobj cxvec, npvec, cb, p
	if (nhost == -1) return
	mcp = 1
	mxs = 0
	cxvec = new Vector(nhost)
	npvec = new Vector(nhost)
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.nsid > mxs) {
			mxs = cb.nsid
			mxsi = i
		}
		if (cb.subtrees.count > mcp) {
			mcp = cb.subtrees.count
			mcpi = i
		}
		if (cb.subtrees.count == 0) {
			cxvec.x[cb.host] += cb.cx
			npvec.x[cb.host] += 1
		}else for j=0, cb.subtrees.count-1 {
			p = cb.subtrees.object(j)
			cxvec.x[p.host] += p.cx
			npvec.x[p.host] += 1
		}
	}
	c = cx()
	cmax = cxvec.max
	printf("total complexity is %g\n", c)
	printf("%d cells\n", bilist.count)
	printf("%d pieces\n", npiece())
	printf("maximum complexity is %g for host %g\n", cmax, cxvec.max_ind)
	printf("load imbalance is %.1f%%\n", 100*cmax*nhost/c - 100)
	printf("maximum of %d pieces on host %g\n", npvec.max, npvec.max_ind)
	if (mcp > 1) {
		printf("at least one cell is broken into %d pieces (bilist[%d], gid %d)\n", mcp, mcpi, bilist.object(mcpi).gid)
		printf("at least one cell has %d sids (bilist[%d], gid %d)\n", mxs, mxsi, bilist.object(mxsi).gid)
	}else{
		printf("no broken cells\n")
	}
}

// write a host file with the format
// max_gid
// nhost
// ihost ngid for that host followed by a list of (index, gid, subtreeindex) triples
// the indices will be in order and the gid is just for error checking
// the idea is to allow the parallel run to not have to save all the info
proc write_balhost() {local i, j, mg  localobj hostgids, cb, p, f, s
	if (nhost == -1) return
	hostgids = new List()
	for i=0, nhost-1 { hostgids.append(new Vector()) }
	mg = -1
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.gid > mg) {
			mg = cb.gid
		}
	}
	for (msgid=10; msgid <= mg; msgid *= 10) {}// encode gid in spgid along with subtree index
	if (msgid < 1e7) {
		msgid = 1e7 // there may be unused gids in BlueBrain
	}
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.subtrees.count == 0) {
			hostgids.object(cb.host).append(i)	
			hostgids.object(cb.host).append(0)
		}else for j=0, cb.subtrees.count-1 {
			p = cb.subtrees.object(j)
			hostgids.object(p.host).append(i)
			hostgids.object(p.host).append(j)
		}
	}
	s = new String()
	sprint(s.s, "%s.%d.dat", basename, nhost)
	f = new File()
	f.wopen(s.s)
	f.printf("msgid %d\n", msgid)
	f.printf("nhost %d\n", hostgids.count)
	for i=0, hostgids.count-1 {
		p = hostgids.object(i)
		f.printf("%d %d", i, p.size/2)
		for (j = 0; j < p.size; j += 2) {
f.printf("  %d %d %d", p.x[j], bilist.object(p.x[j]).gid, p.x[j+1])
		}
		f.printf("\n")
	}
	f.close
}

// derived from write_balhost but output format is suitable for
// the fast_create.hoc file

proc write_colgid() {local i, j, mg  localobj hostgids, cb, p, f, s
	if (nhost == -1) return
	hostgids = new List()
	for i=0, nhost-1 { hostgids.append(new Vector()) }
	mg = -1
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.gid > mg) {
			mg = cb.gid
		}
	}
	for (msgid=10; msgid <= mg; msgid *= 10) {}// encode gid in spgid along with subtree index
	if (msgid < 1e7) {
		msgid = 1e7 // there may be unused gids in BlueBrain
	}
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.subtrees.count == 0) {
			hostgids.object(cb.host).append(i)	
			hostgids.object(cb.host).append(0)
		}else for j=0, cb.subtrees.count-1 {
			p = cb.subtrees.object(j)
			hostgids.object(p.host).append(i)
			hostgids.object(p.host).append(j)
		}
	}
	s = new String()
	sprint(s.s, "%s.%d.dat", basename, nhost)
	f = new File()
	f.wopen(s.s)
	f.printf("ncell %d\nnhost %d\n", bilist.count, hostgids.count)
	for i=0, hostgids.count-1 {
		p = hostgids.object(i)
		for (j = 0; j < p.size; j += 2) {
f.printf("%d %d\n", bilist.object(p.x[j]).gid, i)
		}
	}
	f.close
}

// just enough basename.ncpu.dat to get complete host communication for
// specified host arg. 2nd arg is the nhost for the original basename.nhost.dat
// file
// $3 is the level

proc write_hostcontext() {local i, j, il, ngid, gid, host  \
    localobj f, s, cb, hostgids, gvec, mark, marked, gidvec, hostvec, gi
	//first, read the complete basename.$2.dat file
	f = new File()
	s = new String()
	sprint(s.s, "%s.%d.dat", basename, $2)
	f.ropen(s.s)
	msgid = f.scanvar
	nhost = f.scanvar
	gidvec = new Vector()
	hostvec = new Vector()
	hostgids = new List()
	mark = new Vector(nhost)
	for i=0, nhost-1 {
		if (i != f.scanvar) { execerror("bad format for ", s.s) }
		gvec = new Vector()
		hostgids.append(gvec)
		ngid = f.scanvar
		for j=0, ngid-1 {
			gvec.append(f.scanvar)
			gid = f.scanvar
			gvec.append(gid)
			gvec.append(f.scanvar)
			gidvec.append(gid)
			hostvec.append(i)
		}
	}
	// now mark the hosts we need, ie. every host that has a gid
	// referred to by $1
	mark.x[$1] = 1
	gvec = hostgids.object($1)
	domark(gvec, gidvec, hostvec, mark)
	for il=1, $3-1 {
		// try the second level
		gi = mark.c.indvwhere("==", 1)
		for i=0, gi.size-1 {
			domark(hostgids.object(gi.x[i]), gidvec, hostvec, mark)
		}
	}

	// print the reduced basename.ncpu.dat file
	marked = mark.c.indvwhere("==", 1)
	sprint(s.s, "%s.%d.dat", basename, marked.size)
	f.wopen(s.s)
	f.printf("msgid %d\n", msgid)
	f.printf("nhost %d\n", marked.size)
	for i=0, marked.size-1 {
		j = marked.x[i]
		gvec = hostgids.object(j)
		f.printf("%d %d", i, gvec.size/3)
//		f.printf("%d %d", marked.x[i], gvec.size/3)
		for (j=0; j < gvec.size; j += 3) {
			f.printf("  %d %d %d", gvec.x[j], gvec.x[j+1], gvec.x[j+2])
		}
		f.printf("\n")
	}
print "wrote ", s.s
}

proc domark() {local i, j, gid  localobj gvec, gidvec, hostvec, mark, gi
	gvec = $o1  gidvec = $o2  hostvec=$o3  mark=$o4
	for (i=0; i < gvec.size-1; i += 3) {
		gid = gvec.x[i+1]
		gi = gidvec.c.indvwhere("==", gid)
		for j=0, gi.size-1 {
			mark.x[hostvec.x[gi.x[j]]] = 1
		}
	}
}

func base_gid() {
	return $1 % msgid
}
func thishost_gid() {local i
	if ($1 >= msgid) { return $1 }
	i = gids.indwhere("==", $1)
	if (i < 0) { return -1 }
	return bilist.object(cbindex.x[i]).spgid
}
func gid2cx() {local i localobj cb
	for i=0, bilist.count-1 {
		cb = bilist.object(i)
		if (cb.spgid == $1) {
			return cb.cx
		}
	}
	return 0
}
endtemplate BalanceInfo

proc mymetis2() { localobj bi
	bi = new BalanceInfo($s1)
	bi.mymetis2($2)
	bi.stat()
	bi.write_colgid()
}

proc mymetis3() { localobj bi
	bi = new BalanceInfo($s1)
	bi.mymetis2($2)
	bi.stat()
	bi.write_balhost()
}

proc locality_balance() { localobj bi
	bi = new BalanceInfo($s1)
	bi.locality_balance($2, $3)
	bi.stat()
	bi.write_balhost()
}

/*
objref bi

if (0) {
bi = new BalanceInfo("cx")
bi.mymetis2(8192)
}
*/