hmv.c

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

 *    v and w are vectors of elements of size n x DOF.
 *    smvp is computed starting at row firstrow (zero based)
 *    and for numrows rows.
 * 
 *    K is a block symmetric sparse matrix which is stored
 *    in compressed sparse row format as three vectors:
 *      - A is a coefficient vector of elements of size DOFxDOF. Only 
 *        the upper right triangle is stored.
 *      - Acol[i] is the column number of the coefficient in K[i].
 *      - Row i consists of elements A[Acol[i]]...A[Acol[i+1]-1].
 */
void local_smvp(int nodes, double (*A)[DOF][DOF], int *Acol, 
		int *Aindex, double (*v)[DOF], double (*w)[DOF],
		int firstrow, int numrows) {
  int i;
  int Anext, Alast, col;
  double sum0, sum1, sum2;
  
  for (i = firstrow; i < (firstrow + numrows); i++) {
    Anext = Aindex[i];
    Alast = Aindex[i + 1];

    sum0 = A[Anext][0][0]*v[i][0] + A[Anext][0][1]*v[i][1] + A[Anext][0][2]*v[i][2];
    sum1 = A[Anext][1][0]*v[i][0] + A[Anext][1][1]*v[i][1] + A[Anext][1][2]*v[i][2];
    sum2 = A[Anext][2][0]*v[i][0] + A[Anext][2][1]*v[i][1] + A[Anext][2][2]*v[i][2];

    Anext++;
    while (Anext < Alast) {
      col = Acol[Anext];

      sum0 += A[Anext][0][0]*v[col][0] + A[Anext][0][1]*v[col][1] + A[Anext][0][2]*v[col][2];
      sum1 += A[Anext][1][0]*v[col][0] + A[Anext][1][1]*v[col][1] + A[Anext][1][2]*v[col][2];
      sum2 += A[Anext][2][0]*v[col][0] + A[Anext][2][1]*v[col][1] + A[Anext][2][2]*v[col][2];
      
      w[col][0] += A[Anext][0][0]*v[i][0] + A[Anext][1][0]*v[i][1] + A[Anext][2][0]*v[i][2];
      w[col][1] += A[Anext][0][1]*v[i][0] + A[Anext][1][1]*v[i][1] + A[Anext][2][1]*v[i][2];
      w[col][2] += A[Anext][0][2]*v[i][0] + A[Anext][1][2]*v[i][1] + A[Anext][2][2]*v[i][2];

      Anext++;
    }

    w[i][0] += sum0;
    w[i][1] += sum1;
    w[i][2] += sum2;
  }
}

/*
 * full_assemble - assemble a distributed vector
 */
void full_assemble(double (*v)[DOF], int id, struct gi *gip) {
  int i, j, s, pos;
  int overwrite;

  spark_barrier();
  
  /* copy partial results to our neighbors' receive buffers */
  for (s = 0; s < gip->subdomains; s++) {
    if (gip->commindex[id][s] < gip->commindex[id][s+1]) {
      pos = gip->commindex[s][id];
      for (j = gip->commindex[id][s]; j < gip->commindex[id][s+1]; j++) {
	gip->recvbuf[s][pos][0] = v[gip->comm[id][j]][0];
	gip->recvbuf[s][pos][1] = v[gip->comm[id][j]][1];
	gip->recvbuf[s][pos++][2] = v[gip->comm[id][j]][2];
      }
    }
  }

  spark_barrier();
  
  /* update the copies of my node variables */
  for (i = 0; i < gip->commlen[id]; i++) {
    v[gip->comm[id][i]][0] += gip->recvbuf[id][i][0];
    v[gip->comm[id][i]][1] += gip->recvbuf[id][i][1];
    v[gip->comm[id][i]][2] += gip->recvbuf[id][i][2];
  }

  spark_barrier();
}

void init(int argc, char **argv, struct gi *gip) {
  int i;
  char *sp;    

  init_etime();
  spark_init_threads();

  /* 
   * no need to see the entire 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);

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

  load_pack(gip);

  if (!(gip->recvbuf = (void *)malloc(gip->subdomains * sizeof(double *)))) {
    fprintf(stderr, "%s: couldn't allocate recvbuf\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->recvbuf[i] = (void *)malloc((gip->maxcommlen+1)*DOF*sizeof(double)))) {
      fprintf(stderr, "%s: couldn't allocate recvbuf[i]\n", gip->progname);
      bail(gip);
    }
  }

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

/*
 * load_pack - load the input packfile
 */
void load_pack(struct gi *gip) {
  if (!(gip->packfile = fopen(gip->packfilename, "r"))) {
    fprintf(stderr, "%s: Can't open %s\n", 
	    gip->progname, gip->packfilename);
    bail(gip);
  }
  
  load_global(gip);
  load_nodes(gip);
  load_elems(gip);
  load_matrix(gip);
  load_comm(gip);
}

void load_global(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);

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

#ifdef DEBUGGLOBAL    
  seq_prglobal(gip);
#endif
}

void load_nodes(struct gi *gip) {
  int i, j, k;

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

  for (i = 0; i < gip->subdomains; i++) {
    fscanf(gip->packfile, "%d %d %d\n", &gip->nodes[i], &gip->mine[i], &gip->priv[i]);
  }

  if (!(gip->coord = (void *) malloc(gip->subdomains*sizeof(double *)))) {
    fprintf(stderr, "%s: couldn't allocate coord\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->coord[i] = (void *) malloc(gip->nodes[i]*gip->mesh_dim*sizeof(double)))) {
      fprintf(stderr, "%s: couldn't allocate coord[i]\n", gip->progname);
      bail(gip);
    }
  }

  if (!(gip->globalnode = (void *) malloc(gip->subdomains * sizeof(int *)))) {
    fprintf(stderr, "%s: couldn't alloc globalnode\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->globalnode[i] = (int *) malloc(gip->nodes[i]*sizeof(int)))) {
      fprintf(stderr, "%s: couldn't alloc globalnode[i]\n", gip->progname);
      bail(gip);
    }
  }

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

    /* node data */
    for (j = 0; j < gip->nodes[i]; j++) {
      fscanf(gip->packfile, "%d", &gip->globalnode[i][j]);
      for (k = 0; k < gip->mesh_dim; k++) {
	fscanf(gip->packfile, "%lf", &gip->coord[i][j][k]);
      }
    }
  }

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

#ifdef DEBUGNODES
  seq_prnodes(gip);
#endif
}

void load_elems(struct gi *gip) {
  int i, j, k;

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

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

  if (!(gip->globalelem = (void *) malloc(gip->subdomains * sizeof(int *)))) {
    fprintf(stderr, "%s: couldn't alloc globalelem\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->globalelem[i] = (void *) malloc(gip->elems[i] * sizeof(int)))) {
      fprintf(stderr, "%s: couldn't alloc globalelem\n", gip->progname);
      bail(gip);
    }
  }

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

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

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

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

#ifdef DEBUGELEMS
  seq_prelems(gip);
#endif
}  

void load_matrix(struct gi *gip) {
  int i, k, loop1;
  int oldrow, newrow;

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

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

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

  if (!(gip->matrixcol = (void *) malloc(gip->subdomains * sizeof(int *)))) {
    fprintf(stderr, "%s: couldn't allocate matrixcol\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->matrixcol[i] = (void *) malloc((gip->matrixlen[i]+1) * sizeof(int)))) {
      fprintf(stderr, "%s: couldn't allocate matrixcol[i]\n", gip->progname);
      bail(gip);
    }
  }

  if (!(gip->matrixindex = (void *) malloc(gip->subdomains * sizeof(int *)))){
    fprintf(stderr, "%s: couldn't allocate matrixindex\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->matrixindex[i] = (int *) malloc((gip->nodes[i]+1) * sizeof(int)))){
      fprintf(stderr, "%s: couldn't allocate matrixindex\n", gip->progname);
      bail(gip);
    }
  }
    
  for (i = 0; i < gip->subdomains; i++) {

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

    oldrow = -1; 
    for (loop1 = 0; loop1 < gip->matrixlen[i]; loop1++) {

      fscanf(gip->packfile, "%d", &newrow);
      fscanf(gip->packfile, "%d", &gip->matrixcol[i][loop1]);
      while (oldrow < newrow) {
	if (oldrow+1 >= gip->matrixlen[i]) {
	  fprintf(stderr, "%s: index buffer(1) too small (%d >= %d)\n", 
		  gip->progname, oldrow+1, gip->matrixlen[i]);
	  bail(gip);
	}
	gip->matrixindex[i][++oldrow] = loop1;
      }
    }
    while (oldrow < gip->nodes[i]) {
      if (oldrow+1 >= gip->matrixlen[i]) {
	fprintf(stderr, "%s: index buffer(2) too small (%d >= %d)\n", 
		gip->progname, oldrow+1, gip->matrixlen[i]);
	bail(gip);
      }
      gip->matrixindex[i][++oldrow] = gip->matrixlen[i];
    }

  }

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

#ifdef DEBUGMATRIX
  seq_prmatrix(gip);
#endif
}

void load_comm(struct gi *gip) {
  int i, j, k;
  int count;
  int *commnodes; /*subdomain i shares a total of commnodes[i] nodes */
  int *buf;
  
  if (!gip->quiet) {
    fprintf(stderr, "%s: Reading communication info", gip->progname);
    fflush(stderr);
  }
  
  if (!(gip->friends = (int *) malloc(gip->subdomains * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't allocate friends\n", gip->progname);
    bail(gip);
  }

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

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

  gip->maxcommlen = 0;
  for (i = 0; i < gip->subdomains; i++) {
    fscanf(gip->packfile, "%d\n", &commnodes[i]);
    if (commnodes[i] > gip->maxcommlen)
      gip->maxcommlen = commnodes[i];
  }
  
  if (!(gip->comm = (void *) malloc(gip->subdomains * sizeof(int *)))) {
    fprintf(stderr, "%s: couldn't allocate comm\n", gip->progname);
    bail(gip);
  }
  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->comm[i] = (int *) malloc((commnodes[i]+1) * sizeof(int)))) {
      fprintf(stderr, "%s: couldn't allocate comm[i]\n", gip->progname);
      bail(gip);
    }
  }

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

  for (i=0; i<gip->subdomains; i++) {
    if (!(gip->commindex[i] = (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[i] = 0;
    gip->friends[i] = 0;

    for (j = 0; j < gip->subdomains; j++) {
      
      gip->commindex[i][j] = gip->commlen[i];

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

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