mv.c

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

  exit(0);
}

/*
 * init - initialize the kernel
 */
void init(int argc, char **argv, struct gi *gip) {

  init_etime();

  gip->progname = argv[0];
  parsecommandline(argc, argv, gip);

  if (!gip->quiet) 
    fprintf(stderr, "%s: Reading %s.\n", argv[0], gip->packfilename);
  readpackfile(gip);

  if (gip->locks == 0)
    gip->locks = gip->nodes;

  if (!(firstrow = (int *)malloc((gip->threads+1) * sizeof(int)))) {
    fprintf(stderr, "%s: couldn't malloc firstrow\n", gip->progname); 
    exit(0);
  }

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

#if (defined(LOCK) || defined(REDUCE))
  spark_init_threads();
#endif

}

/*
 * smvpthread - perform a sequence of SMVP pairs 
 */
void *smvpthread(void *a) {
  int i;
  int *argp = (int *)a;
  int id = *argp;
  double mycsecs, mystarttime, mystartctime; 

#if (defined(LOCK) || defined(REDUCE))
  int numrows = firstrow[id+1] - firstrow[id];
  spark_barrier();
#endif

  mycsecs = 0.0;
  mystarttime = get_etime();

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

#if defined(LOCK) /* parallel with one or more locks */

    /* w1 = K1*v1 */
    spark_barrier();
    zero_vector(w1, firstrow[id], numrows);
    spark_barrier();
    local_smvp(gip->nodes, K1, gip->matrixcol, gip->matrixindex, 
	       v1, w1, firstrow[id], numrows, id, gip);

    /* w2 = K2*v2 */
    spark_barrier();
    zero_vector(w2, firstrow[id], numrows);
    spark_barrier();
    local_smvp(gip->nodes, K2, gip->matrixcol, gip->matrixindex, 
	       v2, w2, firstrow[id], numrows, id, gip);

#elif defined(REDUCE) /* parallel with reductions instead of locks */

    /* w1 = K1*v1 */
    spark_barrier();
    zero_vector(tmpw[id], 0, gip->nodes);
    local_smvp(gip->nodes, K1, gip->matrixcol, gip->matrixindex, 
	       v1, tmpw[id], firstrow[id], numrows, id, gip);

    mystartctime = get_etime();
    spark_barrier();
    add_vectors(tmpw, w1, firstrow[id], numrows, gip); 
    mycsecs += (get_etime() - mystartctime);

    /* w2 = K2*v2 */
    spark_barrier();
    zero_vector(tmpw[id], 0, gip->nodes);
    local_smvp(gip->nodes, K2, gip->matrixcol, gip->matrixindex, 
	       v2, tmpw[id], firstrow[id], numrows, id, gip);

    mystartctime = get_etime();
    spark_barrier();
    add_vectors(tmpw, w2, firstrow[id], numrows, gip); 
    mycsecs += (get_etime() - mystartctime);

#else /* sequential */

    /* w1 = K1*v1 */
    zero_vector(w1, 0, gip->nodes);
    local_smvp(gip->nodes, K1, gip->matrixcol, gip->matrixindex, 
	       v1, w1, 0, gip->nodes, id, gip);

    /* w2 = K2*v2 */
    zero_vector(w2, 0, gip->nodes);
    local_smvp(gip->nodes, K2, gip->matrixcol, gip->matrixindex, 
	       v2, w2, 0, gip->nodes, id, gip);
#endif

  }

  if (id == 0) {
    secs = (get_etime() - mystarttime)/(2.0*gip->iters);
    csecs = mycsecs / (2.0 * gip->iters);
  }

#if (defined(LOCK) || defined(REDUCE))
  spark_barrier();
#endif

  return NULL;
}

/*
 * assemble_matrix - assemble the local sparse matrix
 */
void assemble_matrix(double (*K)[DOF][DOF], struct gi *gip) {
  int i, j, k, l, m;
  int temp1, elem;

  for (i = 0; i < gip->matrixlen; i++) {
    for (j = 0; j < DOF; j++) {
      for (k = 0; k < DOF; k++) {
	K[i][j][k] = 0.0;
      }
    }
  }

  /* 
   * add the contribution of each element to K
   */
  for (elem = 0; elem < gip->elems; elem++) {
    for (j = 0; j < gip->corners; j++) {
      for (k = 0; k < gip->corners; k++) {
	if (gip->vertex[elem][j] < gip->vertex[elem][k]) {
	  temp1 = gip->matrixindex[gip->vertex[elem][j]];
	  while (gip->matrixcol[temp1] != gip->vertex[elem][k]) {
	    temp1++;
	    if (temp1 >= gip->matrixindex[gip->vertex[elem][k] + 1]) {
	      fprintf(stderr, "%s: K indexing error in assemble %d %d\n", 
		     gip->progname, temp1, gip->vertex[elem][k]);
	      exit(0);
	    }
	  }
	  for (l = 0; l < 3; l++) {
	    for (m = 0; m < 3; m++) {
	      K[temp1][l][m]++;
	    }
	  }
	}
      }
    }

  } 

#ifdef DEBUGASSEMBLE
  local_printmatrix3(K, gip);
#endif
}

/*
 * assemble_vector - assemble the local v vector
 */
void assemble_vector(double (*v)[DOF], struct gi *gip) {
  int i, j;

  for (i = 0; i < gip->nodes; i++) {
    for (j = 0; j < DOF; j++) {
      v[i][j] = 1.0;
    }
  } 

#ifdef DEBUGASSEMBLE
  local_printvec3(v, gip->nodes, gip);
#endif
}

/*
 * zero_vector - clear a portion of a vector 
 */
void zero_vector(double (*v)[DOF], int firstrow, int numrows) {
  int i;
  for (i=firstrow; i<(firstrow+numrows); i++) {
    v[i][0] = 0.0;
    v[i][1] = 0.0;
    v[i][2] = 0.0;
  }
}

/*
 * local_smvp - local sparse matrix vector product w = K*v
 *    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 id, struct gi *gip) {
  int i;
  int Anext, Alast, col;
  double sum0, sum1, sum2;

#if defined(LOCK)
  int lockid;
#endif  

  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];
      
#if defined(LOCK)      
      lockid = col % gip->locks;
      spark_setlock(lockid);
#endif

      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];

#if defined(LOCK)      
      spark_unsetlock(lockid);
#endif

      Anext++;
    }

#if defined(LOCK)      
    lockid = i % gip->locks;
    spark_setlock(lockid);
#endif

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

#if defined(LOCK)      
    spark_unsetlock(lockid);
#endif
  }
}

/*
 * add_vectors - add the vectors produced by each process
 */
void add_vectors(double (**tmpw)[DOF], double (*w)[DOF], 
		 int firstrow, int numrows, struct gi *gip) {
  int i, j;
  double sum0, sum1, sum2;
  
  for (i = firstrow; i < (firstrow+numrows); i++) {
    sum0 = sum1 = sum2 = 0.0;
    for (j=0; j<gip->threads; j++) {
      sum0 += tmpw[j][i][0];
      sum1 += tmpw[j][i][1];
      sum2 += tmpw[j][i][2];
    }
    w[i][0] = sum0;
    w[i][1] = sum1;
    w[i][2] = sum2;
  }
}

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

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

  /* 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");
	  exit(0);
	}
      }
#if (defined(LOCK) || defined(REDUCE))
      else if (argv[i][1] == 't') {
	gip->threads = atoi(&argv[i][2]);
	if (gip->threads <= 0) {
	  fprintf(stderr, "error: must have at least 1 thread.\n");
          fprintf(stderr, "no spaces allowed after the -t (e.g. -t8).\n");
	  exit(0);
	}
      }
#endif
#if defined(LOCK)
      else if (argv[i][1] == 'l') {
	gip->locks = atoi(&argv[i][2]);
	if (gip->locks < 0) {
	  fprintf(stderr, "error: must have at least 1 lock.\n");
          fprintf(stderr, "no spaces allowed after the -l (e.g. -l8).\n");
	  exit(0);
	}
      }
#endif

      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);
	    exit(0);
	  }
	  else if (argv[i][j] == 'O') {
	    gip->output = 1;
	  }
	  else {
	    usage(gip);
	    exit(0);
	  }

	}
      }
    }
    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 Kernels.\n", gip->progname); 
  printf("Copyright (C) 1998, David O'Hallaron, Carnegie Mellon University.\n");