# psimeq_dbg.py   solve simultaneous equations AX=Y 
# solve real linear equations for X where Y = A * X
# method: Gauss-Jordan elimination using maximum pivot
# usage:  X = psimeq(A,Y)
#    Translated to java by : Jon Squire , 26 March 2003
#    First written by Jon Squire December 1959 for IBM 650, translated to
#    other languages  e.g. Fortran converted to Ada converted to C
#    then converted to java,python, then made parallel November 2008

import threading
import thread
from pybarrier import *
import math

global NP
NP=4
global n
n=4
global bar
bar=barrier(NP+1) # barrier for threads and master

class Worker(threading.Thread): # worker initialization
  def __init__ (self,myid):
    self.myid=myid
    threading.Thread.__init__(self)
    print "worker myid=",
    print self.myid,
    print " init, started"
  def run(self):                # worker running
    for job in range(n):
      print "worker myid=",
      print self.myid,
      print ", job=",
      print job
      print " running, at first barrier"
      bar.wait() #1a  master runs
      bar.wait() #1b  master completed
      print "worker myid=",
      print self.myid,
      print ", job=",
      print job
      print " running, at second barrier"
      bar.wait() #2a  master runs
      bar.wait() #2b  master completed
    # end job loop
    print "worker myid=",
    print self.myid,
    print " ending"
    bar.wait() #3a  master waits for all tasks to end, join

"""
for job in range(0,Njob):
  bar.wait() #1a
  print "master at first barrier, job ",
  print job
  bar.wait() #1b

  bar.wait() #2a
  print "master at second barrier, job ",
  print job
  bar.wait() #2b
# end job loop - must be exactly number of iterations as worker
bar.wait() #3a
print "master ending"
"""

def psimeq(A,Y):
  print ' '
  print 'psimeq_dbg starting'
  NP=4              # number of processors, modify to suit
  n=len(Y)          # number of equations
  bar=barrier(NP+1) # barrier for threads and master
  X=[0.0 for i in range(n)]      # returned array
  prow=[-1 for i in range(n)]    # pivot row if not -1, variable number
  frow=[0 for i in range(n)]     # finished row if not 0
  k_col=0    # column being reduced, initial zero for "start"
  k_row=0    # pivot row
  srow=[0 for i in range(NP+1)]  # starting row for each thread
  smax_pivot=[0.0 for i in range(NP)] # max pivot each thread
  spivot=[0 for i in range(NP)]  # row for max pivot each thread
  part=n/NP          # number of rows each thread processes
  # set up row range for each thread
  for i in range(1,NP):
    srow[i]=srow[i-1]+part
  srow[NP]=n # may be less than  part 
  print 'srow=',
  print srow
  # master control
  for i in range(NP): # start N threads
    Worker(i).start()
  for k_col in range(n): # n
    print 'master at 1a'
    bar.wait() #1a  master runs
    print 'master at 1b'
    bar.wait() #1b  master completed
    print 'master at 2a'
    bar.wait() #2a  master runs
    print 'master at 2b'
    bar.wait() #2b  master completed
    print 'master at 3a'
    bar.wait() #3a  master runs
    print 'psimeq_dbg returns'
  return X


"""
    for(int i=0 i<NP i++) sthread[i] = new MyThread(i,A,Y) # construct

    # threads build working data structures
    try
    {
      barrier = new CyclicBarrier(NP+1) # st
      for(int i=0 i<NP i++) sthread[i].start()
      barrier.await() # st
      for(k_col=0 k_col<n k_col++)
      {
        barrier = new CyclicBarrier(NP+1) # p1
        barrier.await() # p1
        update_row() 

        if(k_col>=n-1) break
        barrier = new CyclicBarrier(NP+1) # p2
        barrier.await() # p2
      } # finished solve
    }
    catch(InterruptedException e)
    {
      System.out.println("InterruptedException in master")
    }
    catch(BrokenBarrierException e)
    {
      System.out.println("BrokenBarrierException in master")
    }

    try
    {
      for(int i=0 i<NP i++) sthread[i].join()
    }
    catch(InterruptedException e)
    {
      System.out.println("InterruptedException in master")
    }

    #  build  X  for return, unscrambling rows
    for(int i=0 i<n i++)
    {
      X[i] = B[prow[i]][n]
    }
  } # end solve constructor, solution computed
  
  class MyThread extends Thread
  {
    int myid
    
    MyThread(int id, double A[][], double Y[])
    {
      myid = id
      # build working data structure
      for(int i=srow[myid] i<srow[myid+1] i++)
      {
        for(int j=0 j<n j++)
        {
          B[i][j] = A[i][j]
        }
        B[i][n] = Y[i]
      }
    } # end MyThread constructor
    
    public void run()
    {
      double pivot
      int rowi

      try
      {
        barrier.await() # st
        sleep(1)
	while(true)
	{
	  # find max pivot in my range
          # set smax_pivot[myid]
          # set spivot[myid]
          spivot[myid] = srow[myid]
          smax_pivot[myid] = 0.0 # singular caught elsewhere
          for(int i=srow[myid] i<srow[myid+1] i++)
          {
            if(frow[i]==0 && Math.abs(B[i][k_col]) > smax_pivot[myid])
            {
              spivot[myid] = i
              pivot = B[i][k_col]
              smax_pivot[myid] = Math.abs(pivot)
            }
          }

          barrier.await() # p1
          sleep(10)
          
          # reduce local rows
          # from   srow[myid] to  < srow[myid+1] 
          # using  k_row, k_col
 
          #  inner reduction loop
          for(int i=srow[myid] i<srow[myid+1] i++)
          {
	    if( i != k_row)
            {
              for(int j=k_col+1 j<n+1 j++)
              {
                B[i][j] = B[i][j] - B[i][k_col] *
                                    B[k_row][j]
              }
            }
          }
          #  finished inner reduction
          if(k_col>=n-1) break # this task may terminate
          barrier.await() # p2
          sleep(10)

	} # do next column
      } # end try
      catch(InterruptedException e)
      {
	System.out.println(myid+" InterruptedException")
      }
      catch(BrokenBarrierException e)
      {
	System.out.println(myid+" BrokenBarrierException")
      }
    } # end Run
  } # end class MyThread

  void update_row() # unify max pivots, update row[], 
  {
    int I_pivot, hold
    double abs_pivot

    # find max of threads
    abs_pivot = smax_pivot[0]
    I_pivot = spivot[0]
    for(int i=1 i<NP i++)
    {
      if(smax_pivot[i]>abs_pivot)
      {
	I_pivot = spivot[i]
	abs_pivot = smax_pivot[i]
      }
    }

    # have pivot, interchange row indicies
    k_row = I_pivot
    prow[k_col] = k_row
    frow[k_row] = 1

    # check for near singular
    if(abs_pivot < 1.0E-10)
    {
      for(int j=k_col+1 j<n+1 j++)
      {
        B[k_row][j] = 0.0
      }
      System.out.println("redundant row (singular) "+k_row)
    } # singular, set row to zero, including solution
    else
    {
      # reduce about pivot
      for(int j=k_col+1 j<n+1 j++)
      {
	B[k_row][j] = B[k_row][j] / B[k_row][k_col]
      }
    }
  } # end update_row

} # end class psimeq
"""