grid.c

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

#include <stdio.h>
#include <math.h>

#ifdef PLOT
#include <sys/types.h>
#include <netinet/in.h>
FILE *plot_file;

#define PLOT_VALUE    1
#define PLOT_ERROR    2
#define PLOT_RESIDUAL 3

static long plot_type = 0;      /* 0 means no plot */
#endif

#include "sndrcv.h"

#if defined(DELTA) || defined(IPSC)
#define htonl(a) (a)
#endif

#if !defined(AIX)
extern char *malloc();
#endif
extern void exit();

#define MAX(a,b) (((a)>(b)) ? (a) : (b))
#define MIN(a,b) (((a)<(b)) ? (a) : (b))
#define ABS(a)   (((a)>=0 ) ? (a) : -(a))

#define LO -3.1415926535       /* Phsyical dimensions   */
#define HI  3.1415926535

double scale;                  /* Physical mesh grain   */

long P, ncol_P, nrow_P;         /* P = No. of processes = ncol_P * nrow_P */
long my_col_P, my_row_P;        /* Co-ords of this proc. in grid of processes*/

double *buffer;                /* Buffer used for exchanging rows */

long north, south, east, west;  /* No. of process in that direction on 
                                  process grid or -1 if no-one there */

long col_low = 0;                /* Grid co-ords of top left corner */
long row_low = 0;

/* Mapping from GLOBAL grid index to physical coordinates */
#define MAPCOL(I) (scale*((double) (I+col_low)) + LO)
#define MAPROW(J) (scale*((double) (J+row_low)) + LO)

long cs_exchange = 0;            /* Timing information */
long cs_global = 0;
long cs_interpolate = 0;
long cs_plot = 0;
long cs_total = 0;

/*
void ZeroGrid(grid, ncols, nrows)
     double **grid;
     long ncols, nrows;
{
  long i,j;

  for (i=0; i<=ncols; i++)
    for (j=0; j<=nrows; j++)
      grid[i][j] = 0.0;
}
*/

double Solution(x,y)
     double x,y;
/*
  Model potential that defines solution and boundary conditions.
  This particular choice makes the discretization error quite apparent.
*/
{
/*
  return cos(x)*exp(-y) + 0.1*sin(2.*x)*exp(-2.*y) + 0.1*cos(3.*x)*exp(-3.*y);
*/

  return 0.5 * (cos(x)*exp(-y) + cos(y)*exp(-x) +
		0.1*(sin(2.*x)*exp(-2.*y) + sin(2.*y)*exp(-2.*x)) +
		0.1*(cos(3.*x)*exp(-3.*y) + cos(3.*y)*exp(-3.*x)));
}

double GridError(grid, ncols, nrows, ngrid)
     long ncols, nrows, ngrid;
     double **grid;
/*
  Compute mean absolute error at grid points relative to
  the analytic solution ... error is due to either lack of
  convergence or discretization error.
*/
{
  long i,j, type=9, length=1;
  long start;
  double error = 0.0;

  for (i=1; i<ncols; i++)
    for (j=1; j<nrows; j++)
      error += fabs(grid[i][j] - Solution(MAPCOL(i), MAPROW(j)));

  start = MTIME_();
  DGOP_(&type, &error, &length, "+");
  cs_global += MTIME_() - start;

  return error / ((ngrid-1)*(ngrid-1));
}
  
void BoundaryConditions(grid, ncols, nrows)
     double **grid;
     long ncols, nrows;
/* 
  Fill in b.c.s on a grid
*/
{
  long i,j;

  for (i=0; i<=ncols; i++) {
    if (my_row_P == 0)
      grid[i][0]  = Solution(MAPCOL(i), LO);
    if (my_row_P == (nrow_P-1))
      grid[i][nrows] = Solution(MAPCOL(i), HI);
  }

  for (j=0; j<=nrows; j++) {
    if (my_col_P == 0)
      grid[0][j]  = Solution(LO, MAPROW(j));
    if (my_col_P == (ncol_P-1))
      grid[ncols][j] = Solution(HI, MAPROW(j));
  }
}

void Initialize(grid, ncols, nrows)
     long ncols, nrows;
     double **grid;
/*
  Fill in boundary values and zero initial guess for interior.
*/
{
  long i,j;

  BoundaryConditions(grid, ncols, nrows);

  for (i=1; i<ncols; i++)
    for (j=1; j<nrows; j++)
      grid[i][j] = 0.0;
}

void Exchange(grid, ncols, nrows)
     double ** grid;
     long ncols, nrows;
/*
  Exchange data with neighboring ndoes. In Operate only need
  to exchange red and black elements separately but do NOT do
  this at the moment ... thus are sending twice as much data
  as necessary. However, are still dominated by the latency
  so 'no big whup'.  Interpolate needs to exchange the full
  boundary information (I think).
*/
{
  long type1=1, type2=2, type3=3, type4=4, type5=5, type6=6, type7=7, type8=8;
  long bncols = (ncols+1)*sizeof(double);
  long bnrows = (nrows+1)*sizeof(double);
  long synch = 1, lenmes, nodefrom, i;
  long start = MTIME_();

#define GATHER(k)  for (i=0; i<=ncols; i++) buffer[i] = grid[i][k]
#define SCATTER(k) for (i=0; i<=ncols; i++) grid[i][k] = buffer[i]

  if (my_col_P%2) {
    if (west >= 0) {
      SND_(&type1, (char *) grid[1], &bnrows, &west, &synch);
      RCV_(&type2, (char *) grid[0], &bnrows, &lenmes, &west, 
	   &nodefrom, &synch);
    }
    if (east >= 0) {
      SND_(&type3, (char *) grid[ncols-1], &bnrows, &east, &synch);
      RCV_(&type4, (char *) grid[ncols], &bnrows, &lenmes, &east, 
	   &nodefrom, &synch);
    }
  }
  else {
    if (east >= 0) {
      RCV_(&type1, (char *) grid[ncols], &bnrows, &lenmes, &east, 
	   &nodefrom, &synch);
      SND_(&type2, (char *) grid[ncols-1], &bnrows, &east, &synch);
    }
    if (west >= 0) {
      RCV_(&type3, (char *) grid[0], &bnrows, &lenmes, &west, 
	   &nodefrom, &synch);
      SND_(&type4, (char *) grid[1], &bnrows, &west, &synch);
    }
  }

  if (my_row_P%2) {
    if (north >= 0) {
      GATHER(1);
      SND_(&type5, (char *) buffer, &bncols, &north, &synch);
      RCV_(&type6, (char *) buffer, &bncols, &lenmes, &north, 
	   &nodefrom, &synch);
      SCATTER(0);
    }
    if (south >= 0) {
      GATHER(nrows-1);
      SND_(&type7, (char *) buffer, &bncols, &south, &synch);
      RCV_(&type8, (char *) buffer, &bncols, &lenmes, &south, 
	   &nodefrom, &synch);
      SCATTER(nrows);
    }
  }
  else {
    if (south >= 0) {
      RCV_(&type5, (char *) buffer, &bncols, &lenmes, &south, 
	   &nodefrom, &synch);
      SCATTER(nrows);
      GATHER(nrows-1);
      SND_(&type6, (char *) buffer, &bncols, &south, &synch);
    }
    if (north >= 0) {
      RCV_(&type7, (char *) buffer, &bncols, &lenmes, &north, 
	   &nodefrom, &synch);
      SCATTER(0);
      GATHER(1);
      SND_(&type8, (char *) buffer, &bncols, &north, &synch);
    }
  }
  cs_exchange += MTIME_() - start;
}

#ifdef PLOT

void ClosePlotFile()
{
  if (!plot_type)
    return;

  if (NODEID_() == 0)
    (void) fclose(plot_file);
}

void OpenPlotFile(maxgrid)
     long maxgrid;
{
  if (!plot_type)
    return;

  if (NODEID_())
    return;              /* Only node 0 needs to do this */

  if (!(plot_file = fopen("plot","w+")))
    Error("OpenPlotFile: failed to open plot file", (long) -1);

  maxgrid = htonl((long) maxgrid); /* For portability */
  if (fwrite((char *) &maxgrid, sizeof(int), 1, plot_file) != 1)
    Error("OpenPlotFile: failed to write maxgrid", (long) -1);

  (void) fflush(plot_file);
}

void InsertPlot(plot_full, ngrid, plot_grid, nrows, ncols, col_low, row_low)
     unsigned char *plot_full, *plot_grid;
     long ngrid, nrows, ncols, col_low, row_low;
{
  long i, j;
  long i_lo = (col_low == 0) ? 0 : 1;
  long j_lo = (row_low == 0) ? 0 : 1;
  unsigned char *temp;

  /* Dink around with i_lo, j_lo so that the interior edges in the
     parallel case are not overwritten with incorrect values */

  if (i_lo)
    plot_grid += nrows;

  for (i=i_lo; i<ncols; i++) {
    temp = plot_full + (i+col_low)*ngrid + row_low + j_lo;
    if (j_lo)
      plot_grid++;
    for (j=j_lo; j<nrows; j++)
      *temp++ = *plot_grid++;
  }
}

void PlotGrid(grid, ngrid, ncols, nrows)
     double **grid;
     long ngrid, ncols, nrows;
{
  unsigned char *plot_grid, *plot_full;

  double d_63_19 = 63.0 / 19.0;
  double rlog2 = 1.0 / log((double) 2.0);
#define LOG2(x)  ((x != 0.0) ? (rlog2 * log(x)) : -1024.0)
#define VALUE(x) (d_63_19 * (LOG2(fabs(x)) + 9.0))
  register long i, j, value;
  register unsigned char *temp;
  long n;
  double residual, factor;
  long net_ngrid = htonl((long) ngrid);
  long msg[4], synch_type=13, anyone = (-1);
  long msg_len=sizeof(msg), msg_type=120, length, node, zero=0, synch=1;
  long start = MTIME_();

  if (!plot_type)
    return;

  n = nrows*ncols;
  
  plot_grid = (unsigned char *) malloc((unsigned) n);
  if (!plot_grid)
    Error("PlotGrid: failed to allocate plot grid", (long) -1);

  temp = plot_grid;

  switch (plot_type) {
  case PLOT_VALUE:
    for (i=0; i<ncols; i++)
      for (j=0; j<nrows; j++) {
	value = VALUE(grid[i][j]);
	value = MIN(value, 63);
	value = MAX(value, 0);
	*temp++ = (unsigned char) value;
      }
    break;

  case PLOT_ERROR:
    for (i=0; i<ncols; i++)
      for (j=0; j<nrows; j++) {
	value = VALUE(grid[i][j] - Solution(MAPCOL(i), MAPROW(j)));
	value = MIN(value, 63);
	value = MAX(value, 0);
	*temp++ = (unsigned char) value;
      }
    break;

  case PLOT_RESIDUAL:

    Exchange(grid, ncols, nrows);  /* This to get the edges up date */

    factor = 1.0/(scale*scale);
    for (j=0; j<nrows; j++)        /* Edge has residual zero */
      *temp++ = (unsigned char) VALUE(0.0);
    for (i=1; i<ncols; i++) {
      *temp++ = (unsigned char) VALUE(0.0);
      for (j=1; j<nrows; j++) {
	residual = grid[i+1][j] + grid[i-1][j] + grid[i][j+1] +
	  grid[i][j-1] - 4.0*grid[i][j];
	value = VALUE(residual*factor);
	value = MIN(value, 63);
	value = MAX(value, 0);
	*temp++ = (unsigned char) value;
      }
    }
    break;
    
  default:
    Error("PlotGrid: unknown plot type requested", (long) plot_type);
  }
  
  /* Now have grid with my portion of the plot. This gets sent
     to node 0 which assembles the full grid before dumping to disc */
  
  if (NODEID_()) {
    msg_type = 120;
    msg[0] = nrows; msg[1] = ncols; msg[2] = col_low; msg[3] = row_low;
    SND_(&msg_type, (char *) msg, &msg_len, &zero, &synch);
    msg_type++;
    SND_(&msg_type, (char *) plot_grid, &n, &zero, &synch);
    msg_type++;
  }
  else {
    plot_full = (unsigned char *) malloc((unsigned) (ngrid*ngrid));
    if (!plot_full)
      Error("PlotGrid: failed to allocate plot full", (long) -1);
    InsertPlot(plot_full, ngrid, plot_grid, nrows, ncols, col_low, row_low);

    for (i=1; i<NNODES_(); i++) {
      msg_type = 120;
      RCV_(&msg_type, (char *) msg, &msg_len, &length, &anyone, &node, &synch);
      msg_type++;
      RCV_(&msg_type, (char *) plot_grid, &n, &length, &node, &node, &synch);
      msg_type++;
      InsertPlot(plot_full, ngrid, plot_grid, msg[0], msg[1], msg[2], msg[3]);
    }

    if (fwrite((char *) &net_ngrid, sizeof(int), 1, plot_file) != 1)
      Error("PlotGrid: failed to write ngrid", (long) -1);
    n = ngrid*ngrid;
    if (fwrite(plot_full, 1, n, plot_file) != n)
      Error("PlotGrid: failed writing to file", (long) n);
    (void) free((char *) plot_full);
    (void) fflush(plot_file);
    (void) fflush(plot_file);
  }
    
  BRDCST_(&synch_type, &zero, &zero, &zero);    /* Synchronization */
  (void) free((char *) plot_grid);

  cs_plot += MTIME_() - start;
}
#endif

void PrintGrid(grid, ncols, nrows)
     double **grid;
     long ncols, nrows;
{
  long i,j,newline;

  for (i=0; i<=ncols; i++) {
    (void) printf("Column %d \n",i+col_low);
    newline = 0;
    for(j=0; j<=nrows; j++) {
      if (grid[i][j] != 0.0) {