/***************************************************************************/
/* renumber - program to convert an n-pack file into a 1-pack file.        */
/*           Girija Narlikar and David O'Hallaron, February, 1998,         */
/*           Carnegie Mellon University.                                   */
/*                                                                         */
/* usage: renumber old-1-packfile p-packfile > new-1-packfile              */
/*                                                                         */  
/* The program renumbers the nodes in a 1-subdomain packfile according     */
/* to the partition specified in a p-subdomain packfile. The purpose is to */
/* decrease the bandwidth of the global matrix and provide better locality */
/* of reference.                                                           */
/*                                                                         */
/* Reorder first reads the node order from the p-packfile and stores       */
/* it in an array called "permutation".  It then reads the old 1-packfile, */
/* relabels the (row,col) tuples according to this permutation, converts   */
/* to upper triangular tuples, sorts the new tuples, and writes            */
/* out the new matrix to standard output.                                  */
/*                                                                         */
/* Distributed as part of the Spark98 Kernels                              */
/* Copyright (c) David O'Hallaron 1998                                     */
/*                                                                         */
/* You are free to use this software without restriction. If you find that */
/* the Spark kernels are helpful to you, it would be very helpful if you   */
/* sent me at droh@cs.cmu.edu letting me know how you are using them.      */
/***************************************************************************/ 

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

/* sparse matrix tuple */
typedef struct {
  int row;
  int col;
} rowcol_t;

/* 
 * program wide constants 
 */
#define STRLEN 128        /* default string length */
#define DOF 3             /* degrees of freedom in underlying simulation */
#define EMPTY -9999       /* uninitialized variable */

/*
 * global variables
 */
struct gi { 
  /* misc */
  char progname[STRLEN];     /* program name */
  FILE *packpfp;             /* file descriptor p-subdomain pack file */
  FILE *pack1fp;             /* file descriptor 1-subdomain pack file */

  /* command line options */
  int quiet;
  char packpfilename[STRLEN];/* p-pack file name */
  char pack1filename[STRLEN];/* 1-pack file name */

  /* 1-packfile problem parameters */
  int globalnodes1;          /* number of global nodes */
  int globalelems1;          /* number of global elements */
  int mesh_dim1;             /* mesh dimension */
  int corners1;              /* nodes per element */
  int subdomains1;           /* number of partition sets */ 
  int processors1;           /* not used */

  int elems1;                /* number of elements on each PE */
  int matrixlen1;            /* number of nonzeros on each PE */

  int *globalnode1;          /* local-to-global node mapping function */
  double (*coord1)[DOF];     /* geometric node coordinates */
  int (*vertex1)[4];         /* nodes associated with each element */

  /* p-packfile problem parameters */
  int globalnodesp;          /* number of global nodes */
  int globalelemsp;          /* number of global elements */
  int mesh_dimp;             /* mesh dimension */
  int cornersp;              /* nodes per element */
  int subdomainsp;           /* number of partition sets */ 
  int processorsp;           /* not used */

  int *nodesp;               /* number of nodes on each PE */
  int *minep;                /* number of nodes that each PE owns */

  /* 
   * global arrays to store the permutation for reordering the nodes, and its 
   * inverse. Permutation[i] is the node number of the ith unique node to
   * appear in the p-pack file, and inv_permutation[i] is the new position of
   * node number i.
   */
  int *permutation, *inv_permutation;

} gi, *gip=&gi;


/* end globals */

/* 
 * function prototypes 
 */
/* initialization and exit routines */
void parsecommandline(int argc, char **argv, struct gi *gip);
void bail(struct gi *gip);
void finalize(struct gi *gip);

/* routines that read and parse the packfile */
void validate(struct gi *gip);
void readpack1file(struct gi *gip); 
void readpackpfile(struct gi *gip); 

/* misc routines */
void usage(struct gi *gip);
int comp_rowcol(const void *a, const void *b);


/*
 * main program 
 */
void main(int argc, char **argv) {
  char *sp;    

  /* 
   * no need to see the entire executable path name, the base will do. 
   */ 
  for (sp = argv[0]+strlen(argv[0]); (sp != argv[0]) && *sp != '/'; sp--)
    ;
  if (*sp == '/')
    strcpy(gip->progname, sp+1);
  else    
    strcpy(gip->progname, sp);

  parsecommandline(argc, argv, gip);

  /* 
   * make sure the 1-packfile and the p-packfile are OK
   */
  if (!(gip->pack1fp = fopen(gip->pack1filename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", gip->progname, gip->pack1filename);
    bail(gip);
  }
  if (!(gip->packpfp = fopen(gip->packpfilename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", gip->progname, gip->packpfilename);
    bail(gip);
  }
  validate(gip);

  fclose(gip->pack1fp);
  fclose(gip->packpfp);

  /*
   * Process the p-packfile
   */
  if (!(gip->packpfp = fopen(gip->packpfilename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", gip->progname, gip->packpfilename);
    bail(gip);
  }
  if (!gip->quiet) 
    fprintf(stderr, "%s: Reading %s.\n", argv[0], gip->packpfilename);
  readpackpfile(gip);

  /*
   * Process the 1-packfile
   */
  if (!(gip->pack1fp = fopen(gip->pack1filename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", gip->progname, gip->pack1filename);
    bail(gip);
  }
  if (!gip->quiet) 
    fprintf(stderr, "%s: Reading %s.\n", argv[0], gip->pack1filename);
  readpack1file(gip);

  fclose(gip->pack1fp);
  fclose(gip->packpfp);

  finalize(gip);
}

/*
 * validate - make sure the two input packfiles are OK
 */
void validate(struct gi *gip) {

  fscanf(gip->pack1fp, "%d %d\n", &gip->globalnodes1, &gip->mesh_dim1);
  fscanf(gip->pack1fp, "%d %d\n", &gip->globalelems1, &gip->corners1);
  fscanf(gip->pack1fp, "%d %d\n", &gip->subdomains1, &gip->processors1);

  fscanf(gip->packpfp, "%d %d\n", &gip->globalnodesp, &gip->mesh_dimp);
  fscanf(gip->packpfp, "%d %d\n", &gip->globalelemsp, &gip->cornersp);
  fscanf(gip->packpfp, "%d %d\n", &gip->subdomainsp, &gip->processorsp);
  
  /* first, make sure that both input files are indeed packfiles */
  if ((gip->subdomains1 < 1) || (gip->subdomains1 > 1024) ||
      ((gip->mesh_dim1 != 2) && (gip->mesh_dim1 != 3)) ||
      ((gip->corners1 != 3) && (gip->corners1 != 4)) ||
      gip->processors1 != gip->subdomains1) {
    fprintf(stderr, "%s: the input file %s doesn't appear to be a packfile\n",
	    gip->progname, gip->pack1filename);
    bail(gip);
  }

  if ((gip->subdomainsp < 1) || (gip->subdomainsp > 1024) ||
      ((gip->mesh_dimp != 2) && (gip->mesh_dimp != 3)) ||
      ((gip->cornersp != 3) && (gip->cornersp != 4)) ||
      gip->processorsp != gip->subdomainsp) {
    fprintf(stderr, "%s: the input file %s doesn't appear to be a packfile\n",
	    gip->progname, gip->packpfilename);
    bail(gip);
  }

  /* the first input packfile should have exactly 1 subdomain */
  if (gip->subdomains1 != 1) {
    fprintf(stderr, "%s: the first argument should be a packfile with exactly 1 subdomain\n", gip->progname);
    bail(gip);
  }

  /* the second input packfile should have more than 1 subdomain */
  if (gip->subdomainsp <= 1) {
    fprintf(stderr, "%s: the second argument should be a packfile with more than 1 subdomain\n", gip->progname);
    bail(gip);
  }

  /* make sure the 1-packfile and the p-packfile describe the same mesh */
  if ((gip->globalnodes1 != gip->globalnodesp) ||
      (gip->globalelems1 != gip->globalelemsp) ||
      (gip->corners1 != gip->cornersp) ||
      (gip->mesh_dim1 != gip->mesh_dimp)) {
    fprintf(stderr, "%s: The input packfiles don't seem to be consistent\n",
	    gip->progname);
    bail(gip);
  }
}

/*
 * readpackpfile - get the node ordering and partition from the p-packfile 
 */
void readpackpfile(struct gi *gip) {
  int i, j, k;
  int p;
  int globalnode;


  if (!(gip->permutation = (int*)malloc(sizeof(int)*gip->globalnodesp))) {
    fprintf(stderr, "%s: couldn't allocate permutation\n", 
	    gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->globalnodesp; i++) 
    gip->permutation[i] = EMPTY;
  
  if (!(gip->minep = (int *) malloc(gip->subdomainsp * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate minep\n", gip->progname);
    bail(gip);
  }

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

  /* read the global mesh info */
  fscanf(gip->packpfp, "%*d %*d %*d %*d %*d %*d\n");

  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading node partition", gip->progname);
    fflush(stderr);
  }

  /* get the number of local nodes (shared and unshared) that each PE owns */
  for (i = 0; i < gip->subdomainsp; i++) {
    fscanf(gip->packpfp, "%d %d %*d\n", &gip->nodesp[i], &gip->minep[i]);
  }

  /* 
   * set up the node permutation mapping 
   */
  p = 0;
  for (i = 0; i < gip->subdomainsp; i++) {
    if (!gip->quiet) {
      fprintf(stderr, ".");
      fflush(stderr);
    }
    for (j = 0; j < gip->nodesp[i]; j++) {
      fscanf(gip->packpfp, "%d", &globalnode);
      globalnode--;
      if ((globalnode < 0)||(globalnode > gip->globalnodesp)) {
	fprintf(stderr, "Invalid node number %d\n", globalnode+1);
	bail(gip);
      }
      for (k = 0; k < gip->mesh_dimp; k++) 
	fscanf(gip->packpfp, "%*lf");      

      if (j < gip->minep[i]) {
	if (p >= gip->globalnodesp) {
	  fprintf(stderr, "%s: packfile problem: too many nodes\n", 
		  gip->progname);
	  bail(gip);
	}
	gip->permutation[p++] = globalnode;
      }
    }
  }

  if (p != gip->globalnodesp) {
    fprintf(stderr, "%s: packfile problem: not enough nodes\n", 
	    gip->progname);
    bail(gip);
  }

  /* 
   * set up the inverse permutation mapping 
   */
  if (!(gip->inv_permutation = (int*)malloc(sizeof(int)*gip->globalnodes1))) {
    fprintf(stderr, "%s: couldn't allocate inv_permutation\n", 
	    gip->progname);
    bail(gip);
  }
  for (j = 0; j < gip->globalnodes1; j++) {
    gip->inv_permutation[j] = EMPTY;
  }
  for (j = 0; j < gip->globalnodes1; j++) {
    if (gip->inv_permutation[gip->permutation[j]] != EMPTY) {
      fprintf(stderr, "%s: packfile problem: duplicate nodes\n",
	      gip->progname);
      bail(gip);
    }
    else 
      gip->inv_permutation[gip->permutation[j]] = j;
  }
  for (j = 0; j < gip->globalnodes1; j++) {
    if (gip->inv_permutation[j] == EMPTY) {
      fprintf(stderr, "%s: packfile problem: missing node\n",
	      gip->progname);
    }
  }

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

/*
 * readpack1file - read the 1-packfile and permute the nodes according 
 *                 to the node partition specified in the p-packfile
 */
void readpack1file(struct gi *gip) {
  int oldrow, newrow;
  int i, j, loop1;
  int temp1, temp2;
  int node;
  int ival;
  rowcol_t *tuples;
  int idx;

  /* read the global mesh info */
  fscanf(gip->pack1fp, "%*d %*d %*d %*d %*d %*d\n");
  
  /* output globlal mesh info and node count */
  printf("%d %d\n", gip->globalnodes1, gip->mesh_dim1);
  printf("%d %d\n", gip->globalelems1, gip->corners1);
  printf("1\n1\n%d %d %d\n", gip->globalnodes1, gip->globalnodes1, 
	 gip->globalnodes1);

  /* read nodes */
  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading and writing nodes.", gip->progname);
  }

  fscanf(gip->pack1fp, "%*d %*d %*d");

  if (!(gip->coord1 = (void *) 
	malloc(gip->globalnodes1 * gip->mesh_dim1 * sizeof(double)))) {
    fprintf(stderr, "%s: couldn't allocate coord1(%d)\n", 
	    gip->progname, gip->globalnodes1);
    exit(0);
  }

  if (!(gip->globalnode1 = (int *) malloc(gip->globalnodes1 * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate globalnode1(%d)\n", 
	    gip->progname, gip->globalnodes1);
    exit(0);
  }
  
  for (i=0; i<gip->globalnodes1; i++) {
    fscanf(gip->pack1fp, "%d", &gip->globalnode1[i]);
    if (gip->globalnode1[i] > gip->globalnodes1 || gip->globalnode1[i] < 1) {
      fprintf(stderr, "%s: bad node number(%d). Sure this is a packfile?\n",
	      gip->progname, gip->globalnode1[i]);
	exit(0);
    }
    for (j=0; j<gip->mesh_dim1; j++) {
      fscanf(gip->pack1fp, "%lf", &gip->coord1[i][j]);
    }
  }

  for (i=0; i<gip->globalnodes1; i++) {
    idx = gip->permutation[i];
    printf("%d  ", idx+1);
    for (j=0; j<gip->mesh_dim1; j++) {
      printf("%g ", gip->coord1[idx][j]);
    }
    printf("\n");
  }

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

  /* read elements */
  if (!gip->quiet)
    fprintf(stderr, "%s: Reading and writing elements.", gip->progname);

  fscanf(gip->pack1fp, "%d", &gip->elems1);
  printf("%d\n", gip->elems1);

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

  for (i=0; i<gip->elems1; i++) {
    fscanf(gip->pack1fp, "%d", &ival);	
    printf("%d  ", ival);
    for (j=0; j<gip->corners1; j++) {
      fscanf(gip->pack1fp, "%d", &node);
      gip->vertex1[i][j] = node;
      printf("%d ", gip->inv_permutation[node]); 
    }
    printf("\n");
  }
  
  if (!gip->quiet) {
    fprintf(stderr, " Done.\n");
    fflush(stderr);
  }

  /* read sparse matrix structure */
  if (!gip->quiet)
    fprintf(stderr, "%s: Reading and writing sparse matrix structure.", gip->progname);
  fscanf(gip->pack1fp, "%d %d", &gip->matrixlen1, &temp1);
  printf("%d   %d\n", gip->matrixlen1, temp1);

  tuples = (rowcol_t*) malloc(gip->matrixlen1*sizeof(rowcol_t));
  
  oldrow = -1; 
  for (loop1 = 0; loop1 < gip->matrixlen1; loop1++) {
    fscanf(gip->pack1fp, "%d", &newrow);
    fscanf(gip->pack1fp, "%d", &node);
    tuples[loop1].row = gip->inv_permutation[newrow];
    tuples[loop1].col = gip->inv_permutation[node];

    /* convert any lower triang. entries to upper triangular entries */
    if (tuples[loop1].row > tuples[loop1].col) {
      node = tuples[loop1].col;
      tuples[loop1].col = tuples[loop1].row;
      tuples[loop1].row = node;
    }

    while (oldrow < newrow) { 
      if (oldrow+1 >= gip->globalnodes1+1) {
	printf("%s: error: (1)idx buffer too small (%d >= %d)\n", 
	       gip->progname, oldrow+1, gip->globalnodes1+1);
	exit(0);
      }
      ++oldrow;
    }
  }
  
  qsort((void*)tuples, gip->matrixlen1, sizeof(rowcol_t), comp_rowcol);
  for (loop1 = 0; loop1 < gip->matrixlen1; loop1++) {
    printf("%d %d\n", tuples[loop1].row, tuples[loop1].col);
  }
  

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

  /* read null communication info (not relevant here) */
  fscanf(gip->pack1fp, "%d %d", &temp1, &temp2);
  printf("0\n0\n0\n");
}


/* 
 * comp_rowcol - comparison function used by qsort to 
 *               compare (row, col) tuples 
 */
int comp_rowcol(const void *a, const void *b) {
  rowcol_t *a1, *b1;
  a1 = (rowcol_t*)a;
  b1 = (rowcol_t*)b;
  
  if (a1->row < b1->row) return -1;
  if (a1->row > b1->row) return 1;
  if (a1->col < b1->col) return -1;
  return 1;
  /* need not check for equal */
}

/*
 * parsecommandline - read and interpret command line arguments
 */
void parsecommandline(int argc, char **argv, struct gi *gip) {
  char *sp;    
  int i, j;

  /* no need to see the entire executable path name, the base will do. */
  for (sp = argv[0]+strlen(argv[0]); (sp != argv[0]) && *sp != '/'; sp--)
    ;
  if (*sp == '/')
    strcpy(gip->progname, sp+1);
  else    
    strcpy(gip->progname, sp);
  

  /* must have an input packfile as an argument */
  if (argc < 3) {
    usage(gip);
    bail(gip);
  }

  /* first set up the defaults */
  gip->quiet = 0;

  /* now see if the user wants to change any of these */
  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;
	}
	else if (argv[i][j] == 'h') {
	  usage(gip);
	  exit(0);
	}
	else {
	  usage(gip);
	  exit(0);
	}
      }
    }
    else {
      if ((argv[i] == NULL) || (argv[i+1] == NULL)) {
	usage(gip);
	exit(0);
      }  
      strcpy(gip->pack1filename, &argv[i][0]);
      strcpy(gip->packpfilename, &argv[i+1][0]);
      break;
    }
  }
}

void usage(struct gi *gip) {
  fprintf(stderr, "\n");
  fprintf(stderr, "Usage: %s [-Q] <old-1-packfile> <p-packfile> > <new-1-packfile>\n", gip->progname);
  fprintf(stderr, "\n");
  fprintf(stderr, "Command line args:\n");
  fprintf(stderr, "  -Q                run quietly\n");
  fprintf(stderr, "  <old-1-packfile>  original 1-subdomain packfile\n");
  fprintf(stderr, "  <p-packfile>      partitioned p-subdomain packfile\n");
  fprintf(stderr, "  <new-1-packfile>  renumbered 1-subdomain packfile\n");
  exit(0);
}


/* orderly exit */
void finalize(struct gi *gip) {
  if (!gip->quiet) {
    fprintf(stderr, "%s: Terminating normally.\n", gip->progname);
    fflush(stderr);
  }

  exit(0);
}

/* emergency exit */
void bail(struct gi *gip) {
  fprintf(stderr, "\n");
  fprintf(stderr, "%s: Something bad happened. Terminating abnormally.\n", 
	  gip->progname);
  fprintf(stderr, "\n");
  fflush(stderr);
  fflush(stdout);
  exit(0);
}