/*****************************************************************************/
/*                                                                           */
/*  Count Me                                                                 */
/*                                                                           */
/*  Parceled Mesh Analyzer                                                   */
/*  (countme.c)                                                              */
/*                                                                           */
/*  Version 0.9 (prerelease)                                                 */
/*  Feb 17, 1998                                                             */
/*                                                                           */
/*  (c) copyright 1998                                                       */
/*  David R. O'Hallaron                                                      */
/*  School of Computer Science                                               */
/*  Carnegie Mellon University                                               */
/*  5000 Forbes Avenue                                                       */
/*  Pittsburgh, Pennsylvania  15213-3891                                     */
/*  droh@cs.cmu.edu                                                          */
/*                                                                           */
/*  This program may be freely redistributed under the condition that the    */
/*    copyright notices (including this entire header and the copyright      */
/*    notice printed when the `-h' switch is selected) are not removed, and  */
/*    no compensation is received.  Private, research, and institutional     */
/*    use is free.  You may distribute modified versions of this code UNDER  */
/*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
/*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
/*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
/*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
/*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
/*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
/*    customer, and you are instead telling them how they can obtain it for  */
/*    free, then you are not required to make any arrangement with me.)      */
/*                                                                           */
/*  This program was created as part of the Archimedes project in the School */
/*    of Computer Science at Carnegie Mellon University.                     */
/*                                                                           */
/*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
/*    whatsoever.  Use at your own risk.                                     */
/*                                                                           */
/*****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAX_NAME 128 
#define MAX_HIST  64
#define MAX_SUBS 512 /* maximum subdomains */

#define BIGINT 1<<30
#define BIGFLOAT 1e30

/* 
 * various data structures 
 */
struct stats {
    int min;
    int max;
    int total;
    float avg;
    int hist[MAX_HIST];
    int histsize;
};

struct node {
    int id;
    int subdomain;
    int interface;
    struct node *nextnp;
};

/* 
 * program name
 */
char progname[] = "countme";

/*
 * global information
 */
struct gi {
    char meshname[MAX_NAME];
    char packfilename[MAX_NAME];

    /* command line options */
    int quiet;
    int help;

    /* global mesh info */
    int globalnodes;          /* number of nodes */
    int globalelems;          /* number of elements */
    int globaledges;          /* number of edges */
    int globalsharednodes;    /* total number of interface nodes */
    int globalnonsharednodes; /* total number of interior nodes */
    int globalsharededges;    /* total number of interface edges */
    int globalnonsharededges; /* total number of interior edges */
    int meshdim;              /* mesh dimension */
    int elemsize;             /* nodes per elem */
    int globalmsgs;           /* communication volume in msgs */
    int globalmsgwords;       /* communication volume in words per dof */
    int subdomains;           /* number of subdomains */
    int globalmatrixentries;  /* total number of nonzeros */
    float beta;               /* max msg/vol balance ratio */

    /* various max/min/avg/hist statistics */
    struct stats degreestats;            /* degree per node */
    struct stats *nzstats;               /* nonzeros per row per subdomain */
    struct stats nodeallstats;           /* total nodes per subdomain */
    struct stats nodesharedstats;        /* shared nodes per subdomain */
    struct stats elemstats;              /* elements per subdomain */
    struct stats edgestats;              /* edges per subdomain */
    struct stats sharededgestats;        /* shared edges per subdomain */
    struct stats localmsgsstats;         /* msgs per subdomain */
    struct stats localmsgwordsstats;     /* msg words per subdomain */
    struct stats *localpermsgwordsstats; /* words per msg per subdomain */
    struct stats globalpermsgwordsstats; /* msg words over all messages */
    struct stats symentriesstats;        /* matrix entries per sub (sym) */
    struct stats nonsymentriesstats;     /* mat entries per sub (nonsym) */

    /* supporting mesh info */
    int maxlocalnodes;   /* max local nodes in any subdomain */
    int dbisectionwords; /* default bisection volume */
    int rbisectionwords; /* random bisection volume */
    int *nodeall;        /* total nodes per subdomain */
    int *nodemine;       /* shared nodes owned by me per subdomain */
    int *nodepriv;       /* non-shared nodes per subdomain */
    int *globalnode[MAX_SUBS];   /* local to global node mapping */
    int *elems;          /* elements per subdomain */
    int *refcnt;         /* node reference count to determine shared nodes */
    int *commeach[MAX_SUBS]; /* i shares commeach[i][j] nodes with j */

    /* mesh adjacency list */
    struct node *nodevec;
} gi;

/* 
 * private function prototypes 
 */
void parsecommandline(int argc, char **argv, struct gi *gip);
void syntax(void);
void info(void);
void readheader(FILE *packfile, struct gi *gip);
void readnodes(FILE *packfile, struct gi *gip);
void readelems(FILE *packfile, struct gi *gip);
void readmatrix(FILE *packfile, struct gi *gip);
void readcomm(FILE *packfile, struct gi *gip);
void emit(struct gi *gip);
void dostats0(int *vec, int n, struct stats *statsp);
int dohist0(int *vec, int n, int *hist, int maxhist);
void dostats1(int *vec, int n, struct stats *statsp);
int dohist1(int *vec, int n, int *hist, int maxhist);
void printstats(struct stats *statsp);
void insertedge(struct node *nodevec, int node1, int node2, int subdomain);
struct node *findedge(struct node *nodevec, int node1, int node2);
int numedges(struct node *np);

void main(int argc, char **argv) {
    struct gi *gip = &gi;
    FILE *packfile;

    parsecommandline(argc, argv, gip);

    if (!(packfile = fopen(gip->packfilename, "r"))) {
	fprintf(stderr, "%s: Can't open %s\n", progname, gip->packfilename);
	exit(0);
    }

    if (!gip->quiet) {
	fprintf(stderr, "Reading pack file (%s).\n", gip->packfilename);
	fflush(stderr);
    }

    readheader(packfile, gip);
    readnodes(packfile, gip);
    readelems(packfile, gip);
    readmatrix(packfile, gip);
    readcomm(packfile, gip);
    fclose(packfile);

    emit(gip);
    
    exit(0);
}

void readheader(FILE *packfile, struct gi *gip)
{
    fscanf(packfile, "%d %d\n", &gip->globalnodes, &gip->meshdim);
    fscanf(packfile, "%d %d\n", &gip->globalelems, &gip->elemsize);
    fscanf(packfile, "%d\n", &gip->subdomains);
    fscanf(packfile, "%*d\n");
}

void readnodes(FILE *packfile, struct gi *gip) 
{
    int i, j, k;
    int *nodeshared;

    if (!gip->quiet) {
	fprintf(stderr, "  Reading nodes.\n");
	fflush(stderr);
    }

    if (!(gip->nodeall = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nodeall\n", progname);
	exit(0);
    }

    if (!(gip->nodemine = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nodemine\n", progname);
	exit(0);
    }

    if (!(gip->nodepriv = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nodepriv\n", progname);
	exit(0);
    }

    if (!(gip->refcnt = (int *) malloc((gip->globalnodes + 1) * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate refcnt\n", progname);
	exit(0);
    }

    gip->maxlocalnodes = 0;
    for (i = 0; i < gip->subdomains; i++) {
	fscanf(packfile, "%d %d %d\n",
	       &gip->nodeall[i], &gip->nodemine[i], &gip->nodepriv[i]);
	if (gip->nodeall[i] > gip->maxlocalnodes)
	  gip->maxlocalnodes = gip->nodeall[i];
    }

    for (i = 0; i < gip->subdomains; i++) {
	if (!(gip->globalnode[i] = 
	      (int *) malloc(gip->maxlocalnodes * sizeof(int)))) {
	    fprintf(stderr, "%s: couldn't allocate globalnode[%d]\n", 
		    progname, i);
	}
    }

    for (i = 1; i <= gip->globalnodes; i++) {
	gip->refcnt[i] = 0;
    }

    for (i = 0; i < gip->subdomains; i++) {
	for (j = 0; j < gip->nodeall[i]; j++) {
	    fscanf(packfile, "%d", &gip->globalnode[i][j]);
	    gip->refcnt[gip->globalnode[i][j]]++;
	    for (k = 0; k < gip->meshdim; k++) {
		fscanf(packfile, "%*f");
	    }
	}
    }

    dostats0(gip->nodeall, gip->subdomains, &gip->nodeallstats);

    if (!(nodeshared = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nodeshared\n", progname);
	exit(0);
    }

    for (i = 0; i < gip->subdomains; i++) {
	nodeshared[i] = gip->nodeall[i] - gip->nodepriv[i];
    }
    dostats0(nodeshared, gip->subdomains, &gip->nodesharedstats);

    gip->globalsharednodes = gip->globalnonsharednodes = 0;
    for (i = 1; i <= gip->globalnodes; i++) {
	if (gip->refcnt[i] > 1) {
	    gip->globalsharednodes++;
	}
	else if (gip->refcnt[i] == 1) {
	    gip->globalnonsharednodes++;
	}
	else {
	    fprintf(stderr, "%s: bogus reference count\n", progname);
	    exit(0);
	}
    }
    if ((gip->globalsharednodes + gip->globalnonsharednodes) != gip->globalnodes) {
	fprintf(stderr, "%s: inconsistent shared node count\n", progname);
	exit(0);
    }

    (void)free(nodeshared);
}

void readelems(FILE *packfile, struct gi *gip) {
    int i, j, k;
    int node;
    int *gn;
    int *degree;

    if (!gip->quiet) {
	fprintf(stderr, "  Reading elements.\n");
	fflush(stderr);
    }

    if (!(gn = (int *) malloc(gip->elemsize * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate little gn\n", progname);
	exit(0);
    }

    if (!(degree = (int *) malloc((gip->globalnodes+1) * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate degree\n", progname);
	exit(0);
    }

    if (!(gip->elems = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate elems\n", progname);
	exit(0);
    }

    if ((gip->nodevec = (struct node *) 
	 malloc((gip->globalnodes+1) * sizeof(struct node))) == NULL) {
	fprintf(stderr, "%s: couldn't allocate nodevec\n", progname);
	exit(0);
    }

    for (i=1; i < gip->globalnodes+1; i++) {
	gip->nodevec[i].id = i;
	gip->nodevec[i].nextnp = NULL;
	degree[i] = 0;
    }

    for (i = 0; i < gip->subdomains; i++) {
	fscanf(packfile, "%d\n", &gip->elems[i]);
    }
    
    if (gip->elemsize != 3 && gip->elemsize != 4) {
	fprintf(stderr, "Error: elements with %d corners are not supported\n",
		gip->elemsize);
	exit(0);
    }

    if (gip->meshdim == 2 && gip->elemsize == 6) {
	fprintf(stderr, "Warning: edge counts for 2D quadratic elements are not supported\n");
	fflush(stderr);
    }
    
    for (i = 0; i < gip->subdomains; i++) {
	for (j = 0; j < gip->elems[i]; j++) {
	    fscanf(packfile, "%*d");
	    for (k = 0; k < gip->elemsize; k++) {
		fscanf(packfile, "%d", &node);
		gn[k] = gip->globalnode[i][node];
	    }

	    /* 3d linear tetrahedra */
	    if (gip->elemsize == 4) {
		insertedge(gip->nodevec, gn[0], gn[1], i);
		insertedge(gip->nodevec, gn[1], gn[0], i);
		insertedge(gip->nodevec, gn[0], gn[2], i);
		insertedge(gip->nodevec, gn[2], gn[0], i);
		insertedge(gip->nodevec, gn[0], gn[3], i);
		insertedge(gip->nodevec, gn[3], gn[0], i);

		insertedge(gip->nodevec, gn[1], gn[2], i);
		insertedge(gip->nodevec, gn[2], gn[1], i);
		insertedge(gip->nodevec, gn[1], gn[3], i);
		insertedge(gip->nodevec, gn[3], gn[1], i);

		insertedge(gip->nodevec, gn[2], gn[3], i);
		insertedge(gip->nodevec, gn[3], gn[2], i);
	    }

	    /* 2d linear triangles */
	    else if (gip->elemsize == 3) {
		insertedge(gip->nodevec, gn[0], gn[1], i);
		insertedge(gip->nodevec, gn[1], gn[0], i);
		insertedge(gip->nodevec, gn[0], gn[2], i);
		insertedge(gip->nodevec, gn[2], gn[0], i);
		insertedge(gip->nodevec, gn[1], gn[2], i);
		insertedge(gip->nodevec, gn[2], gn[1], i);
	    }
	    
	    /* everything else */
	    else 
	      ;
	}
    }
    
    gip->globaledges = 0;
    for (i = 1; i < gip->globalnodes+1; i++) {
	degree[i] = numedges(&gip->nodevec[i]);
	gip->globaledges += degree[i];
    }
    gip->globaledges /= 2; /* each edge appears twice */
    dostats1(degree, gip->globalnodes, &gip->degreestats);
    dostats0(gip->elems, gip->subdomains, &gip->elemstats);

    (void)free(degree);
    (void)free(gn);
}  

void readmatrix(FILE *packfile, struct gi *gip) 
{
    int i, j;
    int node1, node2;
    int *nonzerosperrow;
    int *symentries;
    int *nonsymentries;
    int *edges;
    int *sharededges;
    struct node *np;

    if (!gip->quiet) {
	fprintf(stderr, "  Reading matrix.\n");
	fflush(stderr);
    }

    if (!(edges = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate edges\n", progname);
	exit(0);
    }

    if (!(sharededges = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate sharededges\n", progname);
	exit(0);
    }

    if (!(symentries = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate symentries\n", progname);
	exit(0);
    }

    if (!(nonsymentries = (int *) 
	  malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nonsymentries\n", 
		progname);
	exit(0);
    }

    if (!(nonzerosperrow = (int *)
	  malloc(gip->maxlocalnodes * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate nonzerosper\n", progname);
	exit(0);
    }

    if (!(gip->nzstats = (struct stats *)
	  malloc(gip->subdomains * sizeof(struct stats)))) {
	fprintf(stderr, "%s: couldn't allocate nzstats\n", progname);
	exit(0);
    }

    for (i = 0; i < gip->subdomains; i++) { 
	fscanf(packfile, "%d %*d\n", &symentries[i]); 
    }

    for (i = 0; i < gip->subdomains; i++) {
	for (j = 0; j < gip->maxlocalnodes; j++)
	  nonzerosperrow[j] = 0;

	edges[i] = sharededges[i] = 0;
	for (j = 0 ; j < symentries[i]; j++) {
	    fscanf(packfile, "%d %d", &node1, &node2);

	    /* record the number of edges and interface edges per subdomain */
	    if (node1 != node2) {
		edges[i]++;
		np = findedge(gip->nodevec, 
			      gip->globalnode[i][node1], 
			      gip->globalnode[i][node2]);
		if (np && np->interface) {
		    sharededges[i]++;
		}
	    }

	    /* 
	     * Record the number of nonzeros per row assuming nonsymmetric
	     * storage.
	     */
	    if (node1 == node2) {
		nonzerosperrow[node1]++;
	    }
	    else {
		nonzerosperrow[node1]++;
		nonzerosperrow[node2]++;
	    }
	}
	dostats0(nonzerosperrow, gip->nodeall[i], &gip->nzstats[i]);
    }   
    dostats0(symentries, gip->subdomains, &gip->symentriesstats);
    for (i = 0; i < gip->subdomains; i++) {
	nonsymentries[i] = symentries[i] * 2 - gip->nodeall[i];
    }
    dostats0(nonsymentries, gip->subdomains, &gip->nonsymentriesstats);
    dostats0(edges, gip->subdomains, &gip->edgestats);
    dostats0(sharededges, gip->subdomains, &gip->sharededgestats);

    /* 
     * count the number of global interface edges 
     */
    gip->globalsharededges = gip->globalnonsharededges = 0;
    for (i = 1; i < gip->globalnodes+1; i++) {
	np = gip->nodevec[i].nextnp;
	while (np) {
	    if (np->interface) {
		gip->globalsharededges++;
	    }
	    else {
		gip->globalnonsharededges++;
	    }
	    np = np->nextnp;
	}
    }
    gip->globalsharededges /= 2;
    gip->globalnonsharededges /= 2;
    if ((gip->globalsharededges + gip->globalnonsharededges) !=
	gip->globaledges) {
	fprintf(stderr, "%s: inconsistent shared edge count: %d + %d (%d)\n",
		progname, gip->globalsharededges, gip->globalnonsharededges, 
		gip->globaledges);
	exit(0);
    }

    (void)free(nonzerosperrow);
    (void)free(symentries);
    (void)free(nonsymentries);
    (void)free(edges);
    (void)free(sharededges);
}

void readcomm(FILE *packfile, struct gi *gip) {
    int i, j, k;
    int rn, tmp, sched[MAX_SUBS];
    int globalcount;
    int localcount;
    int *commall;
    int *localmsgs;
    int *localmsgwords;
    int *localpermsgwords;
    int *globalpermsgwords;
    float cmax, bmax, bi, ci, term1, term2, term;
    
    if (!gip->quiet) {
	fprintf(stderr, "  Reading communication.\n");
	fflush(stderr);
    }

    if (!(gip->localpermsgwordsstats = (struct stats *)
	  malloc(gip->subdomains * sizeof(struct stats)))) {
	fprintf(stderr, "%s: couldn't allocate localpermsgwordsstats\n", 
		progname);
	exit(0);
    }

    if (!(commall = (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate commall\n", progname);
	exit(0);
    }

    if (!(localmsgs = 
	  (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate localmsgs\n", progname);
	exit(0);
    }

    if (!(localmsgwords = 
	  (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate localmsgwords\n", progname);
	exit(0);
    }

    if (!(localpermsgwords = 
	  (int *) malloc(gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate localpermsgwords\n", progname);
	exit(0);
    }

    if (!(globalpermsgwords = 
	  (int *) malloc(gip->subdomains * gip->subdomains * sizeof(int)))) {
	fprintf(stderr, "%s: couldn't allocate globalpermsgwords\n", progname);
	exit(0);
    }

    for (i = 0; i < gip->subdomains; i++) {
	fscanf(packfile, "%d\n", &commall[i]);
    }
    
    for (i = 0; i < gip->subdomains; i++) {
	if (!(gip->commeach[i] = 
	      (int *) malloc(gip->subdomains * sizeof(int)))) {
	    fprintf(stderr, "%s: couldn't allocate commeach[%d]\n", 
		    progname, i);
	    exit(0);
	}
    }

    gip->globalmsgs = 0;
    gip->globalmsgwords = 0;
    globalcount = 0;
    for (i = 0; i < gip->subdomains; i++) {
	localmsgs[i] = 0;
	localmsgwords[i] = 0;
	for (j = 0; j < gip->subdomains; j++) {
	    fscanf(packfile, "%d", &gip->commeach[i][j]);
	    gip->globalmsgwords += gip->commeach[i][j];
	    localmsgwords[i] += (gip->commeach[i][j] * 2); 
	    if (gip->commeach[i][j] > 0) {
		gip->globalmsgs++;
		localmsgs[i] += 2;
	    }
	    for (k = 0; k < gip->commeach[i][j]; k++) {
		fscanf(packfile, "%*d");
	    }
	}
	localcount = 0;
	for (j = 0; j < gip->subdomains; j++) {
	    if (gip->commeach[i][j] > 0) {
		localpermsgwords[localcount++] = gip->commeach[i][j]; 
		localpermsgwords[localcount++] = gip->commeach[i][j]; 
		globalpermsgwords[globalcount++] = gip->commeach[i][j];
	    }
	} 
	if (localcount != (localmsgs[i])) {
	    fprintf(stderr, "%s: inconsistent local message count\n", 
		    progname);
	    exit(0);
	}
	dostats0(localpermsgwords, localcount, &gip->localpermsgwordsstats[i]);
    }

    if (gip->globalmsgs != globalcount) {
	fprintf(stderr, "%s: inconsistent global message count\n", progname);
	exit(0);
    }

    dostats0(globalpermsgwords, globalcount, &gip->globalpermsgwordsstats);
    dostats0(localmsgs, gip->subdomains, &gip->localmsgsstats);
    dostats0(localmsgwords, gip->subdomains, &gip->localmsgwordsstats);

    /* compute the default bisection volume */
    gip->dbisectionwords = 0;
    for (i=0; i<gip->subdomains/2; i++) {
      for (j=gip->subdomains/2; j<gip->subdomains; j++) {
	gip->dbisectionwords += gip->commeach[i][j];
      }
    }
    gip->dbisectionwords *= 2;

    /* make a random permutation of the subdomain names */
    for (i=0; i<gip->subdomains; i++)
      sched[i] = i;
    srand(1); /* 1 ensures that sequence is repeatable */
    for (i=gip->subdomains-1; i>0; i--) {
      rn = rand()%(i+1);
      tmp = sched[i];
      sched[i] = sched[rn];
      sched[rn] = tmp;
    }

    /* compute a random bisection volume */
    gip->rbisectionwords = 0;
    for (i=0; i<gip->subdomains/2; i++) {
      for (j=gip->subdomains/2; j<gip->subdomains; j++) {
	gip->rbisectionwords += gip->commeach[sched[i]][sched[j]];
      }
    }
    gip->rbisectionwords *= 2;

    /* compute balance metric */
    if (gip->subdomains > 1) {
      gip->beta = BIGFLOAT;
      for (i=0; i<gip->subdomains; i++) {
	cmax = gip->localmsgwordsstats.max;
	bmax = gip->localmsgsstats.max;
	bi = localmsgs[i];
	ci = localmsgwords[i]; 
	term1 = ((bmax-bi)*cmax)/(bmax*ci);
	term2 = ((cmax-ci)*bmax)/(bi*cmax);
	term = (term1 > term2) ? term1 : term2;
	if (term < gip->beta) 
	  gip->beta = term;
      }
      gip->beta = 1 + gip->beta;
    }

    (void)free(localmsgs); 
    (void)free(localmsgwords);
    (void)free(commall);
    (void)free(globalpermsgwords);
}    

void emit(struct gi *gip)
{
    int i;
    int mini, maxi; 
    int num0;
    char num1;
    float minavg, maxavg; 

    num0 = num1 = 1;
    printf("%1d. Global mesh summary (%s).\n", num0, gip->packfilename);
    printf("   Nodes    : %9d  Edges    : %9d  Elements: %d\n",
	   gip->globalnodes, gip->globaledges, gip->globalelems);
    printf("   Interface: %9d  Interface: %9d\n", 
	   gip->globalsharednodes, gip->globalsharededges);
    printf("   Interior : %9d  Interior : %9d\n", 
	   gip->globalnonsharednodes, gip->globalnonsharededges);
    printf("\n");
    printf("   Dimensions: %d  Subdomains: %d  Nodes/elem: %d\n", 
	   gip->meshdim, gip->subdomains, gip->elemsize);
    printf("\n");

    printf("%d.%d. Node degrees.", num0, num1++);
    printstats(&gip->degreestats);
    printf("\n");

    printf("%d.%d. Nodes per subdomain.", num0, num1++);
    printstats(&gip->nodeallstats);
    printf("\n");

    printf("%d.%d. Interface nodes per subdomain.", num0, num1++);
    printstats(&gip->nodesharedstats);
    printf("\n");

    printf("%d.%d. Elements per subdomain.", num0, num1++);
    printstats(&gip->elemstats);
    printf("\n");

    printf("%d.%d. Edges per subdomain.", num0, num1++);
    printstats(&gip->edgestats);
    printf("\n");

    printf("%d.%d. Interface edges per subdomain.", num0, num1++);
    printstats(&gip->sharededgestats);
    printf("\n");

    printf("%1d. Global matrix summary.\n", ++num0);
    printf("   Nonzero matrix entries (symmetric): %d\n", 
	   gip->globaledges + gip->globalnodes);
    printf("   Nonzero matrix entries (nonsymmetric): %d\n", 
	   gip->globaledges * 2 + gip->globalnodes);
    printf("\n");
    printf("   (Note: Each matrix entry consists of dof*dof words, where dof\n");
    printf("    is the degrees of freedom in the simulation.)\n");
    printf("\n");

    num1 = 1;

    printf("%d.%d. Nonzero matrix entries per subdomain (symmetric).", 
	   num0, num1++);
    printstats(&gip->symentriesstats);
    printf("\n");

    printf("%d.%d. Nonzero matrix entries per subdomain (nonsymmetric).", 
	   num0, num1++);
    printstats(&gip->nonsymentriesstats);
    printf("\n");

    minavg = BIGFLOAT;
    maxavg = 0.0;
    for (i = 0; i < gip->subdomains; i++) {
	if (gip->nzstats[i].avg < minavg) {
	    mini = i;
	    minavg = gip->nzstats[i].avg;
	}
	if (gip->nzstats[i].avg > maxavg) {
	    maxi = i;
	    maxavg = gip->nzstats[i].avg;
	}
    }
    printf("%d.%d. Fewest average entries per row (non-symmetric) [subdomain %d].", num0, num1++, mini);
    printstats(&gip->nzstats[mini]);
    printf("\n");

    printf("%d.%d. Most average entries per row (non-symmetric) [subdomain %d]."
	   , num0, num1++, maxi);
    printstats(&gip->nzstats[maxi]);
    printf("\n");

    if (gip->subdomains == 1)
      return;

    printf("%d. Global communication summary.\n", ++num0);
    printf("   Messages: %d  Volume (words): %d  Beta: %.2f\n", 
	   gip->globalmsgs, gip->globalmsgwords, gip->beta);
    printf("   Default/Random bisection volumes: %d/%d [%2.0f%%/%2.0f%%]\n", 
	   gip->dbisectionwords, gip->rbisectionwords, 
	   ((float)gip->dbisectionwords/(float)gip->globalmsgwords)*100.0,
	   ((float)gip->rbisectionwords/(float)gip->globalmsgwords)*100.0);

    printf("\n");
    printf("   (Note: All message sizes and communication volumes in this\n");
    printf("    section assume one degree of freedom. For a particular\n"); 
    printf("    simulation, multiply by the degrees of freedom to get the\n");
    printf("    actual number of words transferred during an assembly step.)\n");
    printf("\n");

    num1 = 1;
    printf("%d.%d. Messages per subdomain.", num0, num1++);
    printstats(&gip->localmsgsstats);
    printf("\n");

    printf("%d.%d. Communication volume (words) per subdomain.", 
	   num0, num1++);
    printstats(&gip->localmsgwordsstats);
    printf("\n");

    printf("%d.%d. Message sizes (words).", 
	   num0, num1++);
    printstats(&gip->globalpermsgwordsstats);
    printf("\n");

    minavg = BIGFLOAT;
    maxavg = 0.0;
    for (i = 0; i < gip->subdomains; i++) {
	if (gip->localpermsgwordsstats[i].avg < minavg) {
	    mini = i;
	    minavg = gip->localpermsgwordsstats[i].avg;
	}
	if (gip->localpermsgwordsstats[i].avg > maxavg) {
	    maxi = i;
	    maxavg = gip->localpermsgwordsstats[i].avg;
	}
    }
    printf("%d.%d. Smallest average message size (words) per subdomain [subdomain %d].", num0, num1++, mini);
    printstats(&gip->localpermsgwordsstats[mini]);
    printf("\n");
    printf("%d.%d. Largest average message size (words) per subdomain [subdomain %d].", num0, num1++, maxi);
    printstats(&gip->localpermsgwordsstats[maxi]);
    printf("\n");
}

void parsecommandline(int argc, char **argv, struct gi *gip) 
{
    int i, j;

    gip->quiet = gip->help = 0;
    gip->meshname[0] = '\0';
    for (i = 1; i < argc; i++) {
	if (argv[i][0] == '-') {
	    for (j = 1; argv[i][j] != '\0'; j++) {
		if (argv[i][j] == 'Q') {
		    gip->quiet = 1;
		}
		if ((argv[i][j] == 'h') || (argv[i][j] == 'H') ||
		    (argv[i][j] == '?')) {
		    info();
		    exit(0);
		}
	    }
	} else {
	    strcpy(gip->meshname, argv[i]);
	}
    }
    if (gip->meshname[0] == '\0') {
	syntax();
	exit(0);
    }

    if (!strcmp(&gip->meshname[strlen(gip->meshname) - 5], ".pack")) {
	gip->meshname[strlen(gip->meshname) - 5] = '\0';
    }
    strcpy(gip->packfilename, gip->meshname);
    strcat(gip->packfilename, ".pack");
}
 
void dostats1(int *vec, int n, struct stats *statsp) {
    int i;
    
    if (n == 0) {
	statsp->min = 0;
	statsp->max = 0; 
	statsp->avg = 0.0;
	statsp->histsize = dohist1(vec, n, statsp->hist, MAX_HIST);
	return;
    }

    statsp->min = BIGINT;
    statsp->max = 0; 
    statsp->avg = 0.0;

    for (i = 1; i <= n; i++) {
	statsp->avg += vec[i];
	statsp->min = (vec[i] < statsp->min) ? vec[i] : statsp->min;
	statsp->max = (vec[i] > statsp->max) ? vec[i] : statsp->max;
    }
    statsp->total = statsp->avg;
    statsp->avg /= n;
    statsp->histsize = dohist1(vec, n, statsp->hist, MAX_HIST);
}

int dohist1(int *vec, int n, int *histogram, int maxhist)   
{
    int i, maxbin, count, bin;

    for (i = 0; i < maxhist; i++) {
	histogram[i] = 0;
    }
    maxbin = 0;
    for (i = 1; i <= n; i++) {
	bin = 0;
	count = 1;
	while (count < vec[i]) {
	    bin++;
	    count <<= 1;
	}
	if (bin > (maxhist-1)) {
	    fprintf(stderr, "%s: exceeded histogram limit\n", progname);
	    exit(0);
	}
	histogram[bin]++;
	maxbin = (bin > maxbin) ? bin : maxbin;
    }
    return(maxbin+1);
}

void dostats0(int *vec, int n, struct stats *statsp) {
    int i;
    
    if (n == 0) {
	statsp->min = 0;
	statsp->max = 0; 
	statsp->avg = 0.0;
	statsp->histsize = dohist0(vec, n, statsp->hist, MAX_HIST);
	return;
    }

    statsp->min = BIGINT;
    statsp->max = 0; 
    statsp->avg = 0.0;

    for (i = 0; i < n; i++) {
	statsp->avg += vec[i];
	statsp->min = (vec[i] < statsp->min) ? vec[i] : statsp->min;
	statsp->max = (vec[i] > statsp->max) ? vec[i] : statsp->max;
    }
    statsp->total = statsp->avg;
    statsp->avg /= n;
    statsp->histsize = dohist0(vec, n, statsp->hist, MAX_HIST);
}

int dohist0(int *vec, int n, int *histogram, int maxhist)   
{
    int i, maxbin, count, bin;

    for (i = 0; i < maxhist; i++) {
	histogram[i] = 0;
    }
    maxbin = 0;
    for (i = 0; i < n; i++) {
	if (vec[i] > 0) {
	    bin = 0;
	    count = 1;
	    while (count < vec[i]) {
		bin++;
		count <<= 1;
	    }
	    if (bin > (maxhist-1)) {
		fprintf(stderr, "%s: exceeded histogram limit\n", progname);
		exit(0);
	    }
	    histogram[bin]++;
	    maxbin = (bin > maxbin) ? bin : maxbin;
	}
    }
    return(maxbin+1);
}

void printstats(struct stats *statsp) 
{
    int i;
    int count=0;
    int total=0;

    printf("\n");
    printf("       Range          Count          Min: %d Avg: %.1f Max: %d Tot: %d\n", statsp->min, statsp->avg, statsp->max, statsp->total);

    while (!statsp->hist[count])
      count++;
    for (i = count; i < statsp->histsize; i++) {
	total += statsp->hist[i];
	printf("%8d - %-8d : %d\n", ((1 << i) >> 1) + 1, 1 << i,
	       statsp->hist[i]);
    }
    printf("       Total          %-10d\n", total);
}

int numedges(struct node *np) {
    int count=0;

    while (np->nextnp) {
	count++;
	np = np->nextnp;
    }
    return count;
}

void insertedge(struct node *nodevec, int node1, int node2, int subdomain) {
    struct node *currnp, *np;

    currnp = &nodevec[node1];
    while (currnp->nextnp && (node2 > currnp->nextnp->id)) {
	currnp = currnp->nextnp;
    }

    /* duplicate edge */
    if (currnp->nextnp && (node2 == currnp->nextnp->id)) { 
	if (!(currnp->nextnp->interface) && 
	    (subdomain != currnp->nextnp->subdomain)) { /* interface edge */
	    currnp->nextnp->interface = 1;
	}
	return;
    }

    if ((np = (struct node *) malloc(sizeof(struct node))) == NULL) {
	fprintf(stderr, "%s: coudn't allocate node\n", progname);
	exit(0);
    }

    np->id = node2;
    np->subdomain = subdomain;
    np->interface = 0;
    np->nextnp = currnp->nextnp;
    currnp->nextnp = np;
}

struct node *findedge(struct node *nodevec, int node1, int node2) {
    struct node *currnp;

    currnp = &nodevec[node1];
    while (currnp->nextnp && (node2 != currnp->nextnp->id)) {
	currnp = currnp->nextnp;
    }
    return currnp->nextnp;
}

void info(void)
{
  printf("Count Me (parceled mesh analyzer).\n");
  printf("Version 0.9 (prerelease)\n\n");
  printf(
"(c) Copyright January, 1996, David O'Hallaron (droh@cs.cmu.edu).\n");
  printf(
"Carnegie Mellon University, Pittsburgh, Pennsylvania, 15213-3891.\n");
  printf(
"Created as part of the Archimedes project.\n\n");
  printf(
"Count Me reads a parceled mesh (in a .pack file) and prints a variety of\n");
  printf(
"statistics. The statistics summarize structural properties of the mesh such as\n");
 printf(
"the number of nodes, edges, elements, subdomains, element size, and dimension,\n");
  printf(
"the sizes of the global and local sparse matrices induced by the partitioned\n");
  printf(
"mesh, and the volume of communication between subdomains. For each statistic\n");
  printf(
"that consists of a set values, Count Me reports the min, max, and average\n");
  printf(
"values, as well as a histogram. The command syntax is:\n\n");
  printf("countme [-Qh] meshname[.pack]\n\n");
  printf("Command line options:\n\n");
  printf(
"    -Q  Quietly suppresses all explanation of what this program is doing\n");
  printf("        unless an error occurs.\n");
  printf(
"    -h  Prints this message\n");
}

void syntax(void)
{
    printf("Usage: countme [-hQ] meshname[.pack]\n");
}