// simeq_thread.java    multiple threads with await()
//                      construct with desired number of threads
//                      then execute with choice of simeq arguments
// solve real linear equations for X where Y = A * X
// method: Gauss-Jordan elimination using maximum pivot
// usage:  simeq(A,Y,X); or  simeq(n,A,Y,X)
//    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, then threads August 2009
//
// use:   simeq_thread ST = new simeq_thread(8);  number of threads
//        ST.simeq(A,X,Y);  or ST.simeq(n,A,X,Y) etc.

import java.util.concurrent.BrokenBarrierException;
import java.util.concurrent.CyclicBarrier;
 
public class simeq_thread
{
  int nthreads;    // obviously, the number of threads
  Thread t[];      // threads
  double pivots[]; // index is thread id, largest diagonal
  int ipivots[];   // index is thread id, index of largest diagonal
  int i1[];        // index is thread id, starting row in row[]
  int i2[];        // index is thread id, ending row+1 in row[]
  CyclicBarrier barrier; // for nthreads

  int n;           // number of linear equations
  double B[][];    // n rows, n+1 columns, column n+1 is Y
  double X[];      // computed solution
  int row[];       // index of permuted row interchange
  double gpivot;   // global pivot value
  int igpivot;     // index of row that has global pivot

  simeq_thread(int nthreads)
  {
    System.out.println("simeq_thread running with nthreads="+nthreads);
    this.nthreads = nthreads;
    t = new Thread[nthreads];
    pivots = new double[nthreads];
    ipivots = new int[nthreads];
    i1 = new int[nthreads];
    i2 = new int[nthreads];
  } // end simeq_thread that set number of threads

  public void simeq(final double A[][], final double Y[], double X[])
  {
    n = A.length;
    B = new double[n][n+1];  // working matrix
    this.X = X;
    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in simeq, inconsistent array sizes.");
    }
    // build working data structure
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = A[i][j];
      B[i][n] = Y[i];
    }
    solve();
  } // end simeq(A,Y,X)

  public void simeq(int n, final double A[][], final double Y[], double X[])
  {
    this.n = n;
    B = new double[n][n+1];  // working matrix
    this.X = X;
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = A[i][j];
      B[i][n] = Y[i];
    }
    solve();
  } // end simeq(n,A,Y,X)

  public void simeq(int n, final double AA[], final double Y[], double X[])
  {
    this.n = n;
    B = new double[n][n+1];  // working matrix
    this.X = X;
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = AA[i*n+j];
      B[i][n] = Y[i];
    }
    solve();
  } // end simeq(n,AA,Y,X)

  public void psimeq(int NP, int n, final double AA[], final double Y[],
                     double X[])
  {
    this.nthreads = NP;
    this.n = n;
    B = new double[n][n+1];  // working matrix
    this.X = X;
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = AA[i*n+j];
      B[i][n] = Y[i];
    }
    solve();
  } // end psimeq(NP,n,AA,Y,X)

  public void simeq(int n, double B[][], double X[]) // B destroyed
  {
    this.n = n;
    this.B = B;
    this.X = X;
    solve(); // argument B[n][n+1] passed
  } // end simeq(n,B,X)

  private void solve() // n, B, X are set up
  {
    int j = 0;
    int inc = (int)Math.ceil((double)n/(double)nthreads);

    try
    {
      if(n<nthreads) nthreads=n;
      row = new int[n];
      for(int i=0; i<nthreads; i++)
      {
        i1[i] = j;                 // starting row
        if(j+inc>n) i2[i] = n;
        else        i2[i] = j+inc; // < ending row
        j = j+inc;
        if((nthreads-i-1)*(inc-1)+j>=n) inc--; // at most once and last
      }
      // create and start the threads
      Interchanger interchange = new Interchanger();

      barrier = new CyclicBarrier(nthreads, interchange);
      for(int i=0; i<nthreads; i++)
      {
        t[i] = new tsolve(i);
        t[i].start(); // i is thread id
      }
      for(int i=0; i<nthreads; i++)
      {
        t[i].join(); // wait for all to complete
      }
    }
    catch (InterruptedException ex)
    {
      System.out.println("InterruptedException");
    }
  } // end solve

  public class tsolve extends Thread
  { 
    int id;
    int I_pivot;          // pivot indicies
    double abs_pivot;     // gpivot will be global

    tsolve(int id)
    {
      this.id = id;
    } // end tsolve constructor

    public void run()
    {
      try // needed by await()
      {
        //  begin main reduction loop across columns
	for(int k=0; k<n; k++)
        {
          if(k==0)
          {
            // set up row interchange vectors
            for(int i=i1[id]; i<i2[id]; i++)
            {
              row[i] = i;
            }
          } // end one time setup
          // find largest element for pivot
          abs_pivot = -1.0; // not set by this thread
          I_pivot = -1;
          for(int i=i1[id]; i<i2[id]; i++)
          {
            if(i>=k && Math.abs(B[row[i]][k]) > abs_pivot)
            {
              I_pivot = i;
              abs_pivot = Math.abs(B[row[i]][k]);
            }
          }
          // have pivot, interchange row indicies
          pivots[id] = abs_pivot;
          ipivots[id] = I_pivot;
          barrier.await(); // let interchanger do interchange
          //  inner reduction loop
          for(int i=i1[id]; i<i2[id]; i++)
          {
            if(i != k)
            {
              for(int j=k+1; j<n+1; j++)
              {
                B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j];
              }
            }
          } //  finished inner reduction
        } //  end main reduction loop
        pivots[id] = -1.0;
        ipivots[id] = -1; // should never be used
        barrier.await(); // let interchange finish up
      } // try
      catch (InterruptedException ex)
      {
        System.out.println("InterruptedException");
      }
      catch (BrokenBarrierException ex)
      {
        System.out.println("BrokenBarrierException");
      } // end try and catch
    } // end tsolve run()
  } // end class tsolve

  class Interchanger implements Runnable // every time all reach await()
  {
    int hold;
    int k = -1;

    Interchanger()
    {
      // nothing to do
    }

    public void run()
    {
      k++; // pivot column
      if(k<n) // all except last for cleanup
      {
        igpivot = -2;
        gpivot = -2.0;
        for(int i=0; i<nthreads; i++)
        {
	  if(pivots[i]>gpivot)
          {
	    igpivot =ipivots[i]; // from tsolve threads
	    gpivot = pivots[i];  // from tsolve threads
          }
        }
        if(igpivot<0 || gpivot < 0.0)
	{
	  System.out.println("simeq_thread.java internal error");
          System.out.println("n="+n+", nthreads="+nthreads);
          System.out.println("k="+k+", igpivot="+igpivot+", gpivot="+gpivot);
          // System.exit(1);
	}
        hold = row[k];
        row[k] = row[igpivot];
        row[igpivot] = hold;
        // check for near singular
        if(gpivot < 1.0E-15)
        {
	  for(int j=k+1; j<n+1; j++)
	  {
	    B[row[k]][j] = 0.0;
          }
          B[row[k]][k] = 1.0;
        } // singular row zeroed
        // reduce about pivot
        for(int j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
        }
      }
      else if(k==n)
      {
        //  build  X  for return, unscrambling rows
        for(int i=0; i<n; i++)
        {
          X[i] = B[row[i]][n];
        }
      } // end if(k<n)
      else
      {
	  System.out.println("k>n simeq_thread.java internal error");
      }
    } // end interchanger run()
  } // end class Interchanger
} // end class simeq_thread