hmv.c

卡内基梅隆大学开发的并行计算程序

    if (gip->commlen[i] != commnodes[i]) {
      fprintf(stderr, "%s: inconsistent comm lengths\n", gip->progname);
      bail(gip);
    }
  }

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

#ifdef DEBUGCOMM
    seq_prcomm(gip);
#endif
}    

/*
 * parsecommandline - read and interpret command line arguments
 */
void parsecommandline(int argc, char **argv, struct gi *gip) {
  int i, j;
    
  /* must have a file name */
  if (argc < 2) {
    usage(gip);
  }

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

  /* now see if the user wants to change any of these */
  for (i=1; i<argc; i++) {
    if (argv[i][0] == '-') {
      if (argv[i][1] == 'i') {
	gip->iters = atoi(&argv[i][2]);
	if (gip->iters <= 0) {
	  fprintf(stderr, "error: iterations must be greater than zero.\n");
          fprintf(stderr, "no spaces allowed after the -i (e.g. -i100).\n");
	  bail(gip);
	}
      }
      else {
	for (j = 1; argv[i][j] != '\0'; j++) {
	  if (argv[i][j] == 'Q') {
	    gip->quiet = 1;
	  }
	  else if ((argv[i][j] == 'h' ||argv[i][j] == 'H')) {
	    info(gip);
	    finalize(gip);
	  }
	  else if (argv[i][j] == 'O') {
	    gip->output = 1;
	  }
	  else {
	    usage(gip);
	    bail(gip);
	  }
	}
      }
    }
    else {
      strcpy(gip->packfilename, &argv[i][0]);
    }
  }
}

/*
 * info and usage - explain the command line arguments
 */
void info(struct gi *gip) {
  printf("\n");
  printf("You are running the %s kernel from the Spark98 Benchmarks.\n", 
	 gip->progname);
  printf("Copyright (C) 1998, David O'Hallaron, Carnegie Mellon University.\n");
  printf("You are free to use this software without restriction. If you find that\n");
  printf("the suite is helpful to you, it would be very helpful if you sent me a\n");
  printf("note at [e-mail] letting me know how you are using it.\n"); 
  printf("\n");
  printf("%s is a hybrid program based on shared memory but written in a\n", gip->progname);
  printf("message passing style. Each thread performs its own local SMVP\n");
  printf("operation using partially assembled local matrices and vectors\n");
  printf("followed by an assembly phase that combines the partially assembled\n");
  printf("output vectors.\n");
  printf("\n");
  printf("%s [-hOQ] [-i<n>] packfilename\n\n", gip->progname);
  printf("Command line options:\n\n");
  printf("    -h    Print this message and exit.\n");
  printf("    -i<n> Do n iterations of the SMVP pairs (default %d).\n",
	 DEFAULT_ITERS);
  printf("    -O    Print the output vector to stdout.\n");
  printf("    -Q    Quietly suppress all explanations.\n");
  printf("\n");
  printf("Input packfiles are produced using the Archimedes tool chain.\n");
  printf("\n");
  printf("Example: %s -O -i10 sf5.8.pack\n", gip->progname);
}

void usage(struct gi *gip) {
  fprintf(stderr, "\n");
  fprintf(stderr, "usage: %s [-OQh] [-i<n>] packfilename\n\n", gip->progname);
  exit(0);
}

/* 
 * sequential output routines
 */

/* 
 * printnodevector - print one column of a global nodevector
 */
void printnodevector(double (**v)[DOF], int n, struct gi *gip) {
  int i, j, k, found;

  for (k=1; k<=n; k++) {
    found = 0;
    for (i=0; i<gip->subdomains && !found; i++) {
      for (j=0; j<gip->nodes[i]; j++) {
	if (k == gip->globalnode[i][j]) {
	  printf("%d %.0f\n", k, v[i][j][0]);
	  found = 1;
	  break;
	}
      }
    }
  }
}

/* 
 * printvec3 - print the local version of a vector on each subdomain
 */
void printvec3(double (**v)[DOF], int id, struct gi *gip) {
  int i, s;
  
  for (s=0; s<gip->subdomains; s++) {
    for (i = 0; i < gip->nodes[id]; i++) {
      printf("[%d]: %d (%d): %.0f\n", s, i, 
	     gip->globalnode[id][i], v[id][i][0]);
    }
    fflush(stdout);
  }
}


/*
 * parallel output routines
 */

/* emit local matrix in order on each subdomain */
void par_printmatrix3(double (*A)[DOF][DOF], int id, struct gi *gip) {
  int i, j, k, l, row, col, s;

  for (s=0; s<gip->subdomains; s++) {
    spark_barrier();
    if (id == s) {
      for (i = 0; i < gip->nodes[s]; i++) { 
	for (j = 0; j < 3; j++) {
	  for (k = gip->matrixindex[s][i]; k < gip->matrixindex[s][i+1]; k++) {
	    for (l = 0; l < 3; l++) {
	      row =  i*3 + j + 1; /*row*/
	      col =  gip->matrixcol[s][k]*3 + l + 1; /*col*/
	      printf("%d: %d %d %.0f\n", s, row, col, A[k][j][l]);
	      printf("%d: %d %d %.0f\n", s, col, row, A[k][j][l]);
	    }
	  }
	}
      }
    }
  }
}

/* emit local vector in order on each subdomain */
void par_printvec3(double (**v)[DOF], int n, int id, struct gi *gip) {
  int i, s;
  
  for (s=0; s<gip->subdomains; s++) {
    spark_barrier();
    if (id == s) {
      printf("%d: ready to print\n", id);
      for (i = 0; i < n; i++) {
	printf("[%d]: %d (%d): %.0f\n", id, i, 
	       gip->globalnode[s][i], v[s][i][0]);
      }
      fflush(stdout);
    }
  }
}

/*
 * sequential debug routines for the pack file input
 */

/* print global info */
void seq_prglobal(struct gi *gip) {
  int i;

  printf("gnodes=%d dim=%d gelem=%d corners=%d subdomains=%d procs=%d\n",
	 gip->globalnodes, gip->mesh_dim, 
	 gip->globalelems, gip->corners, gip->subdomains, 
	 gip->processors);
  fflush(stdout);
}

/* print finite element mesh nodes on each subdomain */
void seq_prnodes(struct gi *gip) {
  int i, j, k; 
  
  for (i = 0; i < gip->subdomains; i++) {
    printf("[%d]: gip->nodes=%d gip->priv=%d gip->mine=%d\n", 
	   i, gip->nodes[i], gip->priv[i], gip->mine[i]);
    for (j = 0; j < gip->nodes[i]; j++) {
      printf("[%d] %d (%d): ", i, j, gip->globalnode[i][j]);
      for (k = 0; k < gip->mesh_dim; k++) {
	printf("%f ", gip->coord[i][j][k]);
      }
      printf("\n");	
    }
    fflush(stdout);
  }
}

/* print finite element mesh elements on each subdomain */
void seq_prelems(struct gi *gip)
{
  int i, j, k;
    
  for (i = 0; i < gip->subdomains; i++) {
    printf("[%d]: gip->elems=%d\n", i, gip->elems[i]);  
    for (j = 0; j < gip->elems[i]; j++) {
      printf("[%d]: %d (%d): ", i, j, gip->globalelem[i][j]);
      for (k = 0; k < gip->corners; k++) {
	printf("%d ", gip->vertex[i][j][k]);
      }
      printf("\n");
      fflush(stdout);
    }
  }
}

/* print sparse matrix  structure on each subdomain */
void seq_prmatrix(struct gi *gip)
{
  int i, j;
    
  for (i = 0; i < gip->subdomains; i++) {
    fflush(stdout);
    printf("[%d]: gip->matrixcol:\n", i);	
    for (j = 0; j < gip->matrixlen[i]; j++)
      printf("[%d]: %d: %d\n", i, j, gip->matrixcol[i][j]); 
    printf("[%d]: gip->matrixindex:\n", i);	
    for (j = 0; j <= gip->nodes[i]; j++)
      printf("[%d]: %d: %d\n", i, j, gip->matrixindex[i][j]);
    fflush(stdout);
  }
}

/* print communication schedule on each subdomain */
void seq_prcomm(struct gi *gip) {
  int i, j, k;
  
  for (i = 0; i < gip->subdomains; i++) {
    printf("[%d]: comm info (%d shared nodes, %d friends)\n", 
	   i, gip->commlen[i], gip->friends[i]);
    for (j = 0; j < gip->subdomains; j++) {
      printf("[%d]: %d nodes shared with subdomain %d:\n", 
	     i, gip->commindex[i][j+1]-gip->commindex[i][j],j);
      fflush(stdout);
      for (k = gip->commindex[i][j]; k < gip->commindex[i][j+1]; k++) {
	printf("[%d]: %d %d (%d)\n", 
	       i, k, gip->comm[i][k], gip->globalnode[i][gip->comm[i][k]]);
	fflush(stdout);
      }
    }
    fflush(stdout);
  }
}

/* print communication schedule with values of the receive buffer */
void seq_prcommvals(struct gi *gip) {
  int i, j, k;
  
  for (i = 0; i < gip->subdomains; i++) {
    printf("[%d]: comm info (%d shared nodes, %d friends)\n", 
	   i, gip->commlen[i], gip->friends[i]);
    for (j = 0; j < gip->subdomains; j++) {
      printf("[%d]: %d nodes shared with subdomain %d:\n", 
	     i, gip->commindex[i][j+1]-gip->commindex[i][j],j);
      fflush(stdout);
      for (k = gip->commindex[i][j]; k < gip->commindex[i][j+1]; k++) {
	printf("[%d]: %d %d (%d) recvbuf[%d][%d]=%.0f\n", 
	       i, k, gip->comm[i][k], gip->globalnode[i][gip->comm[i][k]],
	       i, k, gip->recvbuf[i][k][0]);
	fflush(stdout);
      }
    }
    fflush(stdout);
  }
}

/*
 * termination routines 
 */

/* 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);
}


/*
 * system dependent timer routines
 */
#define MAX_ETIME 86400   /* 24 hour timer period */
struct itimerval first_u; /* user time */

/* init the timer */
void init_etime(void) {
  first_u.it_interval.tv_sec = 0;
  first_u.it_interval.tv_usec = 0;
  first_u.it_value.tv_sec = MAX_ETIME;
  first_u.it_value.tv_usec = 0;
  setitimer(ITIMER_VIRTUAL, &first_u, NULL);
}

/* return elapsed user seconds since call to init_etime */
double get_etime(void) {
  struct itimerval curr;
  
  getitimer(ITIMER_VIRTUAL, &curr);
  return (double) ((first_u.it_value.tv_sec - curr.it_value.tv_sec) +
		   (first_u.it_value.tv_usec - curr.it_value.tv_usec)*1e-6);
}


/*
 * System dependent thread routines
 */

/* 
 * Posix Threads Routines (pthread)
 */
#if defined(PTHREAD)
void spark_init_threads() {
  int status; 

  status = pthread_mutex_init(&barriermutexhandle, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: couldn't initialize barrier mutex\n", gip->progname);
    exit(0);
  }
  status = pthread_mutex_init(&vectormutexhandle, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: couldn't initialize vector mutex\n", gip->progname);
    exit(0);
  }
  status = pthread_cond_init(&barriercondhandle, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: couldn't initialize barrier cond\n", gip->progname);
    exit(0);
  }
}

void spark_start_threads(int n) {
  int i, status;
  pthread_t threadhandle;
  
  for (i=1; i<=n; i++) {
    ids[i] = i;
    status = pthread_create(&threadhandle, NULL, smvpthread, (void *)&ids[i]);
    if (status != 0) {
      fprintf(stderr, "%s: Could not create thread %d\n", gip->progname, i);
      exit(0);
    }
  }
}

void spark_barrier(void) {
  pthread_mutex_lock(&barriermutexhandle);
  barriercnt++;
  if (barriercnt == gip->subdomains) {
    barriercnt = 0;
    pthread_cond_broadcast(&barriercondhandle);
  }
  else {
    pthread_cond_wait(&barriercondhandle, &barriermutexhandle);
  }
  pthread_mutex_unlock(&barriermutexhandle);
}

void spark_setlock() {
  pthread_mutex_lock(&vectormutexhandle);
}

void spark_unsetlock() {
  pthread_mutex_unlock(&vectormutexhandle);
}
#endif

/*
 * SGI Threads Routines 
 */
#if defined(SGI)
void spark_init_threads() {
  arenahandle = usinit(arenaname);
  if (arenahandle == NULL) {
    fprintf(stderr, "arena init failed\n");
    exit(0);
  }
  barrierhandle = new_barrier(arenahandle);
  if (barrierhandle == NULL) {
    fprintf(stderr, "barrier init failed\n");
    exit(0);
  }
  lockhandle = usnewlock(arenahandle);
  if (lockhandle == NULL) {
    fprintf(stderr, "lock init failed\n");
    exit(0);
  }
  usinitlock(lockhandle);
}

void spark_start_threads(int n) {
  int i;
  
  for (i=1; i<=n; i++) {
    ids[i] = i;
    (void)sproc((void (*)(void*)) smvpthread, PR_SADDR, (void *)&ids[i]);
  }
}

void spark_barrier(void) {
  barrier(barrierhandle, gip->subdomains);
}

void spark_setlock() {
  ussetlock(lockhandle);
}

void spark_unsetlock() {
  usunsetlock(lockhandle);
}
#endif