adding chaco partionner to contrib

contrib/Chaco/Makefile

+DEST_DIR = 	../exec
+DEST=		${DEST_DIR}/chaco
+CC = 		gcc
+IFLAG =		-Imain
+CFLAGS =	-O2
+OFLAGS =	-O2
+
+FILES.c=	main/user_params.c main/interface.c main/main.c \
+		submain/balance.c submain/divide.c submain/submain.c \
+		input/input_assign.c \
+		input/check_input.c input/input.c input/input_geom.c \
+		input/input_graph.c input/read_params.c input/reflect_input.c \
+		input/read_val.c \
+		graph/check_graph.c graph/free_graph.c \
+		graph/reformat.c graph/subgraph.c graph/graph_out.c \
+		inertial/eigenvec2.c inertial/eigenvec3.c inertial/inertial.c \
+		inertial/inertial1d.c inertial/inertial2d.c \
+		inertial/inertial3d.c inertial/make_subgeom.c \
+		klspiff/buckets.c klspiff/buckets_bi.c klspiff/buckets1.c \
+		klspiff/bilistops.c klspiff/coarsen_kl.c klspiff/count_weights.c \
+		klspiff/compress_ewgts.c klspiff/kl_init.c klspiff/kl_output.c \
+		klspiff/klspiff.c klspiff/make_bndy_list.c \
+		klspiff/make_kl_list.c klspiff/nway_kl.c \
+		klvspiff/bpm_improve.c klvspiff/bucketsv.c \
+		klvspiff/clear_dvals.c klvspiff/coarsen_klv.c \
+		klvspiff/countup_vtx_sep.c klvspiff/find_bndy.c klvspiff/flow.c \
+		klvspiff/klv_init.c klvspiff/klvspiff.c klvspiff/make_bpgraph.c \
+		klvspiff/make_sep_list.c klvspiff/matching.c klvspiff/nway_klv.c \
+		klvspiff/flatten.c \
+		coarsen/coarsen.c coarsen/interpolate.c coarsen/makefgraph.c \
+		coarsen/makeccoords.c \
+		coarsen/coarsen1.c coarsen/makev2cv.c \
+		coarsen/maxmatch.c coarsen/maxmatch1.c coarsen/maxmatch2.c \
+		coarsen/maxmatch3.c coarsen/maxmatch4.c coarsen/maxmatch5.c \
+		connect/add_edges.c connect/connected.c connect/find_edges.c \
+		eigen/bidir.c eigen/bisect.c eigen/checkeig.c \
+		eigen/checkeig_ext.c \
+		eigen/checkorth.c eigen/cksturmcnt.c eigen/mkeigvecs.c\
+		eigen/eigensolve.c eigen/get_extval.c eigen/get_ritzvals.c \
+		eigen/lanczos_FO.c eigen/lanczos_SO.c eigen/lanczos_SO_float.c \
+		eigen/lanczos_ext.c eigen/lanczos_ext_float.c eigen/lanc_seconds.c\
+		eigen/lanpause.c eigen/makeorthlnk.c eigen/mkscanlist.c \
+		eigen/orthog1.c eigen/orthogonalize.c eigen/orthogvec.c \
+		eigen/ql.c eigen/rqi.c eigen/rqi_ext.c eigen/scale_diag.c \
+		eigen/scanmax.c eigen/scanmin.c eigen/solistout.c \
+                eigen/sorthog.c eigen/splarax.c eigen/sturmcnt.c \
+		eigen/Tevec.c eigen/tri_solve.c eigen/warnings.c \
+		symmlq/aprod.c symmlq/msolve.c symmlq/pow_dd.c \
+		symmlq/symmlq.c symmlq/symmlqblas.c \
+		tinvit/tinvit.c tinvit/pythag.c tinvit/epslon.c \
+		optimize/determinant.c optimize/func2d.c \
+		optimize/func3d.c optimize/opt2d.c optimize/opt3d.c \
+		assign/assign.c assign/assign_out.c assign/mapper.c \
+		assign/median.c assign/merge_assign.c \
+		assign/rec_median.c assign/rotate.c assign/y2x.c \
+		bpmatch/checkbp.c bpmatch/inits2d.c bpmatch/inits3d.c \
+		bpmatch/genvals2d.c bpmatch/genvals3d.c bpmatch/map2d.c \
+		bpmatch/map3d.c bpmatch/movevtxs.c \
+		bpmatch/sorts2d.c bpmatch/sorts3d.c \
+		refine_map/compute_cube_edata.c refine_map/compute_cube_vdata.c \
+		refine_map/refine_cube.c refine_map/update_cube_edata.c \
+		refine_map/update_cube_vdata.c refine_map/find_edge_cube.c \
+		refine_map/init_cube_edata.c refine_map/compute_mesh_edata.c \
+		refine_map/compute_mesh_vdata.c refine_map/find_edge_mesh.c \
+		refine_map/init_mesh_edata.c refine_map/refine_mesh.c \
+		refine_map/update_mesh_edata.c refine_map/update_mesh_vdata.c \
+		refine_map/refine_map.c refine_map/make_comm_graph.c \
+		refine_part/refine_part.c refine_part/kl_refine.c \
+		refine_part/make_maps_ref.c refine_part/make_terms_ref.c \
+		internal/force_internal.c internal/improve_internal.c \
+		internal/check_internal.c \
+		misc/define_subcubes.c misc/define_submeshes.c \
+		misc/divide_procs.c misc/merge_goals.c misc/make_term_props.c \
+		misc/count.c misc/countup.c misc/countup_cube.c \
+		misc/countup_mesh.c misc/make_subgoal.c \
+		misc/find_maxdeg.c misc/make_maps.c misc/make_setlists.c \
+		misc/sequence.c misc/perturb.c misc/simple_part.c \
+		misc/time_kernels.c misc/timing.c \
+		util/affirm.c util/array_alloc_2D.c util/bit_reverse.c \
+		util/checkpnt.c util/cpvec.c util/dot.c \
+		util/doubleout.c util/input_int.c util/gray.c \
+		util/machine_params.c util/makevwsqrt.c util/mkvec.c util/norm.c \
+                util/normalize.c util/mergesort.c \
+                util/randomize.c util/smalloc.c util/bail.c \
+		util/scadd.c util/seconds.c util/setvec.c util/shell_sort.c \
+		util/strout.c util/tri_prod.c util/true_or_false.c \
+		util/update.c  util/vecout.c util/vecran.c \
+                util/vecscale.c 
+
+FILES.o=	$(FILES.c:.c=.o) 
+
+
+${DEST}:	${FILES.o} Makefile
+		${CC} ${OFLAGS} ${FILES.o} -lm -o ${DEST}
+
+lint:
+		lint ${IFLAG} ${FILES.c} -lm
+
+alint:
+		alint ${IFLAG} ${FILES.c} -lm
+
+clean:
+		rm -f */*.o ${DEST_DIR}/core
+
+.c.o:
+		${CC} -c ${IFLAG} ${CFLAGS} -o $*.o $*.c

contrib/Chaco/README

+This directory contains all the necessary source code files
+for Chaco.  Simply type "make", and the code will be compiled
+and the executable image placed in "../exec/chaco".  
+
+The default values for the user accessible parameters discussed 
+in the user's guide "doc/User_guide.ps" can be modified by 
+editing the file "main/user_params.c".  You must recompile for 
+these changes to take affect. 
+
+You can also change the code parameters at run time (i.e. without 
+recompiling) by including a file "User_Params" in whichever 
+directory you are running the code from (so in the file 
+"exec/User_Params" in our orignal set-up). Parameter changes 
+can be specified in a natural way in this file using the rules 
+described in the user's guide.

contrib/Chaco/assign/assign.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"structs.h"
+
+
+void      assign(graph, yvecs, nvtxs, ndims, cube_or_mesh, nsets, wsqrt, sets,
+		           active, mediantype, goal, vwgt_max)
+struct vtx_data **graph;	/* data structure with vtx weights */
+double  **yvecs;		/* ptr to list of y-vectors (lengths nvtxs+1) */
+int       nvtxs;		/* number of vertices in graph */
+int       ndims;		/* number of vectors for dividing */
+int       cube_or_mesh;		/* 0 => hypercube, d => d-dimensional mesh */
+int       nsets;		/* number of sets to divide into */
+double   *wsqrt;		/* sqrt of vertex weights */
+short    *sets;			/* processor assignment for my vtxs */
+int      *active;		/* space for nvtxs integers */
+int       mediantype;		/* which partitioning strategy to use */
+double   *goal;			/* desired set sizes */
+int       vwgt_max;		/* largest vertex weight */
+{
+    extern int DEBUG_TRACE;	/* trace execution path of code */
+    extern int DEBUG_ASSIGN;	/* turn on debugging in assignment */
+    double    theta, phi, gamma;/* angles for optimal rotation */
+    double    temp;
+    int       using_vwgts;	/* are vertex weights active? */
+    double    tri_prod();
+    double    opt2d();
+    void      y2x(), mapper(), rotate2d(), opt3d(), rotate3d();
+
+    if (DEBUG_TRACE > 0) {
+	printf("<Entering assign, nvtxs = %d, ndims = %d>\n", nvtxs, ndims);
+    }
+
+    using_vwgts = (vwgt_max != 1);
+
+    if (ndims == 1) {
+	/* Unscale yvecs to get xvecs. */
+	y2x(yvecs, ndims, nvtxs, wsqrt);
+
+	mapper(graph, yvecs, nvtxs, active, sets, ndims, cube_or_mesh, nsets,
+	       mediantype, goal, vwgt_max);
+    }
+
+    else if (ndims == 2) {
+	theta = opt2d(graph, yvecs, nvtxs, nvtxs);
+
+	rotate2d(yvecs, nvtxs, theta);
+
+	y2x(yvecs, ndims, nvtxs, wsqrt);
+
+	mapper(graph, yvecs, nvtxs, active, sets, ndims, cube_or_mesh, nsets,
+	       mediantype, goal, vwgt_max);
+    }
+
+    else if (ndims == 3) {
+
+	if (DEBUG_ASSIGN > 0) {
+	    temp = tri_prod(yvecs[1], yvecs[2], yvecs[3], wsqrt, nvtxs);
+	    printf("Before rotation, 3-way orthogonality = %e\n", temp);
+	}
+
+	opt3d(graph, yvecs, nvtxs, nvtxs, wsqrt, &theta, &phi, &gamma, using_vwgts);
+
+	rotate3d(yvecs, nvtxs, theta, phi, gamma);
+
+	if (DEBUG_ASSIGN > 0) {
+	    temp = tri_prod(yvecs[1], yvecs[2], yvecs[3], wsqrt, nvtxs);
+	    printf("After rotation (%f,%f,%f), 3-way orthogonality = %e\n",
+		   theta, phi, gamma, temp);
+	}
+
+	y2x(yvecs, ndims, nvtxs, wsqrt);
+
+	mapper(graph, yvecs, nvtxs, active, sets, ndims, cube_or_mesh, nsets,
+	       mediantype, goal, vwgt_max);
+    }
+}

contrib/Chaco/assign/assign_out.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+
+static void assign_out_normal(), assign_out_inv();
+
+void      assign_out(nvtxs, sets, nsets, outname)
+int       nvtxs;		/* number of vertices to output */
+short    *sets;			/* values to be printed */
+int       nsets;		/* number of sets */
+char     *outname;		/* name of output file */
+{
+    extern int OUT_ASSIGN_INV;	/* print assignment in inverted form? */
+
+    if (OUT_ASSIGN_INV) {
+	assign_out_inv(nvtxs, sets, nsets, outname);
+    }
+    else {
+	assign_out_normal(nvtxs, sets, outname);
+    }
+}
+
+static void assign_out_normal(nvtxs, sets, outname)
+int       nvtxs;		/* number of vertices to output */
+short    *sets;			/* values to be printed */
+char     *outname;		/* name of output file */
+{
+    FILE     *fout;		/* output file */
+    int       i;		/* loop counter */
+
+    /* Print assignment in simple format. */
+
+    if (outname != NULL) {
+        fout = fopen(outname, "w");
+    }
+    else {
+	fout = stdout;
+    }
+
+    for (i = 1; i <= nvtxs; i++) {
+	fprintf(fout, "%d\n", sets[i]);
+    }
+
+    if (outname != NULL) {
+        fclose(fout);
+    }
+}
+
+
+static void assign_out_inv(nvtxs, sets, nsets, outname)
+int       nvtxs;		/* number of vertices to output */
+short    *sets;			/* values to be printed */
+int       nsets;		/* number of sets */
+char     *outname;		/* name of output file */
+{
+    FILE     *fout;		/* output file */
+    int      *size;		/* # vtxs in sets / index into inorder */
+    int      *inorder;		/* list of vtxs in each set */
+    int       i, j;		/* loop counter */
+    double   *smalloc();
+
+    /* Print assignment in inverted format. */
+    /* For each set, print # vertices, followed by list. */
+
+    if (outname != NULL) {
+        fout = fopen(outname, "w");
+    }
+    else {
+	fout = stdout;
+    }
+
+    size = (int *) smalloc((unsigned) (nsets + 1) * sizeof(int));
+    inorder = (int *) smalloc((unsigned) nvtxs * sizeof(int));
+
+    for (j = 0; j < nsets; j++)
+	size[j] = 0;
+    for (i = 1; i <= nvtxs; i++)
+	++size[sets[i]];
+
+    /* Modify size to become index into vertex list. */
+    for (j = 1; j < nsets; j++) {
+	size[j] += size[j - 1];
+    }
+    for (j = nsets - 1; j > 0; j--) {
+	size[j] = size[j - 1];
+    }
+    size[0] = 0;
+
+    for (i = 1; i<= nvtxs; i++) {
+	j = sets[i];
+	inorder[size[j]] = i;
+	++size[j];
+    }
+
+    /* The inorder array now clumps all the vertices in each set. */
+    /* Now reconstruct size array to index into inorder. */
+    for (j = nsets - 1; j > 0; j--) {
+	size[j] = size[j - 1];
+    }
+    size[0] = 0;
+    size[nsets] = nvtxs;
+
+    /* Now print out the sets. */
+    for (j = 0; j < nsets; j++) {
+	fprintf(fout, "  %d\n", size[j + 1] - size[j]);
+	for (i = size[j]; i < size[j + 1]; i++) {
+	    fprintf(fout, "%d\n", inorder[i]);
+	}
+    }
+
+    if (outname != NULL) {
+        fclose(fout);
+    }
+}

contrib/Chaco/assign/mapper.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"params.h"
+#include	"defs.h"
+#include	"structs.h"
+
+
+void      mapper(graph, xvecs, nvtxs, active, sets, ndims, cube_or_mesh, nsets,
+		           mediantype, goal, vwgt_max)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double  **xvecs;		/* continuous indicator vectors */
+int       nvtxs;		/* number of vtxs in graph */
+int      *active;		/* space for nmyvals ints */
+short    *sets;			/* returned processor assignment for my vtxs */
+int       ndims;		/* number of dimensions being divided into */
+int       cube_or_mesh;		/* 0 => hypercube, d => d-dimensional mesh */
+int       nsets;		/* number of sets to divide into */
+int       mediantype;		/* type of eigenvector partitioning to use */
+double   *goal;			/* desired set sizes */
+int       vwgt_max;		/* largest vertex weight */
+{
+    double    temp_goal[2];	/* combined set goals if using option 1. */
+    double    wbelow;		/* weight of vertices with negative values */
+    double    wabove;		/* weight of vertices with positive values */
+    short    *temp_sets;	/* sets vertices get assigned to */
+    int       vweight;		/* weight of a vertex */
+    int       using_vwgts;	/* are vertex weights being used? */
+    int       bits;		/* bits for assigning set numbers */
+    int       i, j;		/* loop counters */
+    int       sfree();
+    double   *smalloc();
+    void      median(), rec_median_k(), rec_median_1(), median_assign();
+    void      map2d(), map3d();
+
+    /* NOTE: THIS EXPECTS XVECS, NOT YVECS! */
+
+    using_vwgts = (vwgt_max != 1);
+
+    if (ndims == 1 && mediantype == 1)
+	mediantype = 3;		/* simpler call than normal option 1. */
+
+    if (mediantype == 0) {	/* Divide at zero instead of median. */
+	bits = 1;
+	temp_sets = (short *) smalloc((unsigned) (nvtxs + 1) * sizeof(short));
+	for (j = 1; j <= nvtxs; j++)
+	    sets[j] = 0;
+
+	for (i = 1; i <= ndims; i++) {
+	    temp_goal[0] = temp_goal[1] = 0;
+	    for (j = 0; j < (1 << ndims); j++) {
+		if (bits & j)
+		    temp_goal[1] += goal[j];
+		else
+		    temp_goal[0] += goal[j];
+	    }
+	    bits <<= 1;
+
+	    wbelow = wabove = 0;
+	    vweight = 1;
+	    for (j = 1; j <= nvtxs; j++) {
+		if (using_vwgts)
+		    vweight = graph[j]->vwgt;
+		if (xvecs[i][j] < 0)
+		    wbelow += vweight;
+		else if (xvecs[i][j] > 0)
+		    wabove += vweight;
+	    }
+
+	    median_assign(graph, xvecs[i], nvtxs, temp_goal, using_vwgts, temp_sets,
+			  wbelow, wabove, (double) 0.0);
+
+	    for (j = 1; j <= nvtxs; j++)
+		sets[j] = (sets[j] << 1) + temp_sets[j];
+	}
+    }
+
+    else if (mediantype == 1) {	/* Divide using min-cost assignment. */
+	if (ndims == 2)
+	    map2d(graph, xvecs, nvtxs, sets, goal, vwgt_max);
+	else if (ndims == 3)
+	    map3d(graph, xvecs, nvtxs, sets, goal, vwgt_max);
+    }
+
+
+    else if (mediantype == 2) {	/* Divide recursively using medians. */
+	rec_median_k(graph, xvecs, nvtxs, active, ndims, cube_or_mesh, goal,
+		     using_vwgts, sets);
+    }
+
+    else if (mediantype == 3) {	/* Cut with independent medians => unbalanced. */
+	bits = 1;
+	temp_sets = (short *) smalloc((unsigned) (nvtxs + 1) * sizeof(short));
+	for (j = 1; j <= nvtxs; j++)
+	    sets[j] = 0;
+
+	for (i = 1; i <= ndims; i++) {
+	    temp_goal[0] = temp_goal[1] = 0;
+	    for (j = 0; j < (1 << ndims); j++) {
+		if (bits & j)
+		    temp_goal[1] += goal[j];
+		else
+		    temp_goal[0] += goal[j];
+	    }
+	    bits <<= 1;
+
+	    median(graph, xvecs[i], nvtxs, active, temp_goal, using_vwgts, temp_sets);
+	    for (j = 1; j <= nvtxs; j++)
+		sets[j] = (sets[j] << 1) + temp_sets[j];
+	}
+	sfree((char *) temp_sets);
+    }
+
+    if (mediantype == 4) {	/* Stripe the domain. */
+	rec_median_1(graph, xvecs[1], nvtxs, active, cube_or_mesh, nsets,
+		     goal, using_vwgts, sets, TRUE);
+    }
+}

contrib/Chaco/assign/median.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"structs.h"
+#include	"defs.h"
+
+
+/* Find the median of set of values. */
+/* Can also find nested medians of several sets of values */
+/* Routine works by repeatedly guessing a value, and discarding those */
+/* values which are on the wrong side of the guess. */
+
+void      median(graph, vals, nvtxs, active, goal, using_vwgts, sets)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double   *vals;			/* values of which to find median */
+int       nvtxs;		/* number of values I own */
+int      *active;		/* space for list of nvtxs ints */
+double   *goal;			/* desired sizes for sets */
+int       using_vwgts;		/* are vertex weights being used? */
+short    *sets;			/* set each vertex gets assigned to */
+{
+    double   *vptr;		/* loops through vals array */
+    double    val;		/* value in vals array */
+    double    maxval;		/* largest active value */
+    double    minval;		/* smallest active value */
+    double    guess;		/* approximate median value */
+    double    nearup;		/* lowest guy above guess */
+    double    neardown;		/* highest guy below guess */
+    double    whigh;		/* total weight of values above maxval */
+    double    wlow;		/* total weight of values below minval */
+    double    wabove;		/* total weight of active values above guess */
+    double    wbelow;		/* total weight of active values below guess */
+    double    wexact;		/* weight of vertices exactly at guess */
+    double    lweight;		/* desired weight of lower values in set */
+    double    uweight;		/* desired weight of upper values in set */
+    double    frac;		/* fraction of values I want less than guess */
+    int      *aptr;		/* loops through active array */
+    int      *aptr2;		/* helps update active array */
+    int       myactive;		/* number of active values I own */
+    double    wfree;		/* weight of vtxs not yet divided */
+    int       removed;		/* number of my values eliminated */
+    /*int npass = 0;*/		/* counts passes required to find median */
+    int       done;		/* check for termination criteria */
+    int       vtx;		/* vertex being considered */
+    int       i;		/* loop counters */
+    void      median_assign();
+
+    /* Initialize. */
+
+    /* Determine the desired weight sums for the two different sets. */
+    lweight = goal[0];
+    uweight = goal[1];
+
+    myactive = nvtxs;
+    whigh = wlow = 0;
+
+    /* Find largest and smallest values in vector, and construct active list. */
+    vptr = vals;
+    aptr = active;
+    minval = maxval = *(++vptr);
+    *aptr++ = 1;
+    for (i = 2; i <= nvtxs; i++) {
+	*aptr++ = i;
+	val = *(++vptr);
+	if (val > maxval)
+	    maxval = val;
+	if (val < minval)
+	    minval = val;
+    }
+
+    /* Loop until all sets are partitioned correctly. */
+    done = FALSE;
+    while (!done) {
+	/*npass++;*/
+
+	/* Select a potential dividing value. */
+	/* Currently, this assumes a linear distribution. */
+	wfree = lweight + uweight - (wlow + whigh);
+	frac = (lweight - wlow) / wfree;
+
+	/* Overshoot a bit to try to cut into largest set. */
+	frac = .5 * (frac + .5);
+
+	guess = minval + frac * (maxval - minval);
+
+	/* Now count the guys above and below this guess. */
+	/* Also find nearest values on either side of guess. */
+	wabove = wbelow = wexact = 0;
+	nearup = maxval;
+	neardown = minval;
+
+	aptr = active;
+	for (i = 0; i < myactive; i++) {
+	    vtx = *aptr++;
+	    val = vals[vtx];
+	    if (val > guess) {
+		if (using_vwgts)
+		    wabove += graph[vtx]->vwgt;
+		else
+		    wabove++;
+		if (val < nearup)
+		    nearup = val;
+	    }
+	    else if (val < guess) {
+		if (using_vwgts)
+		    wbelow += graph[vtx]->vwgt;
+		else
+		    wbelow++;
+		if (val > neardown)
+		    neardown = val;
+	    }
+	    else {
+		if (using_vwgts)
+		    wexact += graph[vtx]->vwgt;
+		else
+		    wexact++;
+	    }
+	}
+
+	/* Select a half to discard. */
+	/* And remove discarded vertices from active list. */
+	removed = 0;
+/*	if (wlow + wbelow - lweight > whigh + wabove - uweight) {*/
+	if (wlow + wbelow - lweight > whigh + wabove - uweight &&
+	    whigh + wabove + wexact < uweight ) {
+	    /* Discard upper set. */
+	    whigh += wabove + wexact;
+	    maxval = neardown;
+	    done = FALSE;
+	    aptr = aptr2 = active;
+	    for (i = 0; i < myactive; i++) {
+		if (vals[*aptr] >= guess) {
+		    ++removed;
+		    if (vals[*aptr] == guess) 
+		        nearup = guess;
+		}
+		else
+		    *aptr2++ = *aptr;
+		aptr++;
+	    }
+	    myactive -= removed;
+	    if (myactive == 0) done = TRUE;
+	}
+/*	else if ( whigh + wabove - uweight > wlow + wbelow - lweight) {*/
+	else if ( whigh + wabove - uweight > wlow + wbelow - lweight &&
+	    wlow + wbelow + wexact < lweight ) {
+	    /* Discard lower set. */
+	    wlow += wbelow + wexact;
+	    minval = nearup;
+	    done = FALSE;
+	    aptr = aptr2 = active;
+	    for (i = 0; i < myactive; i++) {
+		if (vals[*aptr] <= guess) {
+		    ++removed;
+		    if (vals[*aptr] == guess) 
+		        neardown = guess;
+		}
+		else
+		    *aptr2++ = *aptr;
+		aptr++;
+	    }
+	    myactive -= removed;
+	    if (myactive == 0) done = TRUE;
+	}
+	else {			/* Perfect partition! */
+	    wlow += wbelow;
+	    whigh += wabove;
+	    /*
+	    minval = nearup;
+	    maxval = neardown;
+	    myactive = 0;
+	    */
+	    done = TRUE;
+	}
+
+	/* Check for alternate termination criteria. */
+	if (!done && maxval == minval) {
+	    guess = maxval;
+	    done = TRUE;
+	}
+    }
+    median_assign(graph, vals, nvtxs, goal, using_vwgts, sets, wlow, whigh, guess);
+}
+
+
+void      median_assign(graph, vals, nvtxs, goal, using_vwgts, sets,
+			          wlow, whigh, guess)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double   *vals;			/* values of which to find median */
+int       nvtxs;		/* number of values I own */
+double   *goal;			/* desired sizes for sets */
+int       using_vwgts;		/* are vertex weights being used? */
+short    *sets;			/* assigned set for each vertex */
+double    wlow;			/* sum of weights below guess */
+double    whigh;		/* sum of weights above guess */
+double    guess;		/* median value */
+{
+    int       i;		/* loop counter */
+
+    for (i = 1; i <= nvtxs; i++) {
+	if (vals[i] < guess)
+	    sets[i] = 0;
+	else if (vals[i] > guess)
+	    sets[i] = 1;
+	else {
+	    if (goal[0] - wlow > goal[1] - whigh) {
+		sets[i] = 0;
+		if (using_vwgts)
+		    wlow += graph[i]->vwgt;
+		else
+		    wlow++;
+	    }
+	    else {
+		sets[i] = 1;
+		if (using_vwgts)
+		    whigh += graph[i]->vwgt;
+		else
+		    whigh++;
+	    }
+	}
+    }
+}

contrib/Chaco/assign/merge_assign.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+/* Combine the old assignment value with the new partition. */
+void      merge_assignments(assignment, subassign, subsets, subnvtxs, loc2glob)
+short    *assignment;		/* assignment list for graph */
+short    *subassign;		/* subgraph assignment list */
+short    *subsets;		/* mapping from local to global sets */
+int       subnvtxs;		/* number of vtxs in subgraph */
+int      *loc2glob;		/* subgraph -> graph numbering map */
+{
+    int       i;		/* loop counter */
+
+    for (i = 1; i <= subnvtxs; i++) {
+	assignment[loc2glob[i]] = subsets[subassign[i]];
+    }
+}

contrib/Chaco/assign/rec_median.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"structs.h"
+#include	"params.h"
+#include	"defs.h"
+
+/* Recursively apply median to a SINGLE vector of values */
+
+void      rec_median_1(graph, vals, nvtxs, active, cube_or_mesh, nsets,
+		                 goal, using_vwgts, assign, top)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double   *vals;			/* values of which to find median */
+int       nvtxs;		/* number of values I own */
+int      *active;		/* space for list of nmyvals ints */
+int       cube_or_mesh;		/* 0=> hypercube, other=> mesh */
+int       nsets;		/* number of sets to divide into */
+double   *goal;			/* desired sizes for sets */
+int       using_vwgts;		/* are vertex weights being used? */
+short    *assign;		/* set each vertex gets assigned to */
+int       top;			/* is this the top call in the recursion? */
+{
+    struct vtx_data **sub_graph;/* subgraph data structure with vertex weights */
+    double   *sub_vals;		/* subgraph entries in vals vector */
+    double    merged_goal[MAXSETS / 2];	/* combined goal values */
+    double    sub_vwgt_sum;	/* sum of vertex weights in subgraph */
+    int      *loc2glob;		/* mapping from subgraph to graph numbering */
+    int       mapping[MAXSETS];	/* appropriate set number */
+    short    *sub_assign;	/* assignment returned from subgraph */
+    int       sub_nvtxs;	/* number of vertices in subgraph */
+    int       sub_nsets;	/* number of sets in side of first partition */
+    int       setsize[2];	/* number of vertices in two subsets */
+    int       maxsize;		/* number of vertices in larger subset */
+    int       ndims;		/* number of bits defining set */
+    int       mesh_dims[3];	/* shape of fictitious mesh */
+    int       i, j;		/* loop counters */
+    double   *smalloc();
+    int       bit_reverse(), sfree();
+    void      median(), make_subgoal(), make_subvector();
+    void      make_maps2();
+
+    cube_or_mesh = (cube_or_mesh != 0);
+    mesh_dims[1] = mesh_dims[2] = 1;
+
+    if (!cube_or_mesh) {
+	for (i = 0; i < 2; i++) {
+	    merged_goal[i] = 0;
+	    for (j = i; j < nsets; j += 2)
+		merged_goal[i] += goal[j];
+	}
+    }
+    else {
+	merged_goal[0] = merged_goal[1] = 0;
+	for (i = 0; i < (nsets + 1) / 2; i++)
+	    merged_goal[0] += goal[i];
+	for (i = (nsets + 1) / 2; i < nsets; i++)
+	    merged_goal[1] += goal[i];
+    }
+
+    median(graph, vals, nvtxs, active, merged_goal, using_vwgts, assign);
+
+    if (nsets > 2) {
+	/* Find size of largest subproblem. */
+	setsize[0] = setsize[1] = 0;
+	for (i = 1; i <= nvtxs; i++)
+	    ++setsize[assign[i]];
+	maxsize = max(setsize[0], setsize[1]);
+
+	sub_assign = (short *) smalloc((unsigned) (maxsize + 1) * sizeof(short));
+	sub_vals = (double *) smalloc((unsigned) (maxsize + 1) * sizeof(double));
+	loc2glob = (int *) smalloc((unsigned) (maxsize + 1) * sizeof(int));
+	if (!using_vwgts)
+	    sub_graph = NULL;
+	else
+	    sub_graph = (struct vtx_data **)
+	       smalloc((unsigned) (maxsize + 1) * sizeof(struct vtx_data *));
+
+	for (i = 0; i < 2; i++) {
+	    /* Construct subproblem. */
+	    if (i == 0) {
+		sub_nsets = (nsets + 1) / 2;
+	    }
+	    else {
+		sub_nsets = nsets / 2;
+	    }
+	    sub_nvtxs = setsize[i];
+
+	    for (j = 1; j <= sub_nvtxs; j++)
+		sub_assign[j] = 0;
+	    make_maps2(assign, nvtxs, i, (int *) NULL, loc2glob);
+
+	    if (sub_nsets > 1) {
+
+		if (!using_vwgts) {
+		    sub_vwgt_sum = sub_nvtxs;
+		}
+		else {
+		    sub_vwgt_sum = 0;
+		    for (j = 1; j <= sub_nvtxs; j++) {
+			sub_graph[j] = graph[loc2glob[j]];
+			sub_vwgt_sum += sub_graph[j]->vwgt;
+		    }
+		}
+
+		make_subvector(vals, sub_vals, sub_nvtxs, loc2glob);
+
+		mesh_dims[0] = nsets;
+		make_subgoal(goal, merged_goal, 2, cube_or_mesh, nsets, mesh_dims,
+			     i, sub_vwgt_sum);
+
+		rec_median_1(sub_graph, sub_vals, sub_nvtxs, active, cube_or_mesh,
+			     sub_nsets, merged_goal, using_vwgts,
+			     sub_assign, FALSE);
+	    }
+	    /* Merge subassignment with current assignment. */
+	    for (j = 1; j <= sub_nvtxs; j++) {
+		assign[loc2glob[j]] |= (sub_assign[j] << 1);
+	    }
+	}
+
+	/* Now I have the set striped in bit reversed order. */
+	if (top) {
+	    ndims = 0;
+	    for (i = 1; i < nsets; i <<= 1)
+		ndims++;
+
+	    for (i = 0; i < nsets; i++)
+		mapping[i] = bit_reverse(i, ndims);
+
+	    for (i = 1; i <= nvtxs; i++)
+		assign[i] = mapping[assign[i]];
+	}
+
+	if (sub_graph != NULL)
+	    sfree((char *) sub_graph);
+	sfree((char *) loc2glob);
+	sfree((char *) sub_vals);
+	sfree((char *) sub_assign);
+    }
+}
+
+
+
+/* Recursively apply median to a SEVERAL vectors of values */
+/* Divide with first, and use result to divide with second, etc. */
+/* Note: currently only works for power-of-two number of processors. */
+
+void      rec_median_k(graph, vals, nvtxs, active, ndims, cube_or_mesh,
+		                 goal, using_vwgts, assign)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double  **vals;			/* values of which to find median */
+int       nvtxs;		/* number of values I own */
+int      *active;		/* space for list of nmyvals ints */
+int       ndims;		/* number of dimensions to divide */
+int       cube_or_mesh;		/* 0 => hypercube, d => d-dimensional mesh */
+double   *goal;			/* desired sizes for sets */
+int       using_vwgts;		/* are vertex weights being used? */
+short    *assign;		/* set each vertex gets assigned to */
+{
+    struct vtx_data **sub_graph;/* subgraph data structure with vertex weights */
+    double   *sub_vals[MAXDIMS];/* subgraph entries in vals vectors */
+    double    merged_goal[MAXSETS / 2];	/* combined goal values */
+    double    sub_vwgt_sum;	/* sum of vertex weights in subgraph */
+    int      *loc2glob;		/* mapping from subgraph to graph numbering */
+    short    *sub_assign;	/* assignment returned from subgraph */
+    int       sub_nvtxs;	/* number of vertices in subgraph */
+    int       setsize[2];	/* number of vertices in two subsets */
+    int       maxsize;		/* number of vertices in larger subset */
+    int       nsets;		/* number of sets we are dividing into */
+    int       mesh_dims[3];	/* shape of fictitious mesh */
+    int       i, j;		/* loop counters */
+    double   *smalloc();
+    int       sfree();
+    void      median(), make_subgoal(), make_subvector(), make_maps2();
+
+    mesh_dims[1] = mesh_dims[2] = 1;
+
+    /* Note: This is HYPERCUBE/MESH dependent.  We'll want to combine the */
+    /* sizes of different sets on the different architectures. */
+    nsets = 1 << ndims;
+    for (i = 0; i < 2; i++) {
+	merged_goal[i] = 0;
+	for (j = i; j < nsets; j += 2)
+	    merged_goal[i] += goal[j];
+    }
+
+    median(graph, vals[1], nvtxs, active, merged_goal, using_vwgts, assign);
+
+    if (ndims > 1) {
+	/* Find size of largest subproblem. */
+	setsize[0] = setsize[1] = 0;
+	for (i = 1; i <= nvtxs; i++)
+	    ++setsize[assign[i]];
+	maxsize = max(setsize[0], setsize[1]);
+
+	sub_assign = (short *) smalloc((unsigned) (maxsize + 1) * sizeof(short));
+	for (i = 1; i < ndims; i++) {
+	    sub_vals[i] = (double *) smalloc((unsigned) (maxsize + 1) * sizeof(double));
+	}
+	loc2glob = (int *) smalloc((unsigned) (maxsize + 1) * sizeof(int));
+	if (!using_vwgts)
+	    sub_graph = NULL;
+	else
+	    sub_graph = (struct vtx_data **)
+	       smalloc((unsigned) (maxsize + 1) * sizeof(struct vtx_data *));
+
+	for (i = 0; i < 2; i++) {
+	    /* Construct subproblem. */
+
+	    sub_nvtxs = setsize[i];
+
+	    for (j = 1; j <= sub_nvtxs; j++)
+		sub_assign[j] = 0;
+
+	    make_maps2(assign, nvtxs, i, (int *) NULL, loc2glob);
+
+	    if (!using_vwgts)
+		sub_vwgt_sum = sub_nvtxs;
+	    else {
+		sub_vwgt_sum = 0;
+		for (j = 1; j <= sub_nvtxs; j++) {
+		    sub_graph[j] = graph[loc2glob[j]];
+		    sub_vwgt_sum += sub_graph[j]->vwgt;
+		}
+	    }
+
+	    for (j = 2; j <= ndims; j++) {
+		make_subvector(vals[j], sub_vals[j - 1], sub_nvtxs, loc2glob);
+	    }
+
+	    mesh_dims[0] = nsets;
+	    make_subgoal(goal, merged_goal, 2, cube_or_mesh, nsets, mesh_dims,
+			 i, sub_vwgt_sum);
+
+	    rec_median_k(sub_graph, sub_vals, sub_nvtxs, active, ndims - 1,
+			 cube_or_mesh, merged_goal, using_vwgts,
+			 sub_assign);
+
+	    /* Merge subassignment with current assignment. */
+	    for (j = 1; j <= sub_nvtxs; j++) {
+		assign[loc2glob[j]] |= (sub_assign[j] << 1);
+	    }
+	}
+
+	if (sub_graph != NULL)
+	    sfree((char *) sub_graph);
+	sfree((char *) loc2glob);
+	for (i = 1; i < ndims; i++)
+	    sfree((char *) sub_vals[i]);
+	sfree((char *) sub_assign);
+    }
+}

contrib/Chaco/assign/rotate.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	<math.h>
+#include	"defs.h"
+
+
+void      rotate2d(yvecs, nmyvtxs, theta)
+double  **yvecs;		/* ptr to list of y-vectors (rotated) */
+int       nmyvtxs;		/* length of yvecs */
+double    theta;		/* angle to rotate by */
+{
+    double    temp1;		/* hold values for a while */
+    double    c, s;		/* cosine and sine of theta */
+    int       i;		/* loop counter */
+
+    s = sin(theta);
+    c = cos(theta);
+
+    for (i = 1; i <= nmyvtxs; i++) {
+	temp1 = yvecs[1][i];
+	yvecs[1][i] = c * temp1 + s * yvecs[2][i];
+	yvecs[2][i] = -s * temp1 + c * yvecs[2][i];
+    }
+}
+
+
+void      rotate3d(yvecs, nmyvtxs, theta, phi, gamma2)
+double  **yvecs;		/* ptr to list of y-vectors (to be rotated) */
+int       nmyvtxs;		/* length of yvecs */
+double    theta, phi, gamma2;	/* rotational parameters */
+{
+    double    temp1, temp2;	/* hold values for a while */
+    double    ctheta, stheta;	/* cosine and sine of theta */
+    double    cphi, sphi;	/* cosine and sine of phi */
+    double    cgamma, sgamma;	/* cosine and sine of gamma */
+    double    onemcg;		/* 1.0 - cosine(gamma) */
+    double    a1, a2, a3;	/* rotation matrix entries */
+    double    b1, b2, b3;	/* rotation matrix entries */
+    double    c1, c2, c3;	/* rotation matrix entries */
+    int       i;		/* loop counter */
+
+    stheta = sin(theta);
+    ctheta = cos(theta);
+    sphi = sin(phi);
+    cphi = cos(phi);
+    sgamma = sin(gamma2);
+    cgamma = cos(gamma2);
+
+    onemcg = 1.0 - cgamma;
+
+    a1 = cgamma + cphi * ctheta * onemcg * cphi * ctheta;
+    a2 = sgamma * sphi + cphi * stheta * onemcg * cphi * ctheta;
+    a3 = -sgamma * cphi * stheta + sphi * onemcg * cphi * ctheta;
+
+    b1 = -sgamma * sphi + cphi * ctheta * onemcg * cphi * stheta;
+    b2 = cgamma + cphi * stheta * onemcg * cphi * stheta;
+    b3 = sgamma * cphi * ctheta + sphi * onemcg * cphi * stheta;
+
+    c1 = sgamma * cphi * stheta + cphi * ctheta * onemcg * sphi;
+    c2 = -sgamma * cphi * ctheta + cphi * stheta * onemcg * sphi;
+    c3 = cgamma + sphi * onemcg * sphi;
+
+    for (i = 1; i <= nmyvtxs; i++) {
+	temp1 = yvecs[1][i];
+	temp2 = yvecs[2][i];
+
+	yvecs[1][i] = a1 * temp1 + b1 * temp2 + c1 * yvecs[3][i];
+	yvecs[2][i] = a2 * temp1 + b2 * temp2 + c2 * yvecs[3][i];
+	yvecs[3][i] = a3 * temp1 + b3 * temp2 + c3 * yvecs[3][i];
+    }
+}

contrib/Chaco/assign/y2x.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+
+
+void      y2x(xvecs, ndims, nmyvtxs, wsqrt)
+double  **xvecs;		/* pointer to list of x-vectors */
+int       ndims;		/* number of divisions to make (# xvecs) */
+int       nmyvtxs;		/* number of vertices I own (lenght xvecs) */
+double   *wsqrt;		/* sqrt of vertex weights */
+
+/* Convert from y to x by dividing by wsqrt. */
+{
+    double   *wptr;		/* loops through wsqrt */
+    double   *xptr;		/* loops through elements of a xvec */
+    int       i, j;		/* loop counters */
+
+    if (wsqrt == NULL)
+	return;
+
+    for (i = 1; i <= ndims; i++) {
+	xptr = xvecs[i];
+	wptr = wsqrt;
+	for (j = nmyvtxs; j; j--) {
+	    *(++xptr) /= *(++wptr);
+	}
+    }
+}
+
+
+void      x2y(yvecs, ndims, nmyvtxs, wsqrt)
+double  **yvecs;		/* pointer to list of y-vectors */
+int       ndims;		/* number of divisions to make (# yvecs) */
+int       nmyvtxs;		/* number of vertices I own (lenght yvecs) */
+double   *wsqrt;		/* sqrt of vertex weights */
+
+/* Convert from x to y by multiplying by wsqrt. */
+{
+    double   *wptr;		/* loops through wsqrt */
+    double   *yptr;		/* loops through elements of a yvec */
+    int       i, j;		/* loop counters */
+
+    if (wsqrt == NULL)
+	return;
+
+    for (i = 1; i <= ndims; i++) {
+	yptr = yvecs[i];
+	wptr = wsqrt;
+	for (j = nmyvtxs; j; j--) {
+	    *(++yptr) *= *(++wptr);
+	}
+    }
+}

contrib/Chaco/bpmatch/checkbp.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	<math.h>
+#include	"structs.h"
+#include	"params.h"
+#include	"defs.h"
+
+
+/* Confirm that the bipartite match algorithm did the right thing. */
+void      checkbp(graph, xvecs, sets, dists, nvtxs, ndims)
+struct vtx_data **graph;	/* graph data structure for vertex weights */
+double  **xvecs;		/* values to partition */
+short    *sets;			/* set assignments returned */
+double   *dists;		/* distances that separate sets */
+int       nvtxs;		/* number of vertices */
+int       ndims;		/* number of dimensions for division */
+{
+    int       signs[MAXDIMS];	/* signs for coordinates of target points */
+    int       sizes[MAXSETS];	/* size of each set */
+    int       weights[MAXSETS];	/* size of each set */
+    double    setval;		/* value from assigned set */
+    double    val, bestval;	/* value to decide set assignment */
+    double    tol = 1.0e-8;	/* numerical tolerence */
+    int       error = FALSE;	/* are errors encountered? */
+    int       nsets;		/* number of sets */
+    int       bestset;		/* set vtx should be assigned to */
+    int       i, j, k;		/* loop counter */
+    void      checkpnt();
+
+    nsets = 1 << ndims;
+
+    for (i = 0; i < nsets; i++) {
+	sizes[i] = 0;
+	weights[i] = 0;
+    }
+
+    for (i = 1; i <= nvtxs; i++) {
+	/* Is vertex closest to the set it is assigned to? */
+	for (j = 0; j < ndims; j++)
+	    signs[j] = -1;
+	bestval = 0;
+	for (j = 0; j < nsets; j++) {
+	    val = -dists[j];
+	    for (k = 1; k <= ndims; k++) {
+		val += 2 * signs[k - 1] * xvecs[k][i];
+	    }
+	    if (j == sets[i])
+		setval = val;
+	    if (j == 0 || val < bestval) {
+		bestval = val;
+		bestset = j;
+	    }
+	    if (signs[0] == 1 && signs[1] == 1)
+		signs[2] *= -1;
+	    if (signs[0] == 1)
+		signs[1] *= -1;
+	    signs[0] *= -1;
+	}
+	if (fabs(setval - bestval) >= tol * (fabs(setval) + fabs(bestval))) {
+	    printf(" Vtx %d in set %d (%e), but should be in %d (%e)\n",
+		   i, sets[i], setval, bestset, bestval);
+	    error = TRUE;
+	}
+	++sizes[sets[i]];
+	weights[sets[i]] += graph[i]->vwgt;
+    }
+
+    printf(" Sizes:");
+    for (i = 0; i < nsets; i++)
+	printf(" %d(%d)", sizes[i], weights[i]);
+    printf("\n");
+
+    if (error)
+	checkpnt("ERROR in checkbp");
+}

contrib/Chaco/bpmatch/genvals2d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+
+
+void      genvals2d(xvecs, vals, nvtxs)
+/* Create lists of sets of values to be sorted. */
+double  **xvecs;		/* vectors to partition */
+double   *vals[4][MAXSETS];	/* ptrs to lists of values */
+int       nvtxs;		/* number of values */
+{
+    int       nlists = 4;	/* number of lists to generate */
+    double   *temp[4];		/* place holders for vals */
+    int       i;		/* loop counter */
+    double   *smalloc();
+
+    for (i = 0; i < nlists; i++) {
+	temp[i] = (double *) smalloc((unsigned) nvtxs * sizeof(double));
+    }
+
+    for (i = 1; i <= nvtxs; i++) {
+	temp[0][i - 1] = 4 * xvecs[1][i];
+	temp[1][i - 1] = 4 * xvecs[2][i];
+	temp[2][i - 1] = 4 * (xvecs[1][i] + xvecs[2][i]);
+	temp[3][i - 1] = 4 * (xvecs[2][i] - xvecs[1][i]);
+    }
+
+    vals[0][1] = vals[1][0] = vals[2][3] = vals[3][2] = temp[0];
+    vals[0][2] = vals[2][0] = vals[1][3] = vals[3][1] = temp[1];
+    vals[0][3] = vals[3][0] = temp[2];
+    vals[1][2] = vals[2][1] = temp[3];
+}

contrib/Chaco/bpmatch/genvals3d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+
+
+void      genvals3d(xvecs, vals, nvtxs)
+/* Create lists of sets of values to be sorted. */
+double  **xvecs;		/* vectors to partition */
+double   *vals[8][MAXSETS];	/* ptrs to lists of values */
+int       nvtxs;		/* number of values */
+{
+    int       nsets = 8;	/* number of sets */
+    int       nlists = 13;	/* number of lists to generate */
+    double   *temp[13];		/* place holders for vals */
+    int       i, j;		/* loop counter */
+    double   *smalloc();
+
+    for (i = 0; i < nlists; i++) {
+	temp[i] = (double *) smalloc((unsigned) nvtxs * sizeof(double));
+    }
+
+    for (i = 1; i <= nvtxs; i++) {
+	temp[0][i - 1] = 4 * xvecs[1][i];
+	temp[1][i - 1] = 4 * xvecs[2][i];
+	temp[2][i - 1] = 4 * xvecs[3][i];
+	temp[3][i - 1] = 4 * (xvecs[1][i] + xvecs[2][i]);
+	temp[4][i - 1] = 4 * (-xvecs[1][i] + xvecs[2][i]);
+	temp[5][i - 1] = 4 * (xvecs[1][i] + xvecs[3][i]);
+	temp[6][i - 1] = 4 * (-xvecs[1][i] + xvecs[3][i]);
+	temp[7][i - 1] = 4 * (xvecs[2][i] + xvecs[3][i]);
+	temp[8][i - 1] = 4 * (-xvecs[2][i] + xvecs[3][i]);
+	temp[9][i - 1] = 4 * (xvecs[1][i] + xvecs[2][i] + xvecs[3][i]);
+	temp[10][i - 1] = 4 * (-xvecs[1][i] + xvecs[2][i] + xvecs[3][i]);
+	temp[11][i - 1] = 4 * (xvecs[1][i] - xvecs[2][i] + xvecs[3][i]);
+	temp[12][i - 1] = 4 * (-xvecs[1][i] - xvecs[2][i] + xvecs[3][i]);
+    }
+
+    vals[0][1] = vals[2][3] = vals[4][5] = vals[6][7] = temp[0];
+    vals[0][2] = vals[1][3] = vals[4][6] = vals[5][7] = temp[1];
+    vals[0][4] = vals[1][5] = vals[2][6] = vals[3][7] = temp[2];
+    vals[0][3] = vals[4][7] = temp[3];
+    vals[1][2] = vals[5][6] = temp[4];
+    vals[0][5] = vals[2][7] = temp[5];
+    vals[1][4] = vals[3][6] = temp[6];
+    vals[0][6] = vals[1][7] = temp[7];
+    vals[2][4] = vals[3][5] = temp[8];
+    vals[0][7] = temp[9];
+    vals[1][6] = temp[10];
+    vals[2][5] = temp[11];
+    vals[3][4] = temp[12];
+
+    for (i = 0; i < nsets; i++) {
+	for (j = i + 1; j < nsets; j++)
+	    vals[j][i] = vals[i][j];
+    }
+}

contrib/Chaco/bpmatch/inits2d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"params.h"
+#include	"structs.h"
+#include	"defs.h"
+
+
+void      inits2d(graph, xvecs, vals, indices, nvtxs, dist, startvtx, size, sets)
+struct vtx_data **graph;	/* graph data structure for vertex weights */
+double  **xvecs;		/* values to partition with */
+double   *vals[4][MAXSETS];	/* values in sorted lists */
+int      *indices[4][MAXSETS];	/* indices sorting lists */
+int       nvtxs;		/* number of vertices */
+double   *dist;			/* trial separation point */
+int       startvtx[4][MAXSETS];	/* indices defining separation */
+double   *size;			/* size of each set being modified */
+short    *sets;			/* set each vertex gets assigned to */
+{
+    double    xmid, ymid;	/* median x and y values */
+    double    val, bestval;	/* values for determining set preferences */
+    short     bestset;		/* set vertex wants to be in */
+    int       signx, signy;	/* sign values for different target points */
+    int       nsets = 4;	/* number of different sets */
+    int       i, j;		/* loop counters */
+    int       findindex();
+
+/*
+    xmid = .25 * (vals[0][1][indices[0][1][nvtxs / 2]] +
+		  vals[0][1][indices[0][1][nvtxs / 2 - 1]]);
+    ymid = .25 * (vals[0][2][indices[0][2][nvtxs / 2]] +
+		  vals[0][2][indices[0][2][nvtxs / 2 - 1]]);
+*/
+    xmid = .5 * vals[0][1][indices[0][1][nvtxs / 2]];
+    ymid = .5 * vals[0][2][indices[0][2][nvtxs / 2]];
+
+    dist[0] = -xmid - ymid;
+    dist[1] = xmid - ymid;
+    dist[2] = -xmid + ymid;
+    dist[3] = xmid + ymid;
+
+    /* Now initialize startvtx. */
+    startvtx[0][1] = startvtx[2][3] = nvtxs / 2;
+    startvtx[0][2] = startvtx[1][3] = nvtxs / 2;
+    startvtx[1][2] = findindex(indices[1][2], vals[1][2], dist[2] - dist[1], nvtxs);
+    startvtx[0][3] = findindex(indices[0][3], vals[0][3], dist[3] - dist[0], nvtxs);
+
+    for (i = 0; i < nsets; i++)
+	size[i] = 0;
+
+    for (i = 1; i <= nvtxs; i++) {
+	/* Which set is this vertex in? */
+	signx = signy = -1;
+	bestval = 0;
+	for (j = 0; j < nsets; j++) {
+	    val = -dist[j] + 2 * (signx * xvecs[1][i] + signy * xvecs[2][i]);
+	    if (j == 0 || val < bestval) {
+		bestval = val;
+		bestset = (short) j;
+	    }
+	    if (signx == 1)
+		signy *= -1;
+	    signx *= -1;
+	}
+	sets[i] = bestset;
+	size[bestset] += graph[i]->vwgt;
+    }
+}
+
+
+int       findindex(indices, vals, target, nvals)
+int      *indices;		/* indices sorting values */
+double   *vals;			/* values sorted by indices */
+double    target;		/* target value */
+int       nvals;		/* number of values */
+{
+    double    ratio;		/* interpolation parameter */
+    double    vlow, vhigh;	/* values at limits of search range */
+    int       low, high;	/* range left to search */
+    int       new;		/* new index limit */
+
+    if (target <= vals[indices[0]])
+	return (0);
+    if (target >= vals[indices[nvals - 1]])
+	return (nvals - 1);
+
+    low = 0;
+    high = nvals - 1;
+
+    while (high - low > 1) {
+	vlow = vals[indices[low]];
+	vhigh = vals[indices[high]];
+	if (vlow == vhigh)
+	    return ((vlow + vhigh) / 2);
+
+	ratio = (target - vlow) / (vhigh - vlow);
+	new = low + ratio * (high - low);
+	if (new == low)
+	    ++new;
+	else if (new == high)
+	    --new;
+
+	if (vals[indices[new]] < target)
+	    low = new;
+	else
+	    high = new;
+    }
+
+    if (target == vals[indices[high]])
+	return (high);
+    else
+	return (low);
+}

contrib/Chaco/bpmatch/inits3d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"params.h"
+#include	"structs.h"
+#include	"defs.h"
+
+
+void      inits3d(graph, xvecs, vals, indices, nvtxs, dist, startvtx, size, sets)
+struct vtx_data **graph;	/* graph data structure for vertex weights */
+double  **xvecs;		/* values to partition with */
+double   *vals[8][MAXSETS];	/* values in sorted lists */
+int      *indices[8][MAXSETS];	/* indices sorting lists */
+int       nvtxs;		/* number of vertices */
+double   *dist;			/* trial separation point */
+int       startvtx[8][MAXSETS];	/* indices defining separation */
+double   *size;			/* size of each set being modified */
+short    *sets;			/* set each vertex gets assigned to */
+{
+    double    xmid, ymid, zmid;	/* median x, y and z values */
+    double    val, bestval;	/* values for determining set preferences */
+    short     bestset;		/* set vertex wants to be in */
+    int       signx, signy, signz;	/* sign values for different target points */
+    int       nsets = 8;	/* number of different sets */
+    int       i, j;		/* loop counters */
+    int       findindex();
+
+/*
+    xmid = .25 * (vals[0][1][indices[0][1][nvtxs / 2]] +
+		  vals[0][1][indices[0][1][nvtxs / 2 - 1]]);
+    ymid = .25 * (vals[0][2][indices[0][2][nvtxs / 2]] +
+		  vals[0][2][indices[0][2][nvtxs / 2 - 1]]);
+    zmid = .25 * (vals[0][4][indices[0][4][nvtxs / 2]] +
+		  vals[0][4][indices[0][4][nvtxs / 2 - 1]]);
+*/
+    xmid = .5 * vals[0][1][indices[0][1][nvtxs / 2]];
+    ymid = .5 * vals[0][2][indices[0][2][nvtxs / 2]];
+    zmid = .5 * vals[0][4][indices[0][4][nvtxs / 2]];
+
+    dist[0] = -xmid - ymid - zmid;
+    dist[1] = xmid - ymid - zmid;
+    dist[2] = -xmid + ymid - zmid;
+    dist[3] = xmid + ymid - zmid;
+    dist[4] = -xmid - ymid + zmid;
+    dist[5] = xmid - ymid + zmid;
+    dist[6] = -xmid + ymid + zmid;
+    dist[7] = xmid + ymid + zmid;
+
+    /* Now initialize startvtxs. */
+    startvtx[0][1] = startvtx[2][3] = startvtx[4][5] = startvtx[6][7] = nvtxs / 2;
+    startvtx[0][2] = startvtx[1][3] = startvtx[4][6] = startvtx[5][7] = nvtxs / 2;
+    startvtx[0][4] = startvtx[1][5] = startvtx[2][6] = startvtx[3][7] = nvtxs / 2;
+    startvtx[0][3] = startvtx[4][7] =
+       findindex(indices[0][3], vals[0][3], dist[3] - dist[0], nvtxs);
+    startvtx[1][2] = startvtx[5][6] =
+       findindex(indices[1][2], vals[1][2], dist[2] - dist[1], nvtxs);
+    startvtx[0][5] = startvtx[2][7] =
+       findindex(indices[0][5], vals[0][5], dist[5] - dist[0], nvtxs);
+    startvtx[1][4] = startvtx[3][6] =
+       findindex(indices[1][4], vals[1][4], dist[4] - dist[1], nvtxs);
+    startvtx[0][6] = startvtx[1][7] =
+       findindex(indices[0][6], vals[0][6], dist[6] - dist[0], nvtxs);
+    startvtx[2][4] = startvtx[3][5] =
+       findindex(indices[2][4], vals[2][4], dist[4] - dist[2], nvtxs);
+    startvtx[0][7] = findindex(indices[0][7], vals[0][7], dist[7] - dist[0], nvtxs);
+    startvtx[1][6] = findindex(indices[1][6], vals[1][6], dist[6] - dist[1], nvtxs);
+    startvtx[2][5] = findindex(indices[2][5], vals[2][5], dist[5] - dist[2], nvtxs);
+    startvtx[3][4] = findindex(indices[3][4], vals[3][4], dist[4] - dist[3], nvtxs);
+
+    /* Finally, determine the set sizes based on this splitter. */
+
+    for (i = 0; i < nsets; i++)
+	size[i] = 0;
+
+    for (i = 1; i <= nvtxs; i++) {
+	/* Which set is this vertex in? */
+	signx = signy = signz = -1;
+	bestval = 0;
+	for (j = 0; j < nsets; j++) {
+	    val = -dist[j] + 2 * (signx * xvecs[1][i] + signy * xvecs[2][i] + signz * xvecs[3][i]);
+	    if (j == 0 || val < bestval) {
+		bestval = val;
+		bestset = (short) j;
+	    }
+	    if (signx == 1 && signy == 1)
+		signz *= -1;
+	    if (signx == 1)
+		signy *= -1;
+	    signx *= -1;
+	}
+	sets[i] = bestset;
+	size[bestset] += graph[i]->vwgt;
+    }
+}

contrib/Chaco/bpmatch/map2d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+#include	"structs.h"
+
+
+static void free2d();
+
+void      map2d(graph, xvecs, nvtxs, sets, goal, vwgt_max)
+struct vtx_data **graph;	/* data structure with vertex weights */
+double  **xvecs;		/* vectors to partition */
+int       nvtxs;		/* number of vertices */
+short    *sets;			/* set each vertex gets assigned to */
+double   *goal;			/* desired set sizes */
+int       vwgt_max;		/* largest vertex weight */
+{
+    extern int DEBUG_BPMATCH;	/* turn on debuging for bipartite matching */
+    extern int N_VTX_MOVES;	/* total number of vertex moves */
+    extern int N_VTX_CHECKS;	/* total number of moves contemplated */
+    double   *vals[4][MAXSETS];	/* values in sorted lists */
+    double    dist[4];		/* trial separation point */
+    double    size[4];		/* sizes of each set being modified */
+    int      *indices[4][MAXSETS];	/* indices sorting lists */
+    int       startvtx[4][MAXSETS];	/* indices defining separation */
+    int       nsection = 2;	/* number of xvectors */
+    int       nsets = 4;	/* number of sets being divided into */
+    void      genvals2d(), sorts2d(), inits2d(), checkbp(), movevtxs();
+
+    N_VTX_CHECKS = N_VTX_MOVES = 0;
+
+    /* Generate all the lists of values that need to be sorted. */
+    genvals2d(xvecs, vals, nvtxs);
+
+    /* Sort the lists of values. */
+    sorts2d(vals, indices, nvtxs);
+
+    /* Now initialize dists and assign to sets. */
+    inits2d(graph, xvecs, vals, indices, nvtxs, dist, startvtx, size, sets);
+
+    /* Determine the largest and smallest allowed set sizes. */
+    /* (For now, assume all sets must be same size, but can easily change.) */
+
+    if (DEBUG_BPMATCH > 1) {
+	printf(" Calling check before movevtxs\n");
+	checkbp(graph, xvecs, sets, dist, nvtxs, nsection);
+    }
+
+    movevtxs(graph, nvtxs, nsets, dist, indices, vals, startvtx, sets, size,
+	     goal, vwgt_max);
+
+    if (DEBUG_BPMATCH > 0) {
+	printf(" N_VTX_CHECKS = %d, N_VTX_MOVES = %d\n", N_VTX_CHECKS, N_VTX_MOVES);
+	checkbp(graph, xvecs, sets, dist, nvtxs, nsection);
+    }
+
+    free2d(vals, indices);
+}
+
+
+/* Free the space used in the bpmatch routines. */
+static void free2d(vals, indices)
+double   *vals[4][MAXSETS];
+int      *indices[4][MAXSETS];
+{
+    int       sfree();
+
+    sfree((char *) vals[0][1]);
+    sfree((char *) vals[0][2]);
+    sfree((char *) vals[0][3]);
+    sfree((char *) vals[1][2]);
+
+    sfree((char *) indices[0][1]);
+    sfree((char *) indices[0][2]);
+    sfree((char *) indices[0][3]);
+    sfree((char *) indices[1][2]);
+}

contrib/Chaco/bpmatch/map3d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+#include	"structs.h"
+
+
+static void free3d();
+
+void      map3d(graph, xvecs, nvtxs, sets, goal, vwgt_max)
+struct vtx_data **graph;	/* graph data structure */
+double  **xvecs;		/* vectors to partition */
+int       nvtxs;		/* number of vertices */
+short    *sets;			/* set each vertex gets assigned to */
+double   *goal;			/* desired set sizes */
+int       vwgt_max;		/* largest vertex weight */
+{
+    extern int DEBUG_BPMATCH;	/* debug flag for bipartite matching */
+    extern int N_VTX_MOVES;	/* number of vertices moved between sets */
+    extern int N_VTX_CHECKS;	/* number of vertex moves considered */
+    double   *vals[8][MAXSETS];	/* values in sorted lists */
+    double    dist[8];		/* trial separation point */
+    double    size[8];		/* sizes of each set being modified */
+    int      *indices[8][MAXSETS];	/* indices sorting lists */
+    int       startvtx[8][MAXSETS];	/* indices defining separation */
+    int       nsection = 3;	/* number of xvectors */
+    int       nsets = 8;	/* number of sets being divided into */
+    void      genvals3d(), sorts3d(), inits3d(), checkbp(), movevtxs();
+
+    N_VTX_CHECKS = N_VTX_MOVES = 0;
+
+    /* Generate all the lists of values that need to be sorted. */
+    genvals3d(xvecs, vals, nvtxs);
+
+    /* Sort the lists of values. */
+    sorts3d(vals, indices, nvtxs);
+
+    /* Now initialize distances using median values, and assign to sets. */
+    inits3d(graph, xvecs, vals, indices, nvtxs, dist, startvtx, size, sets);
+
+    /* Determine largest and smallest allowed sizes for each set. */
+    /* (For now all are the same, but easily changed.) */
+
+    if (DEBUG_BPMATCH > 1) {
+	printf(" Calling check before movevtxs\n");
+	checkbp(graph, xvecs, sets, dist, nvtxs, nsection);
+    }
+
+    movevtxs(graph, nvtxs, nsets, dist, indices, vals, startvtx, sets, size,
+	     goal, vwgt_max);
+
+    if (DEBUG_BPMATCH > 0) {
+	printf(" N_VTX_CHECKS = %d, N_VTX_MOVES = %d\n", N_VTX_CHECKS, N_VTX_MOVES);
+	checkbp(graph, xvecs, sets, dist, nvtxs, nsection);
+    }
+
+    free3d(vals, indices);
+}
+
+
+static void free3d(vals, indices)
+double   *vals[8][MAXSETS];
+int      *indices[8][MAXSETS];
+{
+    int       sfree();
+
+    sfree((char *) vals[0][1]);
+    sfree((char *) vals[0][2]);
+    sfree((char *) vals[0][4]);
+    sfree((char *) vals[0][3]);
+    sfree((char *) vals[1][2]);
+    sfree((char *) vals[0][5]);
+    sfree((char *) vals[1][4]);
+    sfree((char *) vals[0][6]);
+    sfree((char *) vals[2][4]);
+    sfree((char *) vals[0][7]);
+    sfree((char *) vals[1][6]);
+    sfree((char *) vals[2][5]);
+    sfree((char *) vals[3][4]);
+
+    sfree((char *) indices[0][1]);
+    sfree((char *) indices[0][2]);
+    sfree((char *) indices[0][4]);
+    sfree((char *) indices[0][3]);
+    sfree((char *) indices[1][2]);
+    sfree((char *) indices[0][5]);
+    sfree((char *) indices[1][4]);
+    sfree((char *) indices[0][6]);
+    sfree((char *) indices[2][4]);
+    sfree((char *) indices[0][7]);
+    sfree((char *) indices[1][6]);
+    sfree((char *) indices[2][5]);
+    sfree((char *) indices[3][4]);
+}

contrib/Chaco/bpmatch/movevtxs.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"params.h"
+#include	"structs.h"
+#include	"defs.h"
+
+
+/* Enact a single move from largest (smallest) set in such a way
+   that some small (large) set gets bigger. */
+
+int       N_VTX_CHECKS;		/* number of considered moves */
+int       N_VTX_MOVES;		/* number of actual moves */
+static void nextmove(), undo_coupling(), couple();
+
+void      movevtxs(graph, nvtxs, nsets, dist, indices, vals, startvtx, sets,
+		             size, goal, vwgt_max)
+struct vtx_data **graph;	/* data structure with vertex weights */
+int       nvtxs;		/* number of vertices in graph */
+int       nsets;		/* how many sets am I dividing into? */
+double   *dist;			/* distances defining splitter */
+int      *indices[][MAXSETS];	/* indices that define order in sorted lists */
+double   *vals[][MAXSETS];	/* values in sorted lists */
+int       startvtx[][MAXSETS];	/* index values corresponding to splitter */
+short    *sets;			/* set assignment for each vertex */
+double   *size;			/* sizes of the different sets */
+double   *goal;			/* desired set sizes */
+int       vwgt_max;		/* largest vertex weight */
+{
+    double    largest;		/* largest overshoot from desired size */
+    double    smallest;		/* largest undershoot from desired size */
+    int       active[MAXSETS];	/* flags sets trying to change size */
+    double    delta;		/* amount distances must change */
+    int       vtx;		/* vertex being moved */
+    short     to, from;		/* set vertex is being moved to/from */
+    int       weight;		/* weight of vertex being moved */
+    int       done;		/* have I successfully move a vertex? */
+
+    /* int npass=0; *//* counts passes through main loop */
+    int       badset;		/* most unbalanced set */
+    int       toobig;		/* badset too large or too small? */
+    int       balanced;		/* is balance attained? */
+    double    imbalance;	/* amount of imbalance in badset */
+    int       i;		/* loop counter */
+
+    /* Find most unbalanced set. */
+
+    imbalance = largest = smallest = 0;
+    for (i = 0; i < nsets; i++) {
+	if (size[i] - goal[i] > largest) {
+	    largest = size[i] - goal[i];
+	    if (largest > imbalance) {
+		imbalance = largest;
+		badset = i;
+		toobig = 1;
+	    }
+	}
+	else if (goal[i] - size[i] > smallest) {
+	    smallest = goal[i] - size[i];
+	    if (smallest > imbalance) {
+		imbalance = smallest;
+		badset = i;
+		toobig = -1;
+	    }
+	}
+    }
+    if (largest + smallest <= vwgt_max)
+	balanced = TRUE;
+    else
+	balanced = FALSE;
+
+
+    /* If not balanced, change distances to move vertices between sets. */
+    while (!balanced) {
+	/* npass++; */
+	for (i = 0; i < nsets; i++)
+	    active[i] = FALSE;
+	active[badset] = TRUE;
+
+	done = FALSE;
+	while (!done) {
+	    nextmove(nvtxs, nsets, vals, indices, startvtx, dist, sets,
+		     toobig, active, &vtx, &to, &delta);
+	    from = sets[vtx];
+	    weight = graph[vtx]->vwgt;
+
+	    /* Now adjust all active dists to reflect this move so far. */
+	    for (i = 0; i < nsets; i++)
+		if (active[i])
+		    dist[i] -= toobig * delta;
+
+	    if (toobig > 0) {
+		if (size[to] + weight - goal[to] < largest) {
+		    done = TRUE;
+		    size[from] -= graph[vtx]->vwgt;
+		    size[to] += graph[vtx]->vwgt;
+		    sets[vtx] = to;
+		    undo_coupling(graph, sets, nsets, from, to, toobig, badset,
+				  size);
+		}
+		else {
+		    couple(nsets, from, to, vtx);
+		    active[to] = TRUE;
+		}
+	    }
+	    else {
+		if (goal[from] - (size[from] - weight) < smallest) {
+		    done = TRUE;
+		    size[from] -= graph[vtx]->vwgt;
+		    size[to] += graph[vtx]->vwgt;
+		    sets[vtx] = to;
+		    undo_coupling(graph, sets, nsets, from, to, toobig, badset,
+				  size);
+		}
+		else {
+		    couple(nsets, from, to, vtx);
+		    active[from] = TRUE;
+		}
+	    }
+	}
+
+	/* Find most unbalanced set. */
+	imbalance = largest = smallest = 0;
+	for (i = 0; i < nsets; i++) {
+	    if (size[i] - goal[i] > largest) {
+		largest = size[i] - goal[i];
+		if (largest > imbalance) {
+		    imbalance = largest;
+		    badset = i;
+		    toobig = 1;
+		}
+	    }
+	    else if (goal[i] - size[i] > smallest) {
+		smallest = goal[i] - size[i];
+		if (smallest > imbalance) {
+		    imbalance = smallest;
+		    badset = i;
+		    toobig = -1;
+		}
+	    }
+	}
+	if (largest + smallest <= vwgt_max)
+	    balanced = TRUE;
+	else
+	    balanced = FALSE;
+    }
+}
+
+
+
+
+/* Find the next move from an active to an inactive set, returning */
+/* next_vtx, next_to, and distance. */
+static void nextmove(nvtxs, nsets, vals, indices, startvtx, dist, sets,
+		               toobig, active, next_vtx, next_to, next_delta)
+int       nvtxs;		/* number of vertices in graph */
+int       nsets;		/* how many sets am I dividing into? */
+double   *dist;			/* distances defining splitter */
+int      *indices[][MAXSETS];	/* indices that define order in sorted lists */
+double   *vals[][MAXSETS];	/* values in sorted lists */
+int       startvtx[][MAXSETS];	/* index values corresponding to splitter */
+short    *sets;			/* set assignment for each vertex */
+int       toobig;		/* is bad set too big or too small? */
+int      *active;		/* flags sets trying to change size */
+int      *next_vtx;		/* vertex selected to move next */
+short    *next_to;		/* set vertex should be moved to */
+double   *next_delta;		/* size of change in distances */
+{
+    double    delta;		/* amount distance must change */
+    double    bestdelta;	/* best value see so far */
+    int       good;		/* is this the delta OK? */
+    int       first;		/* is this the first OK delta I've seen? */
+    int       maxset, minset;	/* larger/smaller of two sets */
+    short     bestfrom, bestto;	/* sets best move comes from and goes to */
+    int       bestdir;		/* direction to step in list for best move */
+    int       bestvtx;		/* vertex being moved between sets */
+    short     from, to;		/* sets vertex wants to move from and to */
+    int       index;		/* offset into indices array */
+    int       dir;		/* direction to step through list */
+    int       i, j;		/* loop counter */
+
+    bestdelta = 0;
+    first = TRUE;
+    while (first) {
+	for (i = 0; i < nsets; i++)
+	    if (active[i]) {
+		for (j = 0; j < nsets; j++)
+		    if (!active[j]) {
+			/* Look for next move from active set i to inactive set j. */
+			if (toobig > 0) {
+			    from = (short) i;
+			    to = (short) j;
+			}
+			else {
+			    from = (short) j;
+			    to = (short) i;
+			}
+			minset = min(to, from);
+			maxset = max(to, from);
+
+			index = startvtx[minset][maxset];
+			if (index >= nvtxs || index < 0)
+			    good = FALSE;
+
+			else {
+			    if (j > i) {
+				dir = -1;
+			    }
+			    else {
+				dir = 1;
+			    }
+			    good = TRUE;
+			    delta = -toobig * dir *
+			       vals[minset][maxset][indices[minset][maxset][index]] -
+			       (dist[to] - dist[from]);
+			}
+
+			if (good && (first || delta < bestdelta)) {
+			    /* Is this the best so far? */
+			    first = FALSE;
+			    bestdelta = delta;
+			    bestfrom = from;
+			    bestto = to;
+			    bestdir = dir;
+			    bestvtx = indices[minset][maxset][index] + 1;
+			}
+		    }
+	    }
+
+	/* Only accept a vertex if it's from the right set. */
+	if (sets[bestvtx] != bestfrom || (toobig > 0 && !active[bestfrom]) ||
+		(toobig < 0 && !active[bestto])) {
+	    /* Set rejection flag, and increment startvtx pointer. */
+	    first = TRUE;
+	    minset = min(bestto, bestfrom);
+	    maxset = max(bestto, bestfrom);
+	    startvtx[minset][maxset] -= toobig * bestdir;
+	}
+	++N_VTX_CHECKS;
+    }
+    *next_vtx = bestvtx;
+    *next_to = bestto;
+    *next_delta = bestdelta;
+    ++N_VTX_MOVES;
+}
+
+
+
+static int ncoupled = 0;
+static int coupled_vtxs[MAXSETS];
+static short coupled_sets[MAXSETS][MAXSETS];
+
+static void couple(nsets, from, to, vtx)
+int       nsets;		/* number of sets being divided into */
+short     from, to;		/* sets to be coupled */
+int       vtx;			/* vertex that they share */
+{
+    int       i;		/* loop counter */
+
+    /* Check for degenerate case of vertex shared by more than two sets. */
+    for (i = 0; i < ncoupled; i++) {
+	if (coupled_vtxs[i] == vtx) {
+	    coupled_sets[i][from] = -1;
+	    coupled_sets[i][to] = 1;
+	    return;
+	}
+    }
+
+    coupled_vtxs[ncoupled] = vtx;
+    for (i = 0; i < nsets; i++)
+	coupled_sets[ncoupled][i] = 0;
+    coupled_sets[ncoupled][from] = -1;
+    coupled_sets[ncoupled][to] = 1;
+    ++ncoupled;
+}
+
+
+static void undo_coupling(graph, sets, nsets, from, to, toobig, badset, size)
+struct vtx_data **graph;	/* data structure with vertex weights */
+short    *sets;			/* sets each vertex is in */
+int       nsets;		/* number of sets being divided into */
+short     from, to;		/* set final vertex moved from and to */
+int       toobig;		/* are we shrinking or enlarging a set? */
+short     badset;		/* the set number being shrunk or enlarged */
+double   *size;			/* sizes of the different sets */
+{
+    int       done;		/* have enough vertices been moved? */
+    int       found;		/* have I found the right set? */
+    int       vtx;		/* vertex being moved between sets */
+    int       i, j;		/* loop counter */
+
+
+    if (ncoupled == 0)
+	return;
+
+    if (toobig > 0) {
+	done = FALSE;
+	if (from == badset)
+	    done = TRUE;
+	while (!done) {
+	    found = FALSE;
+	    to = from;
+	    for (i = 0; i < ncoupled && !found; i++) {
+		if (coupled_sets[i][from] == 1) {	/* Found an edge into set. */
+		    found = TRUE;
+		    /* And find other end of edge. */
+		    for (j = 0; j < nsets; j++) {
+			if (coupled_sets[i][j] == -1)
+			    from = (short) j;
+		    }
+
+		    /* Switch the set for this vertex. */
+		    vtx = coupled_vtxs[i];
+		    size[from] -= graph[vtx]->vwgt;
+		    size[to] += graph[vtx]->vwgt;
+		    sets[vtx] = to;
+		    coupled_sets[i][from] = 0;
+
+		    if (from == badset)
+			done = TRUE;
+		}
+	    }
+	}
+    }
+    else {
+	done = FALSE;
+	if (to == badset)
+	    done = TRUE;
+	while (!done) {
+	    found = FALSE;
+	    from = to;
+	    for (i = 0; i < ncoupled && !found; i++) {
+		if (coupled_sets[i][from] == -1) {	/* Found an edge from set. */
+		    found = TRUE;
+		    /* And find other end of edge. */
+		    for (j = 0; j < nsets; j++) {
+			if (coupled_sets[i][j] == 1)
+			    to = (short) j;
+		    }
+
+		    /* Switch the set for this vertex. */
+		    vtx = coupled_vtxs[i];
+		    size[from] -= graph[vtx]->vwgt;
+		    size[to] += graph[vtx]->vwgt;
+		    sets[vtx] = to;
+		    coupled_sets[i][to] = 0;
+
+		    if (to == badset)
+			done = TRUE;
+		}
+	    }
+	}
+    }
+
+    ncoupled = 0;
+}

contrib/Chaco/bpmatch/sorts2d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+
+
+void      sorts2d(vals, indices, nvtxs)
+/* Sort the lists needed to find the splitter. */
+double   *vals[4][MAXSETS];	/* lists of values to sort */
+int      *indices[4][MAXSETS];	/* indices of sorted lists */
+int       nvtxs;		/* number of vertices */
+{
+    int      *space;		/* space for mergesort routine */
+    int      *temp[4];		/* place holders for indices */
+    int       nlists = 4;	/* number of directions to sort */
+    int       i;		/* loop counter */
+    double   *smalloc();
+    int       sfree();
+    void      mergesort();
+
+    space = (int *) smalloc((unsigned) nvtxs * sizeof(int));
+
+    for (i = 0; i < nlists; i++) {
+	temp[i] = (int *) smalloc((unsigned) nvtxs * sizeof(int));
+    }
+
+    mergesort(vals[0][1], nvtxs, temp[0], space);
+    mergesort(vals[0][2], nvtxs, temp[1], space);
+    mergesort(vals[0][3], nvtxs, temp[2], space);
+    mergesort(vals[1][2], nvtxs, temp[3], space);
+
+    sfree((char *) space);
+
+    indices[0][1] = indices[1][0] = indices[2][3] = indices[3][2] = temp[0];
+    indices[0][2] = indices[2][0] = indices[1][3] = indices[3][1] = temp[1];
+    indices[0][3] = indices[3][0] = temp[2];
+    indices[1][2] = indices[2][1] = temp[3];
+}

contrib/Chaco/bpmatch/sorts3d.c

+/* This software was developed by Bruce Hendrickson and Robert Leland   *
+ * at Sandia National Laboratories under US Department of Energy        *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include	<stdio.h>
+#include	"defs.h"
+#include	"params.h"
+
+
+void      sorts3d(vals, indices, nvtxs)
+/* Sort the lists needed to find the splitter. */
+double   *vals[8][MAXSETS];	/* lists of values to sort */
+int      *indices[8][MAXSETS];	/* indices of sorted lists */
+int       nvtxs;		/* number of vertices */
+{
+    int      *space;		/* space for mergesort routine */
+    int       nsets = 8;	/* number of sets */
+    int       nlists = 13;	/* number of directions to sort */
+    int      *temp[13];		/* place holders for indices */
+    int       i, j;		/* loop counter */
+    double   *smalloc();
+    int       sfree();
+    void      mergesort();
+
+    space = (int *) smalloc((unsigned) nvtxs * sizeof(int));
+
+    for (i = 0; i < nlists; i++) {
+	temp[i] = (int *) smalloc((unsigned) nvtxs * sizeof(int));
+    }
+
+    mergesort(vals[0][1], nvtxs, temp[0], space);
+    mergesort(vals[0][2], nvtxs, temp[1], space);
+    mergesort(vals[0][4], nvtxs, temp[2], space);
+    mergesort(vals[0][3], nvtxs, temp[3], space);
+    mergesort(vals[1][2], nvtxs, temp[4], space);
+    mergesort(vals[0][5], nvtxs, temp[5], space);
+    mergesort(vals[1][4], nvtxs, temp[6], space);
+    mergesort(vals[0][6], nvtxs, temp[7], space);
+    mergesort(vals[2][4], nvtxs, temp[8], space);
+    mergesort(vals[0][7], nvtxs, temp[9], space);
+    mergesort(vals[1][6], nvtxs, temp[10], space);
+    mergesort(vals[2][5], nvtxs, temp[11], space);
+    mergesort(vals[3][4], nvtxs, temp[12], space);
+
+    sfree((char *) space);
+
+    indices[0][1] = indices[2][3] = indices[4][5] = indices[6][7] = temp[0];
+    indices[0][2] = indices[1][3] = indices[4][6] = indices[5][7] = temp[1];
+    indices[0][4] = indices[1][5] = indices[2][6] = indices[3][7] = temp[2];
+    indices[0][3] = indices[4][7] = temp[3];
+    indices[1][2] = indices[5][6] = temp[4];
+    indices[0][5] = indices[2][7] = temp[5];
+    indices[1][4] = indices[3][6] = temp[6];
+    indices[0][6] = indices[1][7] = temp[7];
+    indices[2][4] = indices[3][5] = temp[8];
+    indices[0][7] = temp[9];
+    indices[1][6] = temp[10];
+    indices[2][5] = temp[11];
+    indices[3][4] = temp[12];
+
+    for (i = 0; i < nsets; i++) {
+	for (j = i + 1; j < nsets; j++)
+	    indices[j][i] = indices[i][j];
+    }
+}