mmv.c

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

      MPI_Send(idxbuf, gip->nodeall[i]+1, MPI_INT, i, 0, gip->ccomm);
    } 
    else {
      gip->matrixlen = csrwords[i];
      gip->matrixcol = (void *)buf;
      gip->matrixindex = (void *)idxbuf;
    }
  }   

  (void)free(csrwords);

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

#ifdef DEBUGMATRIX
  prmatrix(gip);
#endif
}

void load_matrix_slave(struct gi *gip) {
  MPI_Status status;

  MPI_Barrier(gip->ccomm);

  MPI_Recv(&gip->matrixlen, 1, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);

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

  if (!(gip->matrixindex = (int *) malloc((gip->nodes+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate matrixindex\n", gip->progname);
    bail(gip);
  }
    
  MPI_Recv(gip->matrixcol, gip->matrixlen, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);

  MPI_Recv(gip->matrixindex, gip->nodes+1, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);


#ifdef DEBUGMATRIX
  prmatrix(gip);
#endif
}

void load_comm_master(struct gi *gip) {
  int i, j, k;
  int count;
  int *commnodes; /*subdomain i shares a total of commnodes[i] nodes */
  int *buf;
  
  MPI_Barrier(gip->ccomm);

  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading communication info", gip->progname);
    fflush(stderr);
  }
  
  if (!(commnodes = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate commnodes\n", gip->progname);
    bail(gip);
  }
  
  gip->maxcommnodes = 0;
  for (i = 0; i < gip->subdomains; i++) {
    fscanf(gip->packfile, "%d\n", &commnodes[i]);
    if (commnodes[i] > gip->maxcommnodes) {
      gip->maxcommnodes = commnodes[i];
    }
  }
  
  MPI_Bcast(&gip->maxcommnodes, 1, MPI_INT, gip->ioproc, gip->ccomm);

  if (!(gip->comm = (int *) malloc((gip->maxcommnodes+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate comm\n", gip->progname);
    bail(gip);
  }
  
  if (!(gip->commindex = (int *) malloc((gip->subdomains+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate commindex\n", gip->progname);
    bail(gip);
  }
  

  for (i = 0; i < gip->subdomains; i++) {
    
    if (!gip->quiet) {
      fprintf(stderr, ".");
      fflush(stderr);
    }
    
    gip->commlen = 0;
    gip->friends = 0;
    for (j = 0; j < gip->subdomains; j++) {
      
      gip->commindex[j] = gip->commlen;

      /* subdomain i shares count nodes with subdomain j */
      fscanf(gip->packfile, "%d", &count);
      if (count > gip->maxcommnodes) {
	fprintf(stderr, "%s: count exceeds maxcommnodes\n", gip->progname);
	bail(gip);
      }
      
      if (count > 0)
	gip->friends++;

      for (k = 0; k < count; k++) {
	fscanf(gip->packfile, "%d", &gip->comm[gip->commlen++]);
      }
    }
    gip->commindex[gip->subdomains] = gip->commlen;

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

    if (i < gip->ioproc) {
	MPI_Send(&gip->commlen, 1, MPI_INT, i, 0, gip->ccomm);
	MPI_Send(&gip->friends, 1, MPI_INT, i, 0, gip->ccomm);
	MPI_Send(gip->comm, gip->commlen, MPI_INT, i, 0, gip->ccomm);
	MPI_Send(gip->commindex, gip->subdomains+1, MPI_INT, i, 0, gip->ccomm);
    }
  }

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

#ifdef DEBUGCOMM
    prcomm(gip);
#endif
}    

void load_comm_slave(struct gi *gip) {
  MPI_Status status;
  
  MPI_Barrier(gip->ccomm);
  MPI_Bcast(&gip->maxcommnodes, 1, MPI_INT, gip->ioproc, gip->ccomm);
  MPI_Recv(&gip->commlen, 1, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);
  MPI_Recv(&gip->friends, 1, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);

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

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

  MPI_Recv(gip->comm, gip->commlen, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);
  MPI_Recv(gip->commindex, gip->subdomains+1, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);

#ifdef DEBUGCOMM
    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);
    bail(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) {
  if (gip->rank == 0) {
    printf("\n");
    printf("You are running the %s kernel from the Spark98 Kernels.\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");
    printf("the suite is helpful to you, it would be very helpful if you sent me a ");
    printf("note at [e-mail] letting me know how you are using it.          "); 
    printf("\n");
    printf("%s is a parallel message passing program. Each processor performs\n",gip->progname);
    printf("its own local SMVP operation using partially assembled local matrices\n");
    printf("and vectors, followed by an assembly phase that combines the partially\n");
    printf("assembled 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: mpirun -np 8 %s -O -i10 sf5.8.pack\n", gip->progname);
  }
}

void usage(struct gi *gip) {
  if (gip->rank == 0) {
    printf("\n");
    printf("usage: %s [-hOQ] [-i<n>] packfilename\n\n", gip->progname);
  }
}


/*
 * printnodevector - print one column of a global nodevector
 */
void printnodevector(double (*v)[DOF], int n, struct gi *gip) {
  MPI_Status status;
  int i, j, nitems;
  int cmd = SEND_CMD;
  double *gvec, *fvec;
  int *ivec;


  MPI_Barrier(gip->ccomm);

  /* allocate scratch buffers for print operations */
  if (!(fvec = (double *)malloc((gip->maxnodes) * sizeof(double)))) {
    fprintf(stderr, "%s: couldn't allocate fvec\n", gip->progname);
    bail(gip);
  }


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

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

  /* i/o master */
  if (gip->rank == 0) {
    for (j=1; j<=gip->globalnodes; j++) {
      gvec[j] = INVALID;
    }

    for (j=0; j<gip->mine; j++) {
      gvec[gip->globalnode[j]] = v[j][0];
    }
    for (i=1; i<gip->subdomains; i++) { 
      MPI_Send(&cmd, 1, MPI_INT, i, 0, gip->ccomm);
      MPI_Probe(i, MPI_ANY_TAG, gip->ccomm, &status);

      MPI_Get_count(&status, MPI_DOUBLE, &nitems);
      MPI_Recv(fvec, nitems, MPI_DOUBLE, i, MPI_ANY_TAG, gip->ccomm, &status);
      MPI_Recv(ivec, nitems, MPI_INT, i, MPI_ANY_TAG, gip->ccomm, &status);
      for (j=0; j<nitems; j++) {
	if (gvec[ivec[j]] != INVALID) {
	  fprintf(stderr, 
		  "warning: received duplicate entry in gvec[%d] from %d\n", 
		  ivec[j], i);
	}
	gvec[ivec[j]] = fvec[j];
      }
    }
    for (j=1; j<=n; j++) {
      if (gvec[j] == INVALID) {
	fprintf(stderr, "error: uninitialized value in gvec[%d]\n", j);
	bail(gip);
      }
      printf("%d %.0f\n", j, gvec[j]); 
    }
  }

  /* i/o slave */
  else {
    for (j=0; j<gip->mine; j++) {
      fvec[j] = v[j][0];
    }
    MPI_Send(fvec, gip->mine, MPI_DOUBLE, 0, 0, gip->ccomm);
    MPI_Send(gip->globalnode, gip->mine, MPI_INT, 0, 0, gip->ccomm);
  }

  free(gvec);
  free(fvec);
  free(ivec);
}

/*
 * local output routines for debugging
 */

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

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

/* emit local vector on each subdomain */
void local_printvec3(double (*v)[DOF], int n, struct gi *gip) {
  int i, j;
  
  for (i=0; i < gip->subdomains; i++) {
    MPI_Barrier(gip->ccomm);
    if (i == gip->rank) {
      for (j = 0; j < n; j++) {
	printf("[%d]: %d (%d): %.0f\n", 
	       gip->rank, j, gip->globalnode[j], v[j][0]);
      }
      fflush(stdout);
    }
  }
}

/*
 * debug routines that print out various local data structures one
 * subdomain at a time during the initialization phase.
 */
void prglobal(struct gi *gip)
{
  int i;

  for (i = 0; i < gip->csize; i++) {
    MPI_Barrier(gip->ccomm);
    if (i == gip->rank) {
      printf("[%d]: gnodes=%d dim=%d gelem=%d corners=%d subdomains=%d procs=%d\n",
	     gip->rank, gip->globalnodes, gip->mesh_dim, 
	     gip->globalelems, gip->corners, gip->subdomains, 
	     gip->processors);
    }
    fflush(stdout);
  }
}

void prnodes(struct gi *gip) {
  int i, j, k; 
  
  for (i = 0; i < gip->subdomains; i++) {
    MPI_Barrier(gip->ccomm);
    if (i == gip->rank) {
      printf("[%d]: gip->nodes=%d gip->priv=%d gip->mine=%d\n", 
	     gip->rank, gip->nodes, gip->priv, gip->mine);
      for (j = 0; j < gip->nodes; j++) {
	printf("[%d] %d (%d): ", gip->rank, j, gip->globalnode[j]);
	for (k = 0; k < gip->mesh_dim; k++) {
	  printf("%f ", gip->coord[j][k]);
	}
	printf("\n");	
      }
      fflush(stdout);
    }
  }
}

void prelems(struct gi *gip)
{
  int i, j, k;
    
  for (i = 0; i < gip->csize; i++) {
    MPI_Barrier(gip->ccomm);
    if (i == gip->rank) {
      printf("[%d]: gip->elems=%d\n", gip->rank, gip->elems);  
      for (j = 0; j < gip->elems; j++) {
	printf("\n");    
	printf("[%d]: %d (%d): ", gip->rank, j, gip->globalelem[j]);
	for (k = 0; k < gip->corners; k++) {
	  printf("%d ", gip->vertex[j][k]);
	}
      }
      printf("\n");
      fflush(stdout);
    }
  }
}

void prmatrix(struct gi *gip)
{
  int i, j;
    
  for (i = 0; i < gip->csize; i++) {
    MPI_Barrier(gip->ccomm);
    if (gip->rank == i) {
      fflush(stdout);
      printf("[%d]: gip->matrixcol:\n", gip->rank);	
      for (j = 0; j < gip->matrixlen; j++)
	printf("[%d]: %d: %d\n", gip->rank, j, gip->matrixcol[j]); 
      printf("[%d]: gip->matrixindex:\n", gip->rank);	
      for (j = 0; j <= gip->nodes; j++)
	printf("[%d]: %d: %d\n", gip->rank, j, gip->matrixindex[j]);
      fflush(stdout);
    }
  }
}

void prcomm(struct gi *gip) {
  int i, j, k;
  
  for (i = 0; i < gip->subdomains; i++) {
    MPI_Barrier(gip->ccomm);
    if (i == gip->rank) {
      printf("[%d]: comm info (%d shared nodes, %d friends)\n", 
	     gip->rank, gip->commlen, gip->friends);
      for (j = 0; j < gip->subdomains; j++) {
	printf("[%d]: %d nodes shared with subdomain %d:\n", 
	       gip->rank, gip->commindex[j+1]-gip->commindex[j],j);
	fflush(stdout);
	for (k = gip->commindex[j]; k < gip->commindex[j+1]; k++) {
	  printf("[%d]: %d (%d)\n", gip->rank, 
		 gip->comm[k], gip->globalnode[gip->comm[k]]);
	  fflush(stdout);
	}
      }
      fflush(stdout);
    }
  }
}


/*
 * termination routines 
 */

/* orderly exit */
void finalize(struct gi *gip) 
{
  MPI_Barrier(MPI_COMM_WORLD);
  if ((gip->rank == 0) && !gip->quiet) {
    fprintf(stderr, "%s: Terminating normally.\n", gip->progname);
    fflush(stderr);
  }
  MPI_Finalize();
  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);
  MPI_Abort(MPI_COMM_WORLD, 0);
}