// psimeq_dbg.java  solve simultaneous equations AX=Y 
// solve real linear equations for X where Y = A * X
// method: Gauss-Jordan elimination using maximum pivot
// usage:  psimeq_dbg(A[][],Y[],X[]); or  psimeq_dbg(n,A[][],Y[],X[]);  or
//                                        psimeq_dbg(n,A[],Y[],X[]);
//    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, then made parallel October 2008

import java.util.concurrent.*; // uses barrier.await() and others

public class psimeq_dbg // for threads on a shared memory architecture
{
  final int NP=4; // number of processors, modify to suit
  int n;          // number of equations
  int m;          // number of equations plus 1
  double B[][];   // working array, various processors on various parts
  int prow[];     // pivot row if not -1, variable number
  int frow[];     // finished row if not 0
  int k_col=0;    // column being reduced, initial zero for "start"
  int k_row=0;    // pivot row
  int srow[] = new int[NP+1]; // starting row for each thread
  double smax_pivot[] = new double[NP]; // max pivot each thread
  int spivot[] = new int[NP]; // row for max pivot each thread
  CyclicBarrier barrier;
  long t_start = System.currentTimeMillis();
  double t_end;
  MyThread sthread[] = new MyThread[NP]; // allocate all threads
 
  psimeq_dbg(final double A[][], final double Y[], double X[])
  {
    n=A.length;
    m=n+1;
    B=new double[n][m];  // working matrix
    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in Matrix.solve, inconsistent array sizes.");
    }
    solve(A, Y, X);
  }

  psimeq_dbg(int nn, final double A[][], final double Y[], double X[])
  {
    n=nn;
    m=n+1;
    B=new double[n][m];  // working matrix
    solve(A, Y, X);
  }

  psimeq_dbg(int nn, final double AA[], final double Y[], double X[])
  {
    n=nn;
    m=n+1;
    double A[][]=new double[n][n];  // reformat
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
	A[i][j] = AA[i*n+j];
    B=new double[n][m];  // working matrix
    solve(A, Y, X);
  }

  void solve(final double A[][], final double Y[], double X[])
  {
    int hold , I_pivot;     // pivot indicies
    double pivot;           // pivot element value
    int part=n/NP;          // number of rows each thread processes

    t_end = (double)(System.currentTimeMillis()-t_start)/1000.0;
    System.out.println("solve starting at   "+t_end);
    // debug print, massively
    System.out.println("psimeq_dbg solve inout data, n="+n+", m="+m);
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
	System.out.println("A["+i+"]["+j+"]="+A[i][j]);
      }
      System.out.println("Y["+i+"]="+Y[i]);
    }
    // set up row range for each thread
    prow=new int[n];     // pivot row
    frow=new int[n];     // finished row
    srow[0]=0;
    for(int i=1; i<NP; i++)
    {
      srow[i]=srow[i-1]+part;
      System.out.println("thread "+0+" start="+srow[i-1]+", <"+srow[i]);
    }
    srow[NP]=n; // may be less than  part 
    System.out.println("thread "+(NP-1)+" start="+srow[NP-1]+", <"+srow[NP]);

    // set up pivot row and finished row
    for(int k=0; k<n; k++)
    {
      prow[k] = -1;
      frow[k] = 0;
    }

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

    // threads build working data structures
    System.out.println("new MyThread's constructed");
    try
    {
      barrier = new CyclicBarrier(NP+1); // st
      System.out.println("master starting to run threads");
      for(int i=0; i<NP; i++) sthread[i].start();
      System.out.println("master waiting for first barrier, st");
      barrier.await(); // st
      for(k_col=0; k_col<n; k_col++)
      {
        barrier = new CyclicBarrier(NP+1); // p1
        System.out.println("at master p1 barrier.await(), k_col="+k_col);
        barrier.await(); // p1
        System.out.println("master passed p1, k_col="+k_col);
        update_row(); 
        t_end = (double)(System.currentTimeMillis()-t_start)/1000.0;
        System.out.println(k_row+" row divided      "+t_end);
        if(k_col>=n-1) break;
        barrier = new CyclicBarrier(NP+1); // p2
        System.out.println("master at p2 barrier.await(), k_col="+k_col);
        barrier.await(); // p2
        System.out.println("master passed p2, k_col="+k_col);
        System.out.println("looping to next k_col");
      } // finished solve
    }
    catch(InterruptedException e)
    {
      System.out.println("InterruptedException in master");
    }
    catch(BrokenBarrierException e)
    {
      System.out.println("BrokenBarrierException in master");
    }

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

    System.out.println("master build X for return");
    //  build  X  for return, unscrambling rows
    for(int i=0; i<n; i++)
    {
      X[i] = B[prow[i]][n];
      System.out.println("prow["+i+"]="+prow[i]);
    }

    t_end = (double)(System.currentTimeMillis()-t_start)/1000.0;
    System.out.println("solve constructor end, took "+t_end);
  } // end solve constructor, solution computed
  
  class MyThread extends Thread
  {
    int myid;
    double my_start;
    
    MyThread(int id, double A[][], double Y[])
    {
      myid = id;
      my_start = (double)(System.currentTimeMillis()-t_start)/1000.0;
      System.out.println(myid+" thread constructed at "+my_start);
      // 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];
	  System.out.println(myid+" B["+i+"]["+j+"]="+B[i][j]);
        }
        B[i][n] = Y[i];
	System.out.println(myid+" B["+i+"]["+n+"]="+B[i][n]);
      }
      System.out.println(myid+" built working data structure "+srow[myid]+
                         "<"+srow[myid+1]);
    } // end MyThread constructor
    
    public void run()
    {
      double th_start, th_end;
      double pivot;
      int rowi;

      th_start = (double)(System.currentTimeMillis()-t_start)/1000.0;
      System.out.println(myid+" starting at     "+my_start);

      try
      {
        System.out.println(myid+" at st barrier.await()");
        barrier.await(); // st
        sleep(1);
        System.out.println(myid+" passed st barrier");

	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);
            }
          }
          System.out.println(myid+" spivot="+spivot[myid]+", smax="+
                             smax_pivot[myid]);
          System.out.println(myid+" at p1 pivot barrier.await()");
          barrier.await(); // p1
          sleep(10);
          System.out.println(myid+" passed p1 pivot barrier");
          
          // 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];
	        System.out.println(myid+" reduce B["+i+
                                   "]["+j+"]="+B[i][j]);
              }
            }
          }
          //  finished inner reduction
          System.out.println(myid+" finished reduction for k_col="+k_col);
          System.out.println(myid+" at p2 reduce barrier.await()");
          if(k_col>=n-1) break; // this task may terminate
          barrier.await(); // p2
          sleep(10);

          System.out.println(myid+" looping to find max for pivot");
	} // do next column
      } // end try
      catch(InterruptedException e)
      {
	System.out.println(myid+" InterruptedException");
      }
      catch(BrokenBarrierException e)
      {
	System.out.println(myid+" BrokenBarrierException");
      }

      th_end = (double)(System.currentTimeMillis()-t_start)/1000.0;
      System.out.println(myid+" thread ends at         "+th_end);
    } // end Run
  } // end class MyThread

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

    System.out.println(" ");
    System.out.println("in update_row for k_col="+k_col);

    // 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];
      }
    }
    System.out.println("update_row I_pivot="+I_pivot+", abs_pivot="+abs_pivot);
    // have pivot, interchange row indicies
    k_row = I_pivot;
    prow[k_col] = k_row;
    frow[k_row] = 1;
    System.out.println("update prow["+k_col+"]="+k_row);
    // 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];
	System.out.println("pivrow B["+k_row+"]["+j+"]="+
                           B[k_row][j]);
      }
      System.out.println("update reduced about pivot k_row="+k_row);
    }
  } // end update_row

} // end class psimeq_dbg