Index: /issm/trunk/externalpackages/meshpart/Chaco.c
===================================================================
--- /issm/trunk/externalpackages/meshpart/Chaco.c	(revision 3997)
+++ /issm/trunk/externalpackages/meshpart/Chaco.c	(revision 3997)
@@ -0,0 +1,352 @@
+/*
+  CHACO : Hendrickson/Leland's graph partitioner.
+
+  [part1,part2,chaco_time] = chaco(A) returns a 50/50 vertex partition of the mesh
+  whose symmetric adjacency matrix is A.  
+
+  Optional arguments:
+  map = chaco(A, xy, method, nparts, goal);
+  map = chaco(A, xy, method, inmap, goal);
+
+  A:          Depending on "method", A may contain vertex and edge weights.
+
+  xy:         Each row of xy is the coordinates of a vertex.
+              If xy is non-null and there is no output, draw a picture.
+
+  method:     Scalar or vector describing the desired method.  
+              Default is multilevel Kernighan-Lin; other possibilities below.
+
+  nparts      Number of parts to divide into.  Default is 2.  If nparts is 
+    or        present, the output is a "map vector", see below.  (If method(5) 
+  inmap:      is specified, nparts is interpreted differently; see below.  In
+              any case, the default is to divide into two parts.)
+              If method(1) = 7 (see below), this argument is a map vector
+              specifying an initial 2-way partition, and Chaco refines it.
+
+  goal:       Optionally, a vector of desired sizes (or total vertex weights)
+              for each of the nparts parts.  Default is all sizes equal.
+
+  map:        If nparts and inmap are not present, the output is a vector of 
+              the n/2 vertex numbers in one part of the 2-way partition, for
+              compatibility with geopart and specpart.
+              If nparts or imap is present, the output is a vector of the
+              n part numbers, from 0 to nparts-1, assigned to the vertices.
+
+  This is a Matlab interface to the graph partitioning software described
+  in B. Hendrickson and R. Leland, "The Chaco User's Guide (Version 2.0)",
+  Sandia National Laboratories report SAND94-2692, October 1994.
+  This interface was written by John Gilbert, Xerox PARC, and is
+  Copyright (c) 1994-1996 by Xerox Corporation.  All rights reserved.
+  HELP COPYRIGHT for complete copyright and licensing notice.
+
+  Modified by Tim Davis, for Matlab 5.1.  July 6, 1998.
+  Modified by Tim Davis, July 1999.  Added 2nd output argument,
+  which is the cpu time taken by chaco itself.
+    05/25/10    jes    Reorganization into Matlab-layer and x-layer.
+
+  See also GEOPART, SPECPART.
+
+  "method" is a vector of flags as follows.  Not all combinations are supported.
+  See Section 6.10 of the Chaco manual for more details on all the arguments.
+  If "method" is shorter than 10, we use the defaults for unspecified entries.
+
+  method[0]:  Global partitioning method  ("global_method" in the Chaco manual).
+              1 Multilevel Kernighan-Lin (default)
+              2 Spectral
+              3 Inertial
+              4 Linear
+              5 Random
+              6 Scattered
+              7 Use "inmap" as the global (2-way) partition
+
+  method[1]:  Local refinement method  ("local_method" in the Chaco manual).
+              1 Kernighan-Lin (default)
+              2 None
+
+  method[2]:  Vertex weighting.
+              0 No weights (default)
+              1 Use diag(A) as (positive integer) vertex weights
+
+  method[3]:  Edge weighting.
+              0 No weights (default)
+              1 Use off-diagonals of A as (positive integer) edge weights
+
+  method[4]:  Target architecture  ("architecture" in the Chaco manual).
+              If method[4] = 0, the target is a hypercube, "nparts" is the 
+              number of dimensions, and the partition is into 2^nparts parts.  
+              If method[4] = 1, 2, or 3, the target is a 1-, 2-, or 3-D grid,
+              "nparts" is a vector of the sizes of the grid in each dimension,
+              and the partition is into prod(nparts) parts.
+              Default is method[4] = 1, so nparts is the number of parts.
+
+  method[5]:  Partitioning dimension  ("ndims" in the Chaco manual).
+              1 Bisection (default)
+              2 Quadrisection
+              3 Octasection
+
+  method[6]:  Number of vertices to coarsen to  ("vmax" in the Chaco manual).
+              Default is 50.
+
+  method[7]:  Eigensolver  ("rqi_flag" in the Chaco manual).
+              0 RQI/Symmlq (default)
+              1 Lanczos 
+
+  method[8]:  Eigensolver convergence tolerance  ("eigtol" in the Chaco manual).
+              Default is .001
+
+  method[9]:  Seed for random number generator  ("seed" in the Chaco manual).
+              Default is 7654321.
+
+  Many esoteric details of Chaco's behavior can be changed by placing a file
+  called "User_Params" in the same directory as the executable mlchaco.mex.
+  As always, see the manual for details.
+*/
+
+#include <stdio.h>
+#include "mex.h"
+#include <sys/times.h>
+#include <time.h>      /*  CLOCKS_PER_SEC  */
+
+/* Aliases for input and output arguments */
+#define A_IN         prhs[0]
+#define VWGTS_IN     prhs[1]
+#define EWGTS_IN     prhs[2]
+#define XYZ_IN       prhs[3]
+#define METHOD_IN    prhs[4]
+#define NPARTS_IN    prhs[5]
+#define GOAL_IN      prhs[6]
+#define MAP_OUT      plhs[0]
+#define TIME_OUT     plhs[1]
+
+#define THISFUNCTION "Chaco"
+
+void ChacoUsage( void );
+
+int Chacox(
+	int       nvtxs,          /* number of vertices in full graph */
+	int      *start,          /* start of edge list for each vertex */
+	int      *adjacency,      /* edge list data */
+	int      *vwgts,          /* weights for all vertices */
+	float    *ewgts,          /* weights for all edges */
+	float    *x,
+	float    *y,
+	float    *z,              /* coordinates for inertial method */
+	short    *assignment,     /* set number of each vtx (length n) */
+	double    method[10],     /* vector describing the desired method */
+	int      *nparts,         /* number of parts to divide into */
+	double   *goal            /* desired set sizes for each set */
+);
+
+
+void mexFunction(
+	int           nlhs,           /* number of expected outputs */
+	mxArray       *plhs[],        /* array of pointers to output arguments */
+	int           nrhs,           /* number of inputs */
+	const mxArray *prhs[]         /* array of pointers to input arguments */
+)
+{
+	/* Arguments to Chaco "interface" routine: */
+
+	int       nvtxs;          /* number of vertices in full graph */
+	int      *start;          /* start of edge list for each vertex */
+	int      *adjacency;      /* edge list data */
+	int      *vwgts;          /* weights for all vertices */
+	float    *ewgts;          /* weights for all edges */
+	float    *x;
+	float    *y;
+	float    *z;              /* coordinates for inertial method */
+	short    *assignment;     /* set number of each vtx (length n) */
+	double    method[10]={1,1,0,0,1,1,50,0,.001,7654321};
+							  /* vector describing the desired method */
+	int      *nparts;         /* number of parts to divide into */
+	double   *goal;           /* desired set sizes for each set */
+
+	int    i, nedges, nterms, geodims, failure;
+	double *p;
+	mwIndex *mwstart, *mwadjacency;
+
+	double *Tout ;			/* time: added by Tim Davis */
+	struct tms t1 ;
+	struct tms t2 ;
+
+	/* Check for correct number of arguments passed from Matlab. */
+
+	if      (nrhs == 0 && nlhs == 0) {
+		ChacoUsage();
+		return;
+	} else if (nrhs < 1) {
+		ChacoUsage();
+		mexErrMsgTxt("Chaco mex-layer requires at least one input argument.");
+	} else if (nlhs > 2) {
+		ChacoUsage();
+		mexErrMsgTxt("Chaco mex-layer requires one or two output arguments.");
+	}
+
+	/* Slog through the arguments, converting from Matlab matrices
+	   to the various types Chaco wants.                            */
+
+	if (!mxIsEmpty(A_IN) && mxIsNumeric(A_IN) && mxIsSparse(A_IN)) {
+		nvtxs = mxGetN(A_IN);
+		mwstart = mxGetJc(A_IN);
+		start = mxMalloc((mxGetN(A_IN)+1)*sizeof(int));
+		for (i=0; i<(mxGetN(A_IN)+1); i++)
+			start[i]= (int)mwstart[i];
+		nedges = start[nvtxs];
+		mwadjacency = mxGetIr(A_IN);
+		adjacency = mxMalloc(mxGetNzmax(A_IN)*sizeof(int));
+		for (i=0; i<mxGetNzmax(A_IN); i++)
+			adjacency[i]= (int)mwadjacency[i];
+	}
+	else {
+		mexPrintf("%s -- Adjacency matrix must be numeric and sparse.\n",
+				  THISFUNCTION);
+		mexErrMsgTxt(" ");
+	}
+
+	if (nrhs >= 2 && !mxIsEmpty(VWGTS_IN)) {
+		if (mxIsNumeric(VWGTS_IN) && !mxIsSparse(VWGTS_IN) &&
+			(nterms=mxGetM(VWGTS_IN)*mxGetN(VWGTS_IN)) == nvtxs) {
+			vwgts = (int *) mxCalloc(nvtxs, sizeof(int));
+			p = mxGetPr(VWGTS_IN);
+			for (i = 0; i < nvtxs; vwgts[i++] = (int) *p++);
+		}
+		else {
+			mexPrintf("%s -- Vertex weight vector must be numeric, full, and length=%d.\n",
+					  THISFUNCTION,nvtxs);
+			mexErrMsgTxt(" ");
+		}
+	} else vwgts = NULL;
+
+	if (nrhs >= 3 && !mxIsEmpty(EWGTS_IN)) {
+		if (mxIsNumeric(EWGTS_IN) && mxIsSparse(EWGTS_IN) &&
+			(nterms=mxGetNzmax(EWGTS_IN)) == nedges) {
+			ewgts = (float *) mxCalloc(nedges, sizeof(float));
+			p = mxGetPr(A_IN);
+			for (i = 0; i < nedges; ewgts[i++] = (float) *p++);
+		}
+		else {
+			mexPrintf("%s -- Edge weight matrix must be numeric, sparse, and nonzeroes=%d.\n",
+					  THISFUNCTION,nedges);
+			mexErrMsgTxt(" ");
+		}
+	} else ewgts = NULL;
+
+	if (nrhs >= 4 && (geodims = mxGetN(XYZ_IN)) >= 1) {
+		if (mxIsNumeric(XYZ_IN) && !mxIsSparse(XYZ_IN) &&
+			mxGetM(XYZ_IN) == nvtxs) {
+			x = (float *) mxCalloc(nvtxs, sizeof(float));
+			p = mxGetPr(XYZ_IN);
+			for (i = 0; i < nvtxs; x[i++] = (float) *p++);
+		}
+		else {
+			mexPrintf("%s -- XYZ coordinate vectors must be numeric, full, and length=%d.\n",
+					  THISFUNCTION,nvtxs);
+			mexErrMsgTxt(" ");
+		}
+	} else x = NULL;
+	if (nrhs >= 4 && (geodims                 ) >= 2) {
+		y = (float *) mxCalloc(nvtxs, sizeof(float));
+		for (i = 0; i < nvtxs; y[i++] = (float) *p++);
+	} else y = NULL;
+	if (nrhs >= 4 && (geodims                 ) >= 3) {
+		z = (float *) mxCalloc(nvtxs, sizeof(float));
+		for (i = 0; i < nvtxs; z[i++] = (float) *p++);
+	} else z = NULL;
+
+	if (nrhs >= 5 && !mxIsEmpty(METHOD_IN)) {
+		if (mxIsNumeric(METHOD_IN) && !mxIsSparse(METHOD_IN) &&
+			(nterms=mxGetM(METHOD_IN)*mxGetN(METHOD_IN))) {
+			p = mxGetPr(METHOD_IN);
+			for (i = 0; i < (nterms<10 ? nterms : 10); method[i++] = *p++);
+		}
+		else {
+			mexPrintf("%s -- Method vector must be numeric and full.\n",
+					  THISFUNCTION);
+			mexErrMsgTxt(" ");
+		}
+	}
+
+	if (nrhs >= 6 && !mxIsEmpty(NPARTS_IN)) {
+		if (mxIsNumeric(NPARTS_IN) && !mxIsSparse(NPARTS_IN) &&
+			(nterms=mxGetM(NPARTS_IN)*mxGetN(NPARTS_IN))) {
+			nparts = (int *) mxCalloc(nterms, sizeof(int));
+			p = mxGetPr(NPARTS_IN);
+			for (i = 0; i < nterms; nparts[i++] = (int) *p++);
+		}
+		else {
+			mexPrintf("%s -- Parts vector must be numeric and full.\n",
+					  THISFUNCTION);
+			mexErrMsgTxt(" ");
+		}
+	} else nparts = NULL;
+
+	if (nrhs >= 7 && !mxIsEmpty(GOAL_IN)) {
+		if (mxIsNumeric(GOAL_IN) && !mxIsSparse(GOAL_IN) &&
+			(nterms=mxGetM(GOAL_IN)*mxGetN(GOAL_IN))) {
+			goal = (double *) mxCalloc(nterms, sizeof(double));
+			p = mxGetPr(GOAL_IN);
+			for (i = 0; i < nterms; goal[i++] = *p++);
+		}
+		else {
+			mexPrintf("%s -- Goal vector must be numeric and full.\n",
+					  THISFUNCTION);
+			mexErrMsgTxt(" ");
+		}
+	} else goal = NULL;
+
+	assignment = (short *) mxCalloc(nvtxs, sizeof(short));
+
+	/* Finally, call Chaco */
+
+	(void) times (&t1) ;
+
+	failure = Chacox(nvtxs, start, adjacency, vwgts, ewgts, x, y, z,
+		assignment, method, nparts, goal);
+
+	(void) times (&t2) ;
+
+	/* Copy the mapping vector into a Matlab matrix. */
+
+	if (!failure) {
+		MAP_OUT = mxCreateDoubleMatrix(1,nvtxs,mxREAL);
+		p = mxGetPr(MAP_OUT);
+		for (i = 0; i < nvtxs; *p++ = (double) assignment[i++]);
+	}
+
+	if (nlhs > 1)
+	{
+		TIME_OUT = mxCreateDoubleMatrix (1, 1, mxREAL) ;
+		Tout = mxGetPr (TIME_OUT) ;
+		Tout [0] =
+			((double) (t2.tms_utime + t2.tms_stime - t1.tms_utime - t1.tms_stime)) /
+			((double) CLOCKS_PER_SEC) ;
+	}
+
+	/* Free what we allocated */
+   
+	if (start != NULL) mxFree((void *) start);
+	if (adjacency != NULL) mxFree((void *) adjacency);
+	if (vwgts != NULL) mxFree((void *) vwgts);
+	if (ewgts != NULL) mxFree((void *) ewgts);
+	if (x != NULL) mxFree((void *) x);
+	if (y != NULL) mxFree((void *) y);
+	if (z != NULL) mxFree((void *) z);
+	if (nparts != NULL) mxFree((void *) nparts);
+	if (goal != NULL) mxFree((void *) goal);
+	if (assignment != NULL) mxFree((void *) assignment);
+
+	if (failure)
+		mexErrMsgTxt("The call to C Chaco returned failure.");
+
+	return;
+}
+
+void ChacoUsage( void )
+{
+	mexPrintf("\n");
+	mexPrintf("Usage: [map,chaco_time] = Chaco(A,vwgts,ewgts,xyz,method,nparts,goal);\n");
+	mexPrintf("\n");
+
+	return;
+}
+
Index: /issm/trunk/externalpackages/meshpart/Chaco_m.m
===================================================================
--- /issm/trunk/externalpackages/meshpart/Chaco_m.m	(revision 3997)
+++ /issm/trunk/externalpackages/meshpart/Chaco_m.m	(revision 3997)
@@ -0,0 +1,217 @@
+function  [part1,part2,chaco_time] = chaco(A,xy,method,nparts,goal)
+% CHACO : Hendrickson/Leland's graph partitioner.
+% 
+%  [part1,part2,chaco_time] = chaco(A) returns a 50/50 vertex partition of the mesh
+%  whose symmetric adjacency matrix is A.  
+%
+%  Optional arguments:
+%  map = chaco(A, xy, method, nparts, goal);
+%  map = chaco(A, xy, method, inmap, goal);
+%
+%  A:          Depending on "method", A may contain vertex and edge weights.
+%
+%  xy:         Each row of xy is the coordinates of a vertex.
+%              If xy is non-null and there is no output, draw a picture.
+%
+%  method:     Scalar or vector describing the desired method.  
+%              Default is multilevel Kernighan-Lin; other possibilities below.
+%
+%  nparts      Number of parts to divide into.  Default is 2.  If nparts is 
+%    or        present, the output is a "map vector", see below.  (If method(5) 
+%  inmap:      is specified, nparts is interpreted differently; see below.  In
+%              any case, the default is to divide into two parts.)
+%              If method(1) = 7 (see below), this argument is a map vector
+%              specifying an initial 2-way partition, and Chaco refines it.
+%
+%  goal:       Optionally, a vector of desired sizes (or total vertex weights)
+%              for each of the nparts parts.  Default is all sizes equal.
+%
+%  map:        If nparts and inmap are not present, the output is a vector of 
+%              the n/2 vertex numbers in one part of the 2-way partition, for
+%              compatibility with geopart and specpart.
+%              If nparts or imap is present, the output is a vector of the
+%              n part numbers, from 0 to nparts-1, assigned to the vertices.
+%
+% This is a Matlab interface to the graph partitioning software described
+% in B. Hendrickson and R. Leland, "The Chaco User's Guide (Version 2.0)",
+% Sandia National Laboratories report SAND94-2692, October 1994.
+% This interface was written by John Gilbert, Xerox PARC, and is
+% Copyright (c) 1994-1996 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+%
+% Modified by Tim Davis, for Matlab 5.1.  July 6, 1998.
+%   05/28/10    jes    Reorganization into Matlab-layer and x-layer.
+%
+% See also GEOPART, SPECPART.
+%
+% "method" is a vector of flags as follows.  Not all combinations are supported.
+% See Section 6.10 of the Chaco manual for more details on all the arguments.
+% If "method" is shorter than 10, we use the defaults for unspecified entries.
+%
+% method(1):  Global partitioning method  ("global_method" in the Chaco manual).
+%             1 Multilevel Kernighan-Lin (default)
+%             2 Spectral
+%             3 Inertial
+%             4 Linear
+%             5 Random
+%             6 Scattered
+%             7 Use "inmap" as the global (2-way) partition
+%
+% method(2):  Local refinement method  ("local_method" in the Chaco manual).
+%             1 Kernighan-Lin (default)
+%             2 None
+%
+% method(3):  Vertex weighting.
+%             0 No weights (default)
+%             1 Use diag(A) as (positive integer) vertex weights
+%
+% method(4):  Edge weighting.
+%             0 No weights (default)
+%             1 Use off-diagonals of A as (positive integer) edge weights
+%
+% method(5):  Target architecture  ("architecture" in the Chaco manual).
+%             If method(5) = 0, the target is a hypercube, "nparts" is the 
+%             number of dimensions, and the partition is into 2^nparts parts.  
+%             If method(5) = 1, 2, or 3, the target is a 1-, 2-, or 3-D grid,
+%             "nparts" is a vector of the sizes of the grid in each dimension,
+%             and the partition is into prod(nparts) parts.
+%             Default is method(5) = 1, so nparts is the number of parts.
+%
+% method(6):  Partitioning dimension  ("ndims" in the Chaco manual).
+%             1 Bisection (default)
+%             2 Quadrisection
+%             3 Octasection
+%
+% method(7):  Number of vertices to coarsen to  ("vmax" in the Chaco manual).
+%             Default is 50.
+%
+% method(8):  Eigensolver  ("rqi_flag" in the Chaco manual).
+%             0 RQI/Symmlq (default)
+%             1 Lanczos 
+%
+% method(9):  Eigensolver convergence tolerance  ("eigtol" in the Chaco manual).
+%             Default is .001
+%
+% method(10): Seed for random number generator  ("seed" in the Chaco manual).
+%             Default is 7654321.
+%
+% Many esoteric details of Chaco's behavior can be changed by placing a file
+% called "User_Params" in the same directory as the executable mlchaco.mex.
+% As always, see the manual for details.
+
+DefaultMethod = [1 1 0 0 1 1 50 0 .001 7654321];
+
+% Fill in default arguments.
+if nargin < 2, xy = []; end;
+if nargin < 3, method = DefaultMethod; end;
+if nargin < 4, nparts = []; end;
+if nargin < 5, goal = []; end;
+if length(method) < length(DefaultMethod)
+    method = [method DefaultMethod(length(method)+1 : length(DefaultMethod))];
+end;
+
+% Decide on output and graphics.
+if (isempty (nparts))
+    mapvector = 0;
+else
+    mapvector = 1;
+end;
+picture = (nargout == 0) & (size(xy,2) >= 2);
+
+% Chaco numbers vertices from 1 and the Matlab sparse data structure 
+% numbers rows from 0, so we add an empty first row to make things line up.
+% This code also makes sure the arg to Chaco will be sparse.
+[n,n] = size(A);
+Adiag = diag(diag(A));
+%Aout = [sparse(1,n) ; A-Adiag];
+Aout = sparse(A-Adiag);
+
+% Make sure all args except the adj matrix are full;
+if issparse(xy)
+    xy = full(xy);
+end;
+if issparse(method)
+    method = full(method);
+end;
+if issparse(nparts)
+    nparts = full(nparts);
+end;
+if issparse(goal)
+    goal = full(goal);
+end;
+
+% Decode "method" to get the actual args to Chaco.
+% Note that "nparts" may correspond to any of several Chaco
+% parameters, depending on the method.
+
+if method(3)
+    vwgts = full(diag(A));
+    totalvwgt = sum(vwgts);
+else
+    vwgts = [];
+    totalvwgt = size(A,2);
+end;
+if method(4)
+    ewgts = Aout;
+else
+    ewgts = [];
+end;
+if method(1) == 7
+    % Refine an input partition: "nparts" is the input partition.
+    % This seems to work only for hypercube architecture,
+    % so we force a 1-D hypercube with 2-way partitioning.
+    nsets = 2;
+elseif method(5) == 0
+    % Partition for hypercube: "nparts" is # of dimensions (default 1).
+    if (isempty (nparts))
+    	nsets = 2^1;
+    else
+    	nsets = 2^nparts;
+    end;
+else
+    % Partition for mesh: "nparts" is vector of mesh sizes in each
+    % dimension, default [2 1 ... 1] with "architecture" dimensions.
+    if (isempty (nparts))
+    	nsets = prod(2);
+    else
+    	nsets = prod(nparts);
+    end;
+end;
+if length(goal) ~= nsets
+    goal = totalvwgt/nsets * ones(1,nsets);
+end;
+
+[map,chaco_time]=Chaco(Aout, vwgts, ewgts, xy, method, nparts, goal);
+
+% Draw the picture.
+if picture
+    if size(xy,2) >= 3 && unique(xy(:,3)) == 0
+        xy=xy(:,1:2);
+    end
+    if length(unique(map)) == 2
+        gplotpart(A,xy,find(map==0));
+    else
+        gplotmap(A,xy,map);
+    end;
+    if     method(1)==1, heading = 'Multilevel Kernighan-Lin';
+    elseif method(1)==2, heading = 'Spectral';
+    elseif method(1)==3, heading = 'Inertial';
+    elseif method(1)==4, heading = 'Linear';
+    elseif method(1)==5, heading = 'Random';
+    elseif method(1)==6, heading = 'Scattered';
+    elseif method(1)==7, heading = 'Input'; 
+    end;
+    heading = [heading ' Partition'];
+    if method(2)==1 & method(1) ~= 1 
+        heading =[heading ' Refined by KL'];
+    end;
+    title(heading);
+end;
+
+% Put output in the right form.
+if mapvector
+    part1 = map;
+else
+    part1 = find(map==0);
+    part2 = find(map==1);
+end;
Index: /issm/trunk/externalpackages/meshpart/Chacox.c
===================================================================
--- /issm/trunk/externalpackages/meshpart/Chacox.c	(revision 3997)
+++ /issm/trunk/externalpackages/meshpart/Chacox.c	(revision 3997)
@@ -0,0 +1,262 @@
+/*
+  CHACO : Hendrickson/Leland's graph partitioner.
+
+  [part1,part2,chaco_time] = chaco(A) returns a 50/50 vertex partition of the mesh
+  whose symmetric adjacency matrix is A.  
+
+  Optional arguments:
+  map = chaco(A, xy, method, nparts, goal);
+  map = chaco(A, xy, method, inmap, goal);
+
+  A:          Depending on "method", A may contain vertex and edge weights.
+
+  xy:         Each row of xy is the coordinates of a vertex.
+              If xy is non-null and there is no output, draw a picture.
+
+  method:     Scalar or vector describing the desired method.  
+              Default is multilevel Kernighan-Lin; other possibilities below.
+
+  nparts      Number of parts to divide into.  Default is 2.  If nparts is 
+    or        present, the output is a "map vector", see below.  (If method(5) 
+  inmap:      is specified, nparts is interpreted differently; see below.  In
+              any case, the default is to divide into two parts.)
+              If method(1) = 7 (see below), this argument is a map vector
+              specifying an initial 2-way partition, and Chaco refines it.
+
+  goal:       Optionally, a vector of desired sizes (or total vertex weights)
+              for each of the nparts parts.  Default is all sizes equal.
+
+  map:        If nparts and inmap are not present, the output is a vector of 
+              the n/2 vertex numbers in one part of the 2-way partition, for
+              compatibility with geopart and specpart.
+              If nparts or imap is present, the output is a vector of the
+              n part numbers, from 0 to nparts-1, assigned to the vertices.
+
+  This is a Matlab interface to the graph partitioning software described
+  in B. Hendrickson and R. Leland, "The Chaco User's Guide (Version 2.0)",
+  Sandia National Laboratories report SAND94-2692, October 1994.
+  This interface was written by John Gilbert, Xerox PARC, and is
+  Copyright (c) 1994-1996 by Xerox Corporation.  All rights reserved.
+  HELP COPYRIGHT for complete copyright and licensing notice.
+
+  Modified by Tim Davis, for Matlab 5.1.  July 6, 1998.
+    05/26/10    jes    Reorganization into Matlab-layer and x-layer.
+
+  See also GEOPART, SPECPART.
+
+  "method" is a vector of flags as follows.  Not all combinations are supported.
+  See Section 6.10 of the Chaco manual for more details on all the arguments.
+  If "method" is shorter than 10, we use the defaults for unspecified entries.
+
+  method[0]:  Global partitioning method  ("global_method" in the Chaco manual).
+              1 Multilevel Kernighan-Lin (default)
+              2 Spectral
+              3 Inertial
+              4 Linear
+              5 Random
+              6 Scattered
+              7 Use "inmap" as the global (2-way) partition
+
+  method[1]:  Local refinement method  ("local_method" in the Chaco manual).
+              1 Kernighan-Lin (default)
+              2 None
+
+  method[2]:  Vertex weighting.
+              0 No weights (default)
+              1 Use diag(A) as (positive integer) vertex weights
+
+  method[3]:  Edge weighting.
+              0 No weights (default)
+              1 Use off-diagonals of A as (positive integer) edge weights
+
+  method[4]:  Target architecture  ("architecture" in the Chaco manual).
+              If method[4] = 0, the target is a hypercube, "nparts" is the 
+              number of dimensions, and the partition is into 2^nparts parts.  
+              If method[4] = 1, 2, or 3, the target is a 1-, 2-, or 3-D grid,
+              "nparts" is a vector of the sizes of the grid in each dimension,
+              and the partition is into prod(nparts) parts.
+              Default is method[4] = 1, so nparts is the number of parts.
+
+  method[5]:  Partitioning dimension  ("ndims" in the Chaco manual).
+              1 Bisection (default)
+              2 Quadrisection
+              3 Octasection
+
+  method[6]:  Number of vertices to coarsen to  ("vmax" in the Chaco manual).
+              Default is 50.
+
+  method[7]:  Eigensolver  ("rqi_flag" in the Chaco manual).
+              0 RQI/Symmlq (default)
+              1 Lanczos 
+
+  method[8]:  Eigensolver convergence tolerance  ("eigtol" in the Chaco manual).
+              Default is .001
+
+  method[9]:  Seed for random number generator  ("seed" in the Chaco manual).
+              Default is 7654321.
+
+  Many esoteric details of Chaco's behavior can be changed by placing a file
+  called "User_Params" in the same directory as the executable mlchaco.mex.
+  As always, see the manual for details.
+*/
+
+#include <stdio.h>
+#include "mex.h"
+
+#define THISFUNCTION "Chacox"
+
+int Chacox(
+	int       nvtxs,          /* number of vertices in full graph */
+	int      *start,          /* start of edge list for each vertex */
+	int      *adjacency,      /* edge list data */
+	int      *vwgts,          /* weights for all vertices */
+	float    *ewgts,          /* weights for all edges */
+	float    *x,
+	float    *y,
+	float    *z,              /* coordinates for inertial method */
+	short    *assignment,     /* set number of each vtx (length n) */
+	double    method[10],     /* vector describing the desired method */
+	int      *nparts,         /* number of parts to divide into */
+	double   *goal            /* desired set sizes for each set */
+)
+{
+	/* Arguments to Chaco "interface" routine: */
+
+	char     *outassignname;  /* name of assignment output file */
+	char     *outfilename;    /* output file name */
+	int       architecture;   /* 0 => hypercube, d => d-dimensional mesh */
+	int       ndims_tot=0;    /* total number of cube dimensions to divide */
+	int       mesh_dims[3]={0,0,0};
+							  /* dimensions of mesh of processors */
+	int       global_method;  /* global partitioning algorithm */
+	int       local_method;   /* local partitioning algorithm */
+	int       rqi_flag;       /* should I use RQI/Symmlq eigensolver? */
+	int       vmax;           /* if so, how many vertices to coarsen down to? */
+	int       ndims;          /* number of eigenvectors (2^d sets) */
+	double    eigtol;         /* tolerance on eigenvectors */
+	long      seed;           /* for random graph mutations */
+ 
+	int    i, nsets, failure;
+	double totalvwgt, totalgoal;
+
+/*	Chaco numbers vertices from 1 and the Matlab sparse data structure 
+	numbers rows from 0, so we add an empty first row to make things line up. */
+
+	for (i=0; i<start[nvtxs]; i++)
+		adjacency[i]++;
+
+	outassignname = NULL;
+	outfilename = NULL;
+
+/*	Decode "method" to get the actual args to Chaco.
+	Note that "nparts" may correspond to any of several Chaco
+	parameters, depending on the method.  */
+
+	global_method = (int)method[0];
+	local_method = (int)method[1];
+
+	if ((int)method[2] && vwgts) {
+		printf("%s -- Applying weights for %d vertices.\n",
+			   THISFUNCTION,nvtxs);
+		totalvwgt = 0.;
+		for (i=0; i<nvtxs; i++)
+			totalvwgt += vwgts[i];
+	}
+	else {
+		totalvwgt = (double)nvtxs;
+		if      ( (int)method[2] && !vwgts)
+			printf("%s -- method[2]=%d, but no vertex weights specified.\n",
+				   THISFUNCTION,method[2]);
+		else if (!(int)method[2] &&  vwgts)
+			printf("%s -- method[2]=%d, so specified vertex weights ignored.\n",
+				   THISFUNCTION,method[2]);
+	}
+
+	if ((int)method[3] && ewgts) {
+		printf("%s -- Applying weights for %d edges.\n",
+			   THISFUNCTION,start[nvtxs]/2);
+	}
+	else {
+		if      ( (int)method[3] && !ewgts)
+			printf("%s -- method[3]=%d, but no edge weights specified.\n",
+				   THISFUNCTION,method[3]);
+		else if (!(int)method[3] &&  ewgts)
+			printf("%s -- method[3]=%d, so specified edge weights ignored.\n",
+				   THISFUNCTION,method[3]);
+	}
+
+	architecture = (int)method[4];
+	if (global_method == 7) {
+/*		Refine an input partition: "nparts" is the input partition.
+		This seems to work only for hypercube architecture, 
+		so we force a 1-D hypercube with 2-way partitioning.  */
+		architecture = 0;
+		for (i=0; i<nvtxs; i++)
+			assignment[i] = nparts[i];
+		ndims_tot = 1;
+		ndims = 1;
+		nsets = 2;
+	}
+	else if (architecture == 0) {
+/*		Partition for hypercube: "nparts" is # of dimensions (default 1).  */
+		if (!nparts)
+			ndims_tot = 1;
+		else
+			ndims_tot = nparts[0];
+		ndims = (int)method[5];
+		nsets = 2^ndims_tot;
+	}
+	else {
+/*		Partition for mesh: "nparts" is vector of mesh sizes in each
+		dimension, default [2 1 ... 1] with "architecture" dimensions.  */
+		if (!nparts) {
+			mesh_dims[0] = 2;
+			for (i=1; i<architecture; i++)
+				mesh_dims[i] = 1;
+		}
+		else
+			for (i=0; i<architecture; i++)
+				mesh_dims[i] = nparts[i];
+		ndims = (int)method[5];
+		nsets = 1;
+		for (i=0; i<architecture; i++)
+			nsets *= mesh_dims[i];
+	}
+
+/*	Scale goal to total vertex weight so relative goals can be used.  */
+	if (goal) {
+		printf("%s -- Applying goals for %d sets.\n",
+			   THISFUNCTION,nsets);
+		totalgoal = 0.;
+		for (i=0; i<nsets; i++)
+			totalgoal += goal[i];
+		for (i=0; i<nsets; i++)
+			goal[i] *= totalvwgt/totalgoal;
+	}
+
+	vmax = (int)method[6];
+	rqi_flag = (int)method[7];
+	eigtol = method[8];
+	seed = (long)method[9];
+
+	/* Finally, call Chaco */
+
+	printf("\n%s -- Calling Chaco interface:\n\n",
+		   THISFUNCTION);
+	failure = interface(nvtxs, start, adjacency,
+		((int)method[2] && vwgts ? vwgts : NULL),
+		((int)method[3] && ewgts ? ewgts : NULL), x, y, z,
+		outassignname, outfilename, assignment, architecture, ndims_tot, 
+		mesh_dims, goal, global_method, local_method, rqi_flag, vmax, 
+		ndims, eigtol, seed);
+	printf("\n%s -- Chaco interface returning failure=%d.\n",
+		   THISFUNCTION,failure);
+
+/*	Reset adjacency matrix in case calling function needs it.  */
+
+	for (i=0; i<start[nvtxs]; i++)
+		adjacency[i]--;
+
+	return(failure);
+}
+
Index: /issm/trunk/externalpackages/meshpart/install.sh
===================================================================
--- /issm/trunk/externalpackages/meshpart/install.sh	(revision 3996)
+++ /issm/trunk/externalpackages/meshpart/install.sh	(revision 3997)
@@ -35,7 +35,11 @@
 cd ../..
 
-# Build mlchaco
+# Build mlchaco and Chaco
+cp -p Chacox.c src/chaco
+cp -p Chaco.c src/chaco
+cp -p Chaco_m.m src
 cd src/chaco
-make
+make mlchaco
+make Chaco
 # Clean up, specifically the objects left in the chaco directories by mlchaco
 make clean
Index: /issm/trunk/externalpackages/meshpart/meshpart.patch
===================================================================
--- /issm/trunk/externalpackages/meshpart/meshpart.patch	(revision 3996)
+++ /issm/trunk/externalpackages/meshpart/meshpart.patch	(revision 3997)
@@ -1,5 +1,7 @@
-diff -rc src/chaco/Makefile new/chaco/Makefile
+Only in new2/chaco: Chaco.c
+Only in new2/chaco: Chacox.c
+diff -rc src/chaco/Makefile new2/chaco/Makefile
 *** src/chaco/Makefile	1998-07-06 13:08:26.000000000 -0700
---- new/chaco/Makefile	2010-04-06 15:51:33.913224536 -0700
+--- new2/chaco/Makefile	2010-05-27 10:48:50.043321738 -0700
 ***************
 *** 24,51 ****
@@ -23,9 +25,9 @@
   
   
-  MLFILES.c=	${MLCHACO}/mlchaco.c \
-  		${MLCHACO}/bail.c \
-- 		${MLCHACO}/chaco_check_graph.c \
-- 		${MLCHACO}/check_input.c \
--                 ${MLCHACO}/smalloc.c \
+! MLFILES.c=	${MLCHACO}/mlchaco.c \
+! 		${MLCHACO}/bail.c \
+! 		${MLCHACO}/chaco_check_graph.c \
+! 		${MLCHACO}/check_input.c \
+!                 ${MLCHACO}/smalloc.c \
                   ${MLCHACO}/user_params.c
   
@@ -61,10 +63,10 @@
   
   
-  MLFILES.c=	${MLCHACO}/mlchaco.c \
-  		${MLCHACO}/bail.c \
+! MLFILES.c=	${MLCHACO}/Chacox.c \
                   ${MLCHACO}/user_params.c
-+ #		${MLCHACO}/chaco_check_graph.c \
-+ #		${MLCHACO}/check_input.c \
-+ #                ${MLCHACO}/smalloc.c \
++ #		${MLCHACO}/bail.c \  # redefined original with -DMATLAB
++ #		${MLCHACO}/chaco_check_graph.c \  # redefined original with -DMATLAB
++ #		${MLCHACO}/check_input.c \  # redefined original with -DMATLAB
++ #                ${MLCHACO}/smalloc.c \  # redefined original with -DMATLAB
   
   CHFILES.c=	${CHACO}/main/interface.c \
@@ -141,9 +143,10 @@
 ***************
 *** 232,237 ****
---- 246,252 ----
+--- 246,253 ----
   		${CHACO}/util/normalize.c \
   		${CHACO}/util/mergesort.c \
   		${CHACO}/util/randomize.c \
 + 		${CHACO}/util/smalloc.c \
++ 		${CHACO}/util/bail.c \
   		${CHACO}/util/scadd.c \
   		${CHACO}/util/seconds.c \
@@ -151,39 +154,46 @@
 ***************
 *** 243,248 ****
---- 258,265 ----
+--- 259,265 ----
   		${CHACO}/util/vecout.c \
   		${CHACO}/util/vecran.c \
   		${CHACO}/util/vecscale.c 
 + #		${CHACO}/main/user_params.c \
-+ #		${CHACO}/util/bail.c \
   
   MLFILES.o=      $(MLFILES.c:.c=.o)
   
 ***************
-*** 252,259 ****
+*** 252,261 ****
   #		${MATLAB}/bin/cmex CC='gcc -G' -lm ${OFLAGS} ${MLFILES.o} chaco.a; \
   #                 mv mlchaco.mex* ${DEST_DIR}
   
-  mlchaco:	${MLFILES.c} chaco.a Makefile
+! mlchaco:	${MLFILES.c} chaco.a Makefile
 ! 		mex -V4 -output mlchaco ${MLFILES.c} chaco.a -I${CHACO}/main
   		mv mlchaco.mex* ${DEST_DIR}
   
   chaco.a:        ${CHFILES.o}
---- 269,280 ----
+  		${AR} chaco.a ${CHFILES.o} ; ${RANLIB} chaco.a
+  
+--- 269,286 ----
   #		${MATLAB}/bin/cmex CC='gcc -G' -lm ${OFLAGS} ${MLFILES.o} chaco.a; \
   #                 mv mlchaco.mex* ${DEST_DIR}
   
-+ #mlchaco:	${MLFILES.c} chaco.a Makefile
-+ #		mex -V4 -output mlchaco ${MLFILES.c} chaco.a -I${CHACO}/main
-+ #		mv mlchaco.mex* ${DEST_DIR}
+! #mlchaco:	${MLFILES.c} chaco.a Makefile
+! #		mex -V4 -output mlchaco ${MLFILES.c} chaco.a -I${CHACO}/main
+! #		mv mlchaco.mex* ${DEST_DIR}
+! 
+! mlchaco:	${MLFILES.c} ${CHFILES.c} chaco.a Makefile
+! 		${MATLAB}/bin/mex mlchaco.c -largeArrayDims -DMATLAB ${MLFILES.c} chaco.a -I${CHACO}/main
+  		mv mlchaco.mex* ${DEST_DIR}
+  
++ Chaco:	${MLFILES.c} ${CHFILES.c} chaco.a Makefile
++ 		${MATLAB}/bin/mex Chaco.c -largeArrayDims -DMATLAB ${MLFILES.c} chaco.a -I${CHACO}/main
++ 		mv Chaco.mex* ${DEST_DIR}
 + 
-  mlchaco:	${MLFILES.c} chaco.a Makefile
-! 		${MATLAB}/bin/mex -output mlchaco -largeArrayDims -DMATLAB ${MLFILES.c} chaco.a -I${CHACO}/main
-  		mv mlchaco.mex* ${DEST_DIR}
-  
   chaco.a:        ${CHFILES.o}
-diff -rc src/chaco/mlchaco.c new/chaco/mlchaco.c
+  		${AR} chaco.a ${CHFILES.o} ; ${RANLIB} chaco.a
+  
+diff -rc src/chaco/mlchaco.c new2/chaco/mlchaco.c
 *** src/chaco/mlchaco.c	1999-07-23 12:19:09.000000000 -0700
---- new/chaco/mlchaco.c	2010-04-06 15:51:33.920224545 -0700
+--- new2/chaco/mlchaco.c	2010-06-02 11:35:10.327038815 -0700
 ***************
 *** 32,37 ****
@@ -334,7 +344,7 @@
       if (ewgts != NULL) mxFree((char *) ewgts);
       if (x != NULL) mxFree((char *) x);
-diff -rc src/chaco.m new/chaco.m
+diff -rc src/chaco.m new2/chaco.m
 *** src/chaco.m	1999-07-23 12:13:26.000000000 -0700
---- new/chaco.m	2010-04-06 15:51:33.926224552 -0700
+--- new2/chaco.m	2010-06-02 11:35:10.335038823 -0700
 ***************
 *** 142,148 ****
@@ -355,7 +365,9 @@
   else
       vwgts = [];
-diff -rc src/fiedler.m new/fiedler.m
+Only in new2: Chaco.mexa64
+Only in new2: Chaco_m.m
+diff -rc src/fiedler.m new2/fiedler.m
 *** src/fiedler.m	1996-02-23 10:03:20.000000000 -0800
---- new/fiedler.m	2010-04-06 15:51:33.932224559 -0700
+--- new2/fiedler.m	2010-06-02 11:35:10.342038830 -0700
 ***************
 *** 58,64 ****
@@ -379,7 +391,7 @@
       R = chol(L-shift*I);
       Rt = R';
-diff -rc src/meshdemo.m new/meshdemo.m
+diff -rc src/meshdemo.m new2/meshdemo.m
 *** src/meshdemo.m	2002-02-08 07:38:21.000000000 -0800
---- new/meshdemo.m	2010-06-02 11:20:25.564177135 -0700
+--- new2/meshdemo.m	2010-06-02 11:35:10.349038837 -0700
 ***************
 *** 98,103 ****
@@ -414,7 +426,7 @@
   disp(' ');
   disp(' Next is a multilevel method from the "Metis" package.');
-diff -rc src/metis/metismex.c new/metis/metismex.c
+diff -rc src/metis/metismex.c new2/metis/metismex.c
 *** src/metis/metismex.c	2009-07-13 14:37:08.000000000 -0700
---- new/metis/metismex.c	2010-04-06 15:51:33.944224574 -0700
+--- new2/metis/metismex.c	2010-06-02 11:35:10.357038845 -0700
 ***************
 *** 92,98 ****
@@ -436,7 +448,7 @@
   
       /* Find MATLAB's matrix structure */
-diff -rc src/metis/proto.h new/metis/proto.h
+diff -rc src/metis/proto.h new2/metis/proto.h
 *** src/metis/proto.h	2010-04-13 12:32:38.504256000 -0700
---- new/metis/proto.h	2010-04-06 15:51:33.949224580 -0700
+--- new2/metis/proto.h	2010-06-02 11:35:10.406038893 -0700
 ***************
 *** 230,237 ****
@@ -479,2 +491,3 @@
   
 Only in src: metismex.mexsol
+Only in new2: mlchaco.mexa64
Index: /issm/trunk/externalpackages/meshpart/mlchaco_jes_notes.txt
===================================================================
--- /issm/trunk/externalpackages/meshpart/mlchaco_jes_notes.txt	(revision 3996)
+++ /issm/trunk/externalpackages/meshpart/mlchaco_jes_notes.txt	(revision 3997)
@@ -164,2 +164,15 @@
 >               chaco.a -I${CHACO}/main
 
+5/25/10:
+
+- reorganized chaco.m and mlchaco.c drivers into Chaco.c matlab-layer (independent of chaco) and Chacox.c x-layer (independent of matlab).
+
+5/26/10:
+
+- added mexchaco target to Makefile (and related changes).
+
+5/27/10
+
+- added define for exit(status) to ${CHACO}/main/defs.h so that local bail.c is
+  unnecessary.
+