mmv.c

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

      bail(gip);
    }

    /* now make sure we have a valid packfile and that there are enough */
    /* runtime processes for the problem */
    fscanf(gip->packfile, "%d %d\n", &gip->globalnodes, &gip->mesh_dim);
    fscanf(gip->packfile, "%d %d\n", &gip->globalelems, &gip->corners);
    fscanf(gip->packfile, "%d\n", &gip->subdomains);
    fscanf(gip->packfile, "%d\n", &gip->processors);

    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);
      bail(gip);
    }

    if ((gip->subdomains) > gip->rsize) {
      fprintf(stderr, "%s: The mesh consists of %d subdomains. Please\n",
	      gip->progname, gip->subdomains);
      fprintf(stderr, "rerun the program using at least %d processes.\n", 
	      gip->subdomains);
      bail(gip);
    }
  }

  /* let all processes know how the number of subdomains in the mesh */
  MPI_Bcast(&gip->subdomains, 1, MPI_INT, 0, MPI_COMM_WORLD);

  /* 
   * Create the compute subgroup (one process per subdomain) 
   */
  crange[0][0] = 0;
  crange[0][1] = gip->subdomains-1;
  crange[0][2] = 1;
  MPI_Group_range_incl(gip->rgroup, 1, crange, &(gip->cgroup)); 
  MPI_Group_rank(gip->cgroup, &crank);
  MPI_Group_size(gip->cgroup, &(gip->csize));
  MPI_Comm_create(MPI_COMM_WORLD, gip->cgroup, &(gip->ccomm));
  gip->ioproc = gip->csize-1;
  
  /* 
   * Run as a compute process 
   */
  if (gip->rank < gip->subdomains) {
    if ((gip->rank != crank)) {
      fprintf(stderr, "%s: inconsistent ranks on compute process\n", 
	      gip->progname);
      bail(gip);
    }	    

    /* input the packfile */
    load_pack(gip);

    /* 
     * preallocate some MPI status and request buffers  that are
     * needed for the asynchronous communication phase 
     */
    if (!(gip->statusvec = (MPI_Status *) 
	  malloc(gip->subdomains * sizeof(MPI_Status)))) {
      fprintf(stderr, "%s: couldn't allocate statusvec\n", gip->progname);
      bail(gip);
    }
    
    if (!(gip->requestvec = (MPI_Request *) 
	  malloc(gip->subdomains * sizeof(MPI_Request)))) {
      fprintf(stderr, "%s: couldn't allocate request\n", gip->progname);
      bail(gip);
    }
    
    /*
     * preallocate send and receive buffers for the communication phase 
     */
    if (gip->maxcommnodes > 0) {
      if (!(gip->sendbuf = (double *)
	    malloc(gip->maxcommnodes * 3 * sizeof(double)))) {
	fprintf(stderr, "%s: couldn't allocate sendbuf\n", gip->progname);
	bail(gip);
      }
      
      if (!(gip->recvbuf = (double *) 
	    malloc(gip->maxcommnodes * 3 * sizeof(double)))) {
	fprintf(stderr, "%s: couldn't allocate recvbuf\n", gip->progname);
	bail(gip);
      }
    }
  }
    
  /* 
   * Run as a dummy process, wait for compute processes, 
   * then terminate. 
   */
  else {
    finalize(gip);
  }
}

/*
 * load_pack - load and distribute the packfile
 */
void load_pack(struct gi *gip) {
  struct stat st;
  
  if (stat(gip->packfilename, &st) < 0) {
    fprintf(stderr, "%s: Can't find %s \n", 
	    gip->progname, gip->packfilename);
    bail(gip);
  }
  
  if (!(gip->packfile = fopen(gip->packfilename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", 
	    gip->progname, gip->packfilename);
    bail(gip);
  }
  
  if (gip->rank == gip->ioproc) {
    load_global_master(gip);
    load_nodes_master(gip);
    load_elems_master(gip);
    load_matrix_master(gip);
    load_comm_master(gip);
  }
  else {
    load_global_slave(gip);
    load_nodes_slave(gip);
    load_elems_slave(gip);
    load_matrix_slave(gip);
    load_comm_slave(gip);
  }
}

void load_global_master(struct gi *gip) {
  int ibuf[6];
  
  fscanf(gip->packfile, "%d %d\n", &gip->globalnodes, &gip->mesh_dim);
  fscanf(gip->packfile, "%d %d\n", &gip->globalelems, &gip->corners);
  fscanf(gip->packfile, "%d\n", &gip->subdomains);
  fscanf(gip->packfile, "%d\n", &gip->processors);
  ibuf[0] = gip->globalnodes;
  ibuf[1] = gip->mesh_dim;
  ibuf[2] = gip->globalelems;
  ibuf[3] = gip->corners;
  ibuf[4] = gip->subdomains;
  ibuf[5] = gip->processors;

  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);
    bail(gip);
  }

  MPI_Bcast(ibuf, 6, MPI_INT, gip->ioproc, gip->ccomm);
  
#ifdef DEBUGGLOBAL    
  prglobal(gip);
#endif
}

void load_global_slave(struct gi *gip) {
  int ibuf[6];
  
  MPI_Bcast(ibuf, 6, MPI_INT, gip->ioproc, gip->ccomm);
  gip->globalnodes = ibuf[0];
  gip->mesh_dim = ibuf[1];
  gip->globalelems = ibuf[2];
  gip->corners = ibuf[3];
  gip->subdomains = ibuf[4];
  gip->processors = ibuf[5];
  
#ifdef DEBUGGLOBAL    
  prglobal(gip);
#endif
}

void load_nodes_master(struct gi *gip) {
  int i, j, k;
  int *nodepriv;
  int *nodemine;
  double *fbuf, *fbufp;
  int *buf, *bufp;
  int wordint;
  double wordfloat;

  MPI_Barrier(gip->ccomm);

  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading nodes", gip->progname);
    fflush(stderr);
  }
  if (!(gip->nodeall = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate nodeall\n", 
	    gip->progname);
    bail(gip);
  }
  if (!(nodepriv = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate nodepriv\n", gip->progname);
    bail(gip);
  }
  if (!(nodemine = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate nodemine\n", gip->progname);
    bail(gip);
  }

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

  MPI_Bcast(&gip->maxnodes, 1, MPI_INT, gip->ioproc, gip->ccomm);
  MPI_Bcast(gip->nodeall, gip->subdomains, MPI_INT, gip->ioproc, gip->ccomm);
    
  if (!(fbuf = (double *)malloc(gip->maxnodes*gip->mesh_dim*sizeof(double)))) {
    fprintf(stderr, "%s: couldn't allocate local node buffer\n", 
	    gip->progname);
    bail(gip);
  }

  if (!(buf = (int *) malloc(gip->maxnodes * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't alloc global node index buffer\n", 
	    gip->progname);
    bail(gip);
  }
    
  for (i = 0; i < gip->subdomains; i++) {
    if (!gip->quiet) {
      fprintf(stderr, ".");
      fflush(stderr);
    }

    /* node header */
    bufp = buf;
    *bufp++ = gip->nodeall[i];
    *bufp++ = nodemine[i];
    *bufp++ = nodepriv[i];
    MPI_Send(buf, (int)(bufp-buf), MPI_INT, i, 0, gip->ccomm);

    /* node data */
    bufp = buf;
    fbufp = fbuf;
    for (j = 0; j < gip->nodeall[i]; j++) {

      /* global node number */
      fscanf(gip->packfile, "%d", &wordint);
      *bufp++ = wordint;

      /* node data */
      for (k = 0; k < gip->mesh_dim; k++) {
	fscanf(gip->packfile, "%lf", &wordfloat);
	if ((fbufp-fbuf) > gip->maxnodes * gip->mesh_dim) {
	  fprintf(stderr, "%s: node buffers too small (%d > %d)\n", 
		  gip->progname, (int)(fbufp-fbuf),
		  gip->maxnodes * gip->mesh_dim);
	  bail(gip);
	}
	*fbufp++ = wordfloat;
      }
    }
    if (i < gip->ioproc) {
      MPI_Send(buf, (int)(bufp-buf), MPI_INT, i, 0, gip->ccomm);
      MPI_Send(fbuf, (int)(fbufp-fbuf), MPI_DOUBLE, i, 0, gip->ccomm);
    }
    else {
      gip->nodes = gip->nodeall[i];
      gip->priv = nodepriv[i];
      gip->mine = nodemine[i];
      gip->globalnode = (int *)buf;
      gip->coord = (void *)fbuf;
    }
  }

  (void)free(nodepriv);
  (void)free(nodemine);

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

#ifdef DEBUGNODES
  prnodes(gip);
#endif
}

void load_nodes_slave(struct gi *gip) {
  int ibuf[3];
  MPI_Status status;
  
  MPI_Barrier(gip->ccomm);

  MPI_Bcast(&gip->maxnodes, 1, MPI_INT, gip->ioproc, gip->ccomm);

  if (!(gip->nodeall = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate nodeall\n", 
	    gip->progname);
    bail(gip);
  }
  MPI_Bcast(gip->nodeall, gip->subdomains, MPI_INT, gip->ioproc, gip->ccomm);
  MPI_Recv(ibuf, 3, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);
  
  gip->nodes = ibuf[0];
  gip->mine = ibuf[1]; 
  gip->priv = ibuf[2]; 
  
  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);
    bail(gip);
  }
  
  if (!(gip->globalnode = (int *) malloc(gip->nodes * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't alloc globalnode\n", 
	    gip->progname);
    bail(gip);
  }
  
  MPI_Recv(gip->globalnode, gip->nodes, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);

  MPI_Recv(gip->coord, gip->nodes * gip->mesh_dim, MPI_DOUBLE, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);
  

#ifdef DEBUGNODES
  prnodes(gip);
#endif 
}

void load_elems_master(struct gi *gip) {
  int i, j, k;
  int wordint;
  int *elemcnt;
  int maxelems;
  int *buf, *bufp;
  int *buf2, *buf2p;

  MPI_Barrier(gip->ccomm);

  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading elements", gip->progname);
    fflush(stderr);
  }
  if (!(elemcnt = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate elemcnt\n", gip->progname);
    bail(gip);
  }

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

  if (!(buf2 = (int *) malloc(maxelems * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't alloc global element buf\n", 
	    gip->progname);
    bail(gip);
  }

  for (i = 0; i < gip->subdomains; i++) {

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

    /* element header */
    bufp = buf;
    *bufp++ = elemcnt[i];
    MPI_Send(buf, 1, MPI_INT, i, 0, gip->ccomm);
	
    /* element data */
    bufp = buf;
    buf2p = buf2;
    for (j = 0; j < elemcnt[i]; j++) {

      /* global elem number */
      fscanf(gip->packfile, "%d", &wordint);
      *buf2p++ = wordint;
    
      /* elem data */
      for (k = 0; k < gip->corners; k++) {
	fscanf(gip->packfile, "%d", &wordint);
	if ((bufp-buf) > maxelems*gip->corners) {
	  fprintf(stderr, "%s: elem buffers too small\n", 
		  gip->progname);
	  bail(gip);
	}
	*bufp++ = wordint;
      }
    }

    if (i < gip->ioproc) {
      MPI_Send(buf2, (int)(buf2p-buf2), MPI_INT, i, 0, gip->ccomm);
      MPI_Send(buf, (int)(bufp-buf), MPI_INT, i, 0, gip->ccomm);
    }
    else {
      gip->elems = elemcnt[i];
      gip->globalelem = (int *)buf2;
      gip->vertex = (void *)buf;
    }
  }

  (void)free(elemcnt);

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

#ifdef DEBUGELEMS
  prelems(gip);
#endif
}  

void load_elems_slave(struct gi *gip) {
  MPI_Status status;
  
  MPI_Barrier(gip->ccomm);

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

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

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

  MPI_Recv(gip->globalelem, gip->elems, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status);
  
  MPI_Recv(gip->vertex, gip->elems*gip->corners, MPI_INT, 
	   gip->ioproc, MPI_ANY_TAG, gip->ccomm, &status); 
  
#ifdef DEBUGELEMS
  prelems(gip);
#endif
}

void load_matrix_master(struct gi *gip) {
  int i, k, loop1;
  int oldrow, newrow;
  int maxcsrwords;
  int *csrwords;
  int *buf;
  int *idxbuf;

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

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

  /* get the number of i,j pairs on each subdomain */
  maxcsrwords = 0;
  for (i = 0; i < gip->subdomains; i++) { 
    fscanf(gip->packfile, "%d %*d\n", &csrwords[i]); 
    if (csrwords[i] > maxcsrwords)
      maxcsrwords = csrwords[i];
  }

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

  if (!(idxbuf = (int *) malloc((gip->maxnodes+1) * sizeof(int)))){
    fprintf(stderr, "%s: couldn't allocate index buffer\n", gip->progname);
    bail(gip);
  }
    
  /* build the CSR format for each subdomain and send it */
  MPI_Barrier(gip->ccomm);
  for (i = 0; i < gip->subdomains; i++) {

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

    MPI_Send(&csrwords[i], 1, MPI_INT, i, 0, gip->ccomm);
    oldrow = -1; 
    for (loop1 = 0; loop1 < csrwords[i]; loop1++) {

      fscanf(gip->packfile, "%d", &newrow);
      fscanf(gip->packfile, "%d", &buf[loop1]);
      while (oldrow < newrow) {
	if (oldrow+1 >= maxcsrwords) {
	  fprintf(stderr, 
		  "%s: index buffer(1) too small (%d >= %d)\n", 
		  gip->progname, oldrow+1, maxcsrwords);
	  bail(gip);
	}
	idxbuf[++oldrow] = loop1;
      }
    }
    while (oldrow < gip->nodeall[i]) {
      if (oldrow+1 >= maxcsrwords) {
	fprintf(stderr, "%s: index buffer(2) too small (%d >= %d)\n", 
		gip->progname, oldrow+1, maxcsrwords);
	bail(gip);
      }
      idxbuf[++oldrow] = csrwords[i];
    }
    if (i < gip->ioproc) {
      MPI_Send(buf, csrwords[i], MPI_INT, i, 0, gip->ccomm);