# Gauss-Seidel                          Qiang A. Zhao, October 92
#
# Usage: a.out Size Workers North East South West < Color-file
#
#        Color-file:
#               <# of colors>
#               <r> <g> <b>
#               ...
#               <from violet to red>
#               ...

resource gauss()

  import MPDWin

  ### Var declarations/initializations

  const real EPSILON = 0.0001
  const int MAXCV = 65535
  const int MINS = 100
  const int MINWK = 2
  const int MINC = 2
  const int SCALE = 5
  const int TIMES = 100

  int  SIZE = MINS;             getarg(1, SIZE)
  int  WORKERS = MINWK;         getarg(2, WORKERS)
  real North = 0.0;             getarg(3, North)
  real East = 0.333;            getarg(4, East)
  real South = 0.666;           getarg(5, South)
  real West = 1.0;              getarg(6, West)
  int  passes = 0

  if ((SIZE < MINS) or (WORKERS < MINWK)) {
    write(stderr, "Invalid arguments")
    stop(1)
  }

  if (SIZE mod WORKERS != 0) {
    write(stderr, "Number of workers must evenly divide image size.")
    stop(1)
  }

  const int  STRIP = SIZE / WORKERS
  int  COLORS = MINC;
  bool stay[WORKERS] = ([WORKERS]true)
  bool changed[WORKERS] = ([WORKERS]false)
  sem  done = 0, continue[WORKERS] = ([WORKERS]0)
  real image[0:SIZE+1, 0:SIZE+1]

  winWindow  mainwin
  winImage  wim 
  op winEventChannel ec
  winEvent  ev

  ### open a window, create an image

  mainwin = WinOpen("", "Gauss-Seidel", ec, UseDefault, SIZE*SCALE, SIZE*SCALE)
  if (mainwin == null) {
    write(stderr, "Ouch, can't open window")
    stop(1)
  }
  WinSetEventMask(mainwin, Ev_KeyUp | Ev_ButtonUp)
  wim = WinCreateImage(mainwin, UseDefault, SIZE*SCALE, SIZE*SCALE)
  if (wim == null) {
    write(stderr, "Ouch, can't create image")
    stop(1)
  }
  WinFlush(mainwin)

  # how many colors?
  read(COLORS)
  if (COLORS < MINC) {
    write(stderr, "Invalid number of colors")
    WinClose(mainwin)
    stop(1)
  }

  # additional vars
  winWindow cwins[0:COLORS-1]
  winPixel cpvs[0:COLORS-1]

  # initialize the "constant" part of the image
  for [ i=0 to SIZE+1 ] {
    image[i, 0] = North
    image[i, SIZE+1] = South
    image[0, i] = West
    image[SIZE+1, i] = East
  }

  # initialize colors
  const int  CSTRIP = SIZE*SCALE/COLORS
  for [ i= 0 to COLORS-1 ] {
    cwins[i] = WinNewContext(mainwin)
    if (cwins[i] == null) {
      write(stderr, "No enough memory... give up")
      WinClose(mainwin)
      stop(1)
    }

    ### grey-scale
    # int cv = MAXCV * (i+1) / COLORS
    # string[5] str
    # sprintf(str, "%04x", cv)
    # WinSetForeground(cwins[i], "#" || str || str || str)

    ### read-in colors
    int r, g, b
    string[5] sr, sg, sb
    read(r, g, b)
    sprintf(sr, "%04x", r*255)
    sprintf(sg, "%04x", g*255)
    sprintf(sb, "%04x", b*255)
    cpvs[i] = WinSetForeground(cwins[i], "#" || sr || sg || sb)

    WinFillRectangle(cwins[i], winRectangle(CSTRIP*i, 0, CSTRIP, SIZE*SCALE))
  }
  WinFlush(mainwin)

  ### Worker processes to compute averages

  # using checkerboard partitioning
  process worker [ w = 1 to WORKERS ] {
    const int R = (w - 1) * STRIP+1                # starting row of strip
    const int EP = R + STRIP - 1
    real  tmp 

    int  phase = w mod 2
    int  ptmp1, ptmp2, wk, i, j 

    for [ i= R to R+STRIP-1 ] {
      for [ j=1 to SIZE ] {
	image[i, j] = 0.0
      }
    }

    P(continue[w])

    while (stay[w]) {
      for [ wk = phase to WORKERS-1 by 2 ] {
	ptmp1 = wk * STRIP + 1
	ptmp2 = ptmp1 + STRIP - 1
	for [ i= R to EP ] {
	  for [ j= ptmp1 to ptmp2 ] {
	    tmp = (image[i+1, j] + image[i-1, j] +
			image[i, j+1] + image[i, j-1])/4.0
	    if (abs(tmp - image[i, j]) > EPSILON) {
	      changed[w] = true
	      if (passes == 0 or passes == 1) {
	      int x = (i-1)*SCALE, y = (j-1)*SCALE
	      for [ l = 0 to SCALE-1, k = 0 to SCALE-1 ] {
		WinPutPixel(wim, winPoint(x+l, y+k), cpvs[int(tmp*COLORS)])
	      }
	      }
	    }
	    image[i, j] = tmp
	  }
	}
      }        # wk

      V(done)
      phase = (phase+1) mod 2
      P(continue[w])
    }

    V(done)             # signal leave
  }# end of worker

  ### Coordinator

  const char  Rolling[0:3] = ('-', '\\', '|', '/')
  int  end_time 
  int  start_time 

  winPixel  lightYellow = WinSetForeground(mainwin, "LightYellow")
  for [ i = 0 to SIZE * SCALE - 1 ] {
    for [ j = 0 to SIZE * SCALE - 1 ] {
      WinPutPixel(wim, winPoint(i, j), lightYellow)
    }
  }

  writes(stderr, "Press a key in the window to start...")
  receive ec(ev)
  write(stderr)
  start_time = age()
  for [ mm = 1 to WORKERS ] { V(continue[mm]) }

  # wait to quit
  bool  ac = true
  while (ac) {
    # check workers
    for [ k = 1 to WORKERS ] { P(done) }
    #writes(stderr, "\r", Rolling[passes mod 4 + 1])
    if (passes == 1) {
      WinPutImage(mainwin, wim, winRectangle(0, 0, SIZE*SCALE, SIZE*SCALE), winPoint(0, 0))
      WinFlush(mainwin)
    }
    if (ac) {            # check convergence
      ac = false
      for [ k = 1 to WORKERS ] { ac = ac or changed[k]; changed[k] = false }
    }
    if (not ac) { stay = ([WORKERS] false) }
    passes = (passes + 1) mod TIMES
    for [ k = 1 to WORKERS ] { V(continue[k]) }
  }

  end_time = age()
  for [ k = 1 to WORKERS ] { P(done) }
  write(stderr)
  writes(stderr, "Press a key in the window to quit...")
  receive ec(ev)
  stop(0)

  # final code -- close connection...
  final {
    write(stderr)
    write(stderr, "Computation time: ", end_time - start_time)
    WinClose(mainwin)
  }

end gauss