grid.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

	(void) printf("(%3d,%10.4f) ",j+row_low,grid[i][j]);
	if (++newline == 4) {
	  (void) printf("\n"); newline = 0;
	}
      }
    }
    if (newline) 
      (void) printf("\n");
    (void) printf("\n");
  }
}

double Operate(grid, ncols, nrows, ngrid, do_sums)
     long ncols, nrows, ngrid, do_sums;
     double **grid;
/*
  Update grid in place according to the simple rule

     new[i][j] = 0.25 * 
       (old[i-1][j] + old[i+1][j] + old[i][j-1] + old[i][j+1]);

  for 0<i<ncols, 0<j<nrows.

  To improve convergence rate for large grids use red/black
  checkerboard -> update red then update black using new reds
  (gauss-seidel red-black w-relaxation).

  Return the mean abs. error (value of the laplacian) on the old grid.
*/
{
  double residual = 0.0;
  double omega = 0.94;          /* Over relaxation parameter */
  long i,j, type=10, length=1, jlo;
  double *gg, *ggm, *ggp;

  Exchange(grid, ncols, nrows);
  
  /* Update red */
  for (i=1; i<ncols; i++) {
    jlo = 1+((i-1+col_low+row_low)%2);
    gg = grid[i]; ggm = grid[i-1]; ggp = grid[i+1];
#ifdef ALLIANT
#pragma safe gg
#pragma safe ggp
#pragma safe ggm
#endif
#ifdef ARDENT
#pragma IVDEP
#endif
    for (j=jlo; j<nrows; j+=2) {
      double new, diff;
      new = 0.25 * (ggp[j] + ggm[j] + gg[j-1] + gg[j+1]);
      new = new + omega*(new - gg[j]);
      diff = new - gg[j];
      residual += ABS(diff);
      gg[j] = new;
    }
  }

  Exchange(grid, ncols, nrows);

  /* Update black */
  for (i=1; i<ncols; i++) {
    jlo = 1+((i+col_low+row_low)%2);
    gg = grid[i]; ggm = grid[i-1]; ggp = grid[i+1];
#ifdef ALLIANT
#pragma safe gg
#pragma safe ggp
#pragma safe ggm
#endif
#ifdef ARDENT
#pragma IVDEP
#endif
    for (j=jlo; j<nrows; j+=2) {
      double new, diff;
      new = 0.25 * (ggp[j] + ggm[j] + gg[j-1] + gg[j+1]);
      new = new + omega*(new - gg[j]);
      diff = new - gg[j];
      residual += ABS(diff);
      gg[j] = new;
    }
  }

  if (do_sums) {
    long start = MTIME_();
    DGOP_(&type, &residual, &length, "+");
    cs_global += MTIME_() - start;

    return 4.0 * residual / (scale*scale*(ngrid-1)*(ngrid-1));
  }
  else
    return 999999.0;
}

double **Make2d(ncols, nrows)
     long ncols, nrows;
/*
  C does not allow variably dimensioned multi-dim. arrays.
  Allocate memory to simulate 2-D array with array of pointers.
*/
{
  double **temp, **array;

  array = temp = (double **) malloc((unsigned) (ncols*sizeof(double *)));

  if (!temp)
    Error("Make2d: malloc 1 failed",(long) (ncols*sizeof(double *)));

  while (ncols--) {
    if (!(*temp++ = (double *) malloc((unsigned) (nrows*sizeof(double)))))
      Error("Make2d: malloc 2 failed", (long) (nrows*sizeof(double)));
  }
  return array;
}

void Interpolate(old_grid, old_ncols, old_nrows, old_col_low, old_row_low,
		 new_grid, new_ncols, new_nrows, new_col_low, new_row_low)
     double **old_grid, **new_grid;
     long old_ncols, old_nrows, old_col_low, old_row_low;
     long new_ncols, new_nrows, new_col_low, new_row_low;
/*
  Interpolate down from a solution on grid ncols/2 to an initial
  guess on a grid dimension ncols.
  The parallel version is complicated by the fact that rows may
  have shifted between processes on scaling up. Ugh.
*/
{
  long col_shift = 2*old_col_low - new_col_low;
  long row_shift = 2*old_row_low - new_row_low;
  long start = MTIME_();
  long i, i1, i2, j, j1, j2;

  for (i=0; i<=old_ncols; i++)
    for (j=0; j<=old_nrows; j++) {
      i1 = 2*i+col_shift; j1 = 2*j+row_shift;
      if ( (i1>=0) && (i1<=new_ncols) && (j1>=0) && (j1<=new_nrows) )
	new_grid[i1][j1] = old_grid[i][j];
    }

  for (i=1; i<=old_ncols; i++)
    for (j=1; j<=old_nrows; j++) {
      i2 = 2*i-1+col_shift; j2 = 2*j-1+row_shift;
      if ( (i2>=0) && (i2<=new_ncols) && (j2>=0) && (j2<=new_nrows) )
	new_grid[i2][j2] = 0.25 * (old_grid[i  ][j  ] + old_grid[i-1][j-1] +
				   old_grid[i-1][j  ] + old_grid[i  ][j-1]);
    }
  
  BoundaryConditions(new_grid, new_ncols, new_nrows);
  Exchange(new_grid, new_ncols, new_nrows);

  for (i=1; i<new_ncols; i++)
    for (j=1+((i+new_col_low+new_row_low)%2); j<new_nrows; j+=2)
      new_grid[i][j] = 0.25 * (new_grid[i+1][j] + 
			       new_grid[i-1][j] + 
			       new_grid[i][j-1] + 
			       new_grid[i][j+1]);

  Exchange(new_grid, new_ncols, new_nrows);

  cs_interpolate += MTIME_() - start;
  
}

void Solve(grid, ncols, nrows, ngrid, niter, nprint, thresh)
     double **grid, thresh;
     long ncols, nrows, ngrid, niter, nprint;
/*
  Apply iterative procedure to solve current grid.

  Parallel version only does global sums every nsums iterations.
*/
{
  long iter, do_sums, nsums;
  double residual, error;
#ifdef PLOT
  long nplots = MAX(5, ngrid/10);  /* Only plot when have changed a lot */
#endif

  nsums = MIN(10, nprint);     /* Need sums whenever we print */
  if (nprint%nsums)
    nprint= nprint + nsums - (nprint%nsums); /* Make nprint a multiple 
					        of nsums */
  for (iter=0; iter<niter; iter++) {

    /* For efficiency only do global sums every 10 iters */
    do_sums = !(iter%nsums);

    /* Actually do the work */
    residual = Operate(grid, ncols, nrows, ngrid, do_sums);

    /* Print the results every now and again or if converged */
    if ((NODEID_()==0) && ((iter%nprint == 0) || (residual < thresh))) {
      (void) printf("ngrid=%d iter=%d residual=%f\n",
		    ngrid, iter+1, residual);
      (void) fflush(stdout);
    }
    
    /* Are we converged ? */
    if (do_sums && (residual < thresh)) {
      if (NODEID_() == 0) {
	(void) printf("Converged!\n");
	(void) fflush(stdout);
      }
      break;
    }

#ifdef PLOT
    if (!(iter%nplots))
      PlotGrid(grid, ngrid, ncols, nrows);
#endif

  }
  
  /* Have to do a final exchange to get edges correct */
  Exchange(grid, ncols, nrows);

  error = GridError(grid, ncols, nrows, ngrid);
  if (NODEID_() == 0) {
    (void) printf("Mean abs. error to exact soln. = %f, ngrid=%d\n\n",
		  error, ngrid);
    (void) fflush(stdout);
  }

#ifdef PLOT
  PlotGrid(grid, ngrid, ncols, nrows);
#endif
}

void ParseArguments(argc, argv, pngrid, pniter, pnprint, pnlevel, pthresh)
     long argc, *pngrid, *pniter, *pnprint, *pnlevel;
     double *pthresh;
     char **argv;
{
  argc--;
  argv++;
  while (argc--) {
    if (strcmp(*argv, "-niter") == 0)
      *pniter = atoi(*++argv);
    else if (strcmp(*argv, "-ngrid") == 0)
      *pngrid = atoi(*++argv);
    else if (strcmp(*argv, "-nprint") == 0)
      *pnprint = atoi(*++argv);
    else if (strcmp(*argv, "-nlevel") == 0)
      *pnlevel = atoi(*++argv);
    else if (strcmp(*argv, "-thresh") == 0)
      (void) sscanf(*++argv, "%lf", pthresh);
#ifdef PLOT
    else if (strcmp(*argv, "-plot") == 0) {
      argv++;
      if (strcmp(*argv,"value") == 0)
	plot_type = PLOT_VALUE;
      else if (strcmp(*argv, "error") == 0)
	plot_type = PLOT_ERROR;
      else if (strcmp(*argv, "residual") == 0)
	plot_type = PLOT_RESIDUAL;
      else
	Error("Unknown plot type - use error|value|residual",(long) -1);
    }	
#endif
    else if (strcmp(*argv, "-help") == 0) {
      (void) fprintf(stderr,"gridtest [-ngrid #] [-nprint #] [-niter #]\n");
      (void) fprintf(stderr,"         [-thresh #] [-nlevel #] [-help]\n");
#ifdef PLOT
      (void) fprintf(stderr,"         [-plot value|error|residual]\n");
#endif
      PEND_();
      exit(1);
    }
    argv++; argc--;
  }
}

void Factor(N, pn, pm)
  long N, *pn, *pm;
/*
  Factor N into two integers that are as close together as possible
*/
{
  long n = ceil(sqrt((double) N));
  long m = floor(sqrt((double) N));
 
  while (n*m != N) {
    if (n*m < N)
      n++;
    else
      m--;
  }
  *pn = n; *pm = m;
}

void Partition(ngrid, pncols, pnrows)
    long ngrid, *pncols, *pnrows;
/*
  The square mesh actually has (ngrid-1)*(ngrid-1) interior points.
  Divide these up equitably between all the processors.
*/
{
  long col_chunk = (ngrid-1) / ncol_P;
  long col_extra = (ngrid-1) - col_chunk*ncol_P;
  long row_chunk = (ngrid-1) / nrow_P;
  long row_extra = (ngrid-1) - row_chunk*nrow_P;

  col_low = my_col_P*col_chunk + MIN(my_col_P,col_extra); /* Col of top LHS */

  row_low = my_row_P*row_chunk + MIN(my_row_P,row_extra); /* Row of top LHS */

  *pncols = col_chunk + ( (my_col_P<col_extra) ? 1 : 0 ) + 1;
  *pnrows = row_chunk + ( (my_row_P<row_extra) ? 1 : 0 ) + 1;
}

int main(argc, argv)
     int argc;
     char **argv;
/*
  Grid test code. Solve Laplaces eqn. on a square grid subject
  to b.c.s on the boundary.  Use 5 point discretization of the
  operator and a heirarchy of grids with red/black gauss seidel
  w-relaxation (omega=0.94 as 1.0 diverges).

  command arguments:

  -ngrid  # = dimension of largest grid (default = 128)
  -niter  # = max. no. of interations   (default = 1000)
  -nprint # = print residual every nprint iterations (default = 100)
  -nlevel # = no. of grid levels        (default = 4)
  -thresh # = convergence criterion on the residual (default = 0.1)
  -help     = print usage and exit with error
  -plot value|error|residual = if compiled with -DPLOT generate plots
                               displaying 
  
*/
{
  double **grid, thresh=0.1,**new_grid;
  long niter=1000, nprint=100, ngrid, ncols=0, nrows=0, maxgrid=128;
  long nlevel=4, level, ngrid1, on=0;
  long old_ncols, old_nrows, old_col_low, old_row_low;
  long start = MTIME_();

  PBEGIN_(argc, argv);            /* Initialize parallel environment */
  SETDBG_(&on);

  P = NNODES_();
  Factor(P, &ncol_P, &nrow_P);    /* Arrange processes into a grid   */
  my_col_P = NODEID_()/nrow_P;
  my_row_P = NODEID_() - my_col_P*nrow_P; 

  /* Who are my neighbors on the process grid */

  north = (my_row_P > 0)          ? NODEID_()-1      : -1;
  south = (my_row_P < (nrow_P-1)) ? NODEID_()+1      : -1;
  east  = (my_col_P < (ncol_P-1)) ? NODEID_()+nrow_P : -1;
  west  = (my_col_P > 0)          ? NODEID_()-nrow_P : -1;

  ParseArguments(argc, argv, 
		 &maxgrid, &niter, &nprint, &nlevel, &thresh);

  ngrid1 = MAX(ncol_P,maxgrid>>(nlevel-1));
  ngrid1 = MAX(nrow_P,ngrid1);                 /* Size of first grid */
  maxgrid = ngrid1<<(nlevel-1);                /* Actual size of final grid */

  if (!(buffer = (double *) malloc((unsigned) (sizeof(double)*(maxgrid+1)))))
    Error("failed to allocate comms buffer", (long) (maxgrid+1));

#ifdef PLOT
  OpenPlotFile(maxgrid);
#endif

  /* Loop from coarse to fine grids */

  for (level=0; level<nlevel; level++) {
    
    ngrid = ngrid1<<level;        /* Grid dimension */
    scale = (HI - LO) / ngrid;

    /* Partition grid between processors */

    old_ncols = ncols; old_nrows = nrows;        /* Save info for interp. */
    old_col_low = col_low; old_row_low = row_low;
    Partition(ngrid, &ncols, &nrows);

    if (level == 0) {
      grid = Make2d(ncols+1, nrows+1);    /* Allocate first grid */
      Initialize(grid, ncols, nrows);     /* Intitial guess */
    }
    else {
      new_grid = Make2d(ncols+1, nrows+1);
      Interpolate(grid, old_ncols, old_nrows, old_col_low, old_row_low,
		  new_grid, ncols, nrows, col_low, row_low);
      grid = new_grid;
    }

    Solve(grid, ncols, nrows, ngrid, niter, nprint, thresh);
  }

#ifdef PLOT
  ClosePlotFile();
#endif

  cs_total = MTIME_() - start;
  if (NODEID_() == 0) {
    (void) printf("\n  Plot    Exchange   Global-sum   Interpolate    Total ");
    (void) printf("\n ------   --------   ----------   -----------   -------");
    (void) printf("\n %6d    %6d     %6d       %6d      %7d\n",
		  cs_plot, cs_exchange, cs_global,
		  cs_interpolate, cs_total);
  }
		  
  PEND_();                               /* Terminate parallel env. */
  return 0;
}