mv.c

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

  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");
#if defined(LOCK)
  printf("%s is a parallel shared memory program based on locks. The threads update\n", 
	 gip->progname);
  printf("a single output vector, using one or more lock to synchronize the updates.\n");
#elif defined(REDUCE)
  printf("%s is a parallel shared memory program based on reductions.\n",
	 gip->progname);
  printf("Each thread updates its own privatized output vector, which are\n");
  printf("later merged into a single output vector. There are no locks,\n");
  printf("and barriers are the only form of synchronization between threads.\n");
  printf("The style is similar to the Fx DO&MERGE model.\n");
#else
  printf("%s is the baseline sequential kernel.\n", gip->progname);
#endif
  printf("\n");

#if defined(LOCK)
  printf("%s [-hOQ] [-i<n>] [-l<n>] [-t<n>] packfilename\n\n", gip->progname);
#elif defined(REDUCE)
  printf("%s [-hOQ] [-i<n>] [-t<n>] packfilename\n\n", gip->progname);
#else
  printf("%s [-hOQ] [-i<n>] packfilename\n\n", gip->progname);
#endif
  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);
#if defined(LOCK)
  printf("    -l<n> Use n locks to sync the updates (default 1).\n"); 
#endif
  printf("    -O    Print the output vector to stdout.\n");
  printf("    -Q    Quietly suppress all explanations.\n");
#if (defined(LOCK) || defined(REDUCE))
  printf("    -t<n> Use n threads (default 1).\n");
#endif
  printf("\n");
  printf("Input packfiles are produced using the Archimedes tool chain. Packfiles\n");
  printf("for this program must consist of exactly 1 subdomain.\n");
  printf("\n");
#if defined(LOCK)
  printf("Example: %s -O -i10 -l256 -t8 sf5.1.pack\n", gip->progname);
#elif defined(REDUCE)
  printf("Example: %s -O -i10 -t8 sf5.1.pack\n", gip->progname);
#else
  printf("Example: %s -O -i10 sf5.1.pack\n", gip->progname);
#endif
}

void usage(struct gi *gip) {
  fprintf(stderr, "\n");
#if defined(LOCK)
  fprintf(stderr, "usage: %s [-hOQ] [-i<n>] [-l<n>] [-t<n>] packfilename\n\n", gip->progname);
#elif defined(REDUCE)
  fprintf(stderr, "usage: %s [-hOQ] [-i<n>] [-t<n>] packfilename\n\n", gip->progname);
#else
  fprintf(stderr, "usage: %s [-hOQ] [-i<n>] packfilename\n\n", gip->progname);
#endif
}


/*
 * readpackfile - reads and parses Archimedes packfile from disk
 */
void readpackfile(struct gi *gip) {
  int oldrow, newrow;
  int i, j, loop1;
  int temp1, temp2;
  int ival;
  FILE *packfile;
  char packfilename[STRLEN];

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

  /* read the basic size constants */
  fscanf(packfile, "%d", &gip->globalnodes);
  fscanf(packfile, "%d", &gip->mesh_dim);
  fscanf(packfile, "%d", &gip->globalelems);
  fscanf(packfile, "%d", &gip->corners);
  fscanf(packfile, "%d", &gip->subdomains);
  fscanf(packfile, "%d", &gip->processors);
  
  /* check for a valid packfile */
  if ((gip->subdomains < 1) || (gip->subdomains > 1024) ||
      ((gip->mesh_dim != 2) && (gip->mesh_dim != 3)) ||
      ((gip->corners != 3) && (gip->corners != 4)) ||
      gip->processors != gip->subdomains) {
    fprintf(stderr, "%s: the input file doesn't appear to be a packfile\n",
	    gip->progname);
    exit(0);
  }
  if (gip->subdomains != 1) {
    fprintf(stderr, "%s: You are using a packfile with %d subdomains. It should have exactly 1.\n",
	    gip->progname, gip->subdomains);
    exit(0);
  }
  
  /* read nodes */
  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading nodes.", gip->progname);
  }

  fscanf(packfile, "%d %d %d", &gip->nodes, &gip->mine, &gip->priv);

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

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

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

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

  fscanf(packfile, "%d", &gip->elems);

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

  for (i=0; i<gip->elems; i++) {
    fscanf(packfile, "%d", &ival);	
    for (j=0; j<gip->corners; j++) {
      fscanf(packfile, "%d", &gip->vertex[i][j]);
    }
  }

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

  /* read sparse matrix structure and convert from tuples to CSR */
  if (!gip->quiet)
    fprintf(stderr, "%s: Reading sparse matrix structure.", gip->progname);
  fscanf(packfile, "%d %d", &gip->matrixlen, &gip->choleskylen);
  
  /* allocate sparse matrix structures */
  if (!(gip->matrixindex = (int *) malloc((gip->nodes+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate matrixindex(%d)\n", 
	    gip->progname, gip->nodes+1);
    exit(0);
  }
  if (!(gip->matrixcol = (int *) malloc((gip->matrixlen+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate matrixcol(%d)\n", 
	    gip->progname, gip->matrixlen + 1);
    exit(0);
  }
  oldrow = -1; 
  for (loop1 = 0; loop1 < gip->matrixlen; loop1++) {
    fscanf(packfile, "%d", &newrow);
    fscanf(packfile, "%d", &gip->matrixcol[loop1]);
    while (oldrow < newrow) { 
      if (oldrow+1 >= gip->globalnodes+1) {
	printf("%s: error: (1)idx buffer too small (%d >= %d)\n", 
	       gip->progname, oldrow+1, gip->globalnodes+1);
	exit(0);
      }
      gip->matrixindex[++oldrow] = loop1;
    }
  }
  while (oldrow < gip->globalnodes) {
    gip->matrixindex[++oldrow] = gip->matrixlen;
  }
  
  if (!gip->quiet) {
    fprintf(stderr, " Done.\n");
    fflush(stderr);
  }

  /* read null communication info (not relevant here) */
  fscanf(packfile, "%d %d", &temp1, &temp2);
}

/*
 * utility routines 
 */
void printnodevector(double (*v)[3], int n) {
  int i;
  for (i=0; i<n; i++)
    printf("%d %.0f\n", i+1, v[i][0]); 
  fflush(stdout);
}


/*
 * 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 thread (pthread) interface
 */
#if defined(PTHREAD)

void spark_init_threads() {
  int i, 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_cond_init(&barriercondhandle, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: couldn't initialize barrier cond\n", gip->progname);
    exit(0);
  }

  if (!(vectormutexhandle = (pthread_mutex_t *) 
	malloc(gip->locks * sizeof(pthread_mutex_t)))) {
    fprintf(stderr, "%s: couldn't allocate vectormutexhandle\n", gip->progname);
    exit(0);
  }
  for (i=0; i<gip->locks; i++) {
    status = pthread_mutex_init(&vectormutexhandle[i], NULL);
    if (status != 0) {
      fprintf(stderr, "%s: couldn't initialize vector mutex %d\n", gip->progname, i);
      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->threads) {
    barriercnt = 0;
    pthread_cond_broadcast(&barriercondhandle);
  }
  else {
    pthread_cond_wait(&barriercondhandle, &barriermutexhandle);
  }
  pthread_mutex_unlock(&barriermutexhandle);
}

void spark_setlock(int lockid) {
  pthread_mutex_lock(&vectormutexhandle[lockid]);
}

void spark_unsetlock(int lockid) {
  pthread_mutex_unlock(&vectormutexhandle[lockid]);
}
#endif

/*
 * SGI Threads Routines 
 */
#if defined(SGI)
void spark_init_threads() {
  int i;

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

  if (!(lockhandle = (ulock_t *) malloc(gip->locks * sizeof(ulock_t)))) {
    fprintf(stderr, "%s: couldn't allocate lockhandle\n", gip->progname);
    exit(0);
  }
  for (i=0; i<gip->locks; i++) {
    lockhandle[i] = usnewlock(arenahandle);
    if (lockhandle[i] == NULL) {
      fprintf(stderr, "Couldn't create lock %d\n", i);
      exit(0);
    }
    usinitlock(lockhandle[i]);
  }
}

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->threads);
}

void spark_setlock(int lockid) {
  ussetlock(lockhandle[lockid]);
  
}

void spark_unsetlock(int lockid) {
  usunsetlock(lockhandle[lockid]);
}
#endif