using System;

public class HPLAlgorithm {
 movable hplDistributed( int c,
   double[,] a, double[] b, int n, int p, int q,
   int yStart, int yEnd, int xStart, int xEnd,
   BDChannel top, BDChannel left, BDChannel current, BDChannel bottom, BDChannel right, 
   Channel(double[], int) xChannel ) {
  double[,] l = new double [yEnd - yStart, xEnd];
  double[,] u = new double [yEnd, xEnd - xStart];
  double[] y = new double [yEnd - yStart];
  double[] ySum = new double [yEnd - yStart];
  double[] x = new double [yEnd - yStart];
  int i = 0, j = 0, k = 0; // counters

  if ( b != null ) // if it is the first column in the row then copy part of vector b to ySum
   for ( i = 0; i < b.Length; i++ ) ySum [i] = b [i];

  // Phase 1: Calculate vector y
  int nTimes = 0; // How many arrays we can receive from neighbour processes?
  if ( left != null && top != null ) nTimes = 2;
  else if ( left != null || top != null ) nTimes = 1;

  for ( k = 0; k < nTimes; k++ ) {
   // Receiving part of matrixes L or U from left or top processes
   object[] o = current.Receive();
   int t = (int) o [0];
   int h = (int) o [2];
   int w = (int) o [3];
   double[,] prev = (double[,]) o [1];
   if ( t == 0 ) { // Received part of matrix L from the left process
    // Part of matrix L needed already has been calculated -
    // pass it to right processor in the row with the calculated part of matrix y[k]
    ySum = (double[]) o[4];
    if ( xStart > yStart && xEnd > yEnd && right != null )
     right.Send( 0, prev, h, w, ySum ); // "0" here means that value was sent from the left

    // copying prev to l
    for ( i = 0; i < h; i++ )
     for ( j = 0; j < w; j++ ) l [i,j] = prev [i,j];
   }
   else if ( t == 1 ) { // Received part of matrix L from top process
    y = (double[]) o[4];
    // Part of matrix L needed already has been calculated - pass it to bottom processor in the column
    if ( yStart > xStart && yEnd > xEnd && bottom != null ) 
     bottom.Send( 1, prev, h, w, y ); // "1" here means that value was sent from the top process

    // copying prev to u
    for ( i = 0; i < h; i++ )
     for ( j = 0; j < w; j++ ) u [i,j] = prev [i,j];
   }
  }

  // Calculate parts of matrixes L, U
  for ( i = 0; i < yEnd - yStart; i++ ) {
   int max = xEnd;
   if ( max > n ) max = n;
   for ( j = 0; j < max - xStart; j++ ) {
    int iPos = i + yStart;
    int jPos = j + xStart;
    if ( iPos == jPos ) {
     // Diagonal
     u [iPos,j] = 1;
     l [i,jPos] = a [i,j];
     for ( k = 0; k < jPos; k++ ) l [i,jPos] -= l [i,k] * u [k,j];
    }
    else if ( iPos < jPos ) {
     // Above diagonal
     l [i,jPos] = 0;
     u [iPos,j] = a [i,j];
     for ( k = 0; k < iPos; k++ ) u [iPos,j] = u [iPos,j] - l [i,k] * u [k,j];
     u [iPos,j] = u [iPos,j] / l [i,iPos];
    }
    else {
     // Beyond diagonal
     u [iPos,j] = 0;
     l [i,jPos] = a [i,j];
     for ( k = 0; k < jPos; k++ ) l [i,jPos] -= l [i,k] * u [k,j];
    }
    if ( xStart < yStart ) ySum [i] -= l [i, jPos] * y [j];
   }
  }
  
  // Calculating y[i]
  if ( xStart == yStart )
   for ( i = 0; i < yEnd - yStart; i++ ) {
    y[i] = ySum [i];
    for ( j = 0; j < i; j++ ) y [i] -= y[j] * l [i,j+xStart];
    y [i] = y[i] / l [i,i+xStart];
   }

  // Sending L matrix to right channel if it hasn't been sent before + partial sum ySum
  if ( right != null && (xStart <= yStart || xEnd <= yEnd) ) {
   right.Send( 0, l, yEnd - yStart, xEnd, ySum ); // "0" means that value was sent from the left process
  }

  // Sending U matrix to bottom channel if it hasn't been sent before + part of vector y
  if ( bottom != null && yStart <= xStart )
   bottom.Send( 1, u, yEnd, xEnd - xStart, y ); // "1" means that value was sent from the top process

  // STEP 2: Backward substitution - vector x calculation
  nTimes = 0; 
  if ( xStart == yStart && bottom != null && right != null ) nTimes = 1;
  else if ( xStart > yStart && (bottom != null && right != null) ) nTimes = 2;
  else if ( xStart > yStart && (bottom != null || right != null) ) nTimes = 1;

  for ( i = 0; i < yEnd - yStart; i++ ) ySum [i] = 0;
  for ( i = 0; i < nTimes; i++ ) {
   object[] o = current.Receive();
   int t = (int) o [0];
   if ( t == 2 ) { // Received part of vector x from the bottom process
    x = (double[]) o [1];
    if ( top != null ) top.Send( 2, x );
   }
   else if ( t == 3 ) // Received partial sum of vector x from the right process
    ySum = (double[]) o [1];
  }

  // Calculate x vector and pass it to the top channel if necessary
  if ( xStart == yStart ) {
   for ( i = yEnd - yStart - 1; i >= 0; i-- ) {
    x [i] = y [i] - ySum [i];
    for ( j = xEnd - xStart - 1; j > i; j-- ) x [i] -= x[j] * u[i + yStart, j];
    x [i] = x [i] / u [i +yStart, i];
   }
   if ( top != null ) top.Send( 2, x ); // "2" means that the value was sent from the bottom process
  }
  else if ( xStart > yStart ) {
   if ( left != null ) {
    for ( i = yEnd - yStart - 1; i >= 0; i-- )
     for ( j = xEnd - xStart - 1; j >= 0; j-- ) {
      ySum [i] += u [i,j] * x [j];
     }
    left.Send( 3, ySum ); // "3" means that the value was sent from the right process
   }
  }
  if ( xStart == yStart ) xChannel.Send( x, yStart );
 }

 public double[] Solve( double[,] a, double[] b, int n, int p, int q ) {
  BDChannel[,] bdc = new BDChannel [p,q];
  // Creating n * n bdchannels, where n - number of processors
  for ( int i = 0; i < p; i++ )
   for ( int j = 0; j < q; j++ ) bdc [i, j] = new BDChannel();
  int diffP = n / p;
  int diffQ = n / q;
  int c = 0; // just needed for console output - remove this variable later...
  for ( int i = 0; i < p; i++ ) {
   int h = diffP;
   if ( i == p - 1 ) h += n % p;
   double[] partOfB = GetSubVector( b, i * diffP, i * diffP + h );

   for ( int j = 0; j < q; j++ ) {
    int w = diffQ;
    if ( j == q - 1 ) w += n % q;
    double[,] partOfA = GetSubMatrix( a, i * diffP, i * diffP + h, j * diffQ, j * diffQ + w );
    BDChannel top = null, left = null, bottom = null, right = null;
    if ( i != 0 ) top = bdc [i - 1, j];
    if ( j != 0 ) left = bdc [i, j - 1];
    if ( i < p - 1 ) bottom = bdc [i + 1, j];
    if ( j < q - 1 ) right = bdc [i, j+1];
    c++;
    if ( j != 0 ) partOfB = null; // we need vector b only in the first in processor column
    hplDistributed( c, partOfA, partOfB, n, p, q, i * diffP, i * diffP + h, j * diffQ, j * diffQ + w, top, left, bdc [i,j], bottom, right, xChannel );
   }
  }
  double[] x = new double [n];
  for ( int i = 0; i < p; i++ ) Get( ref x );
  return x;
 }

 void Get( ref double[] x ) & Channel xChannel( double[] partOfX, int wFrom ) {
  for ( int i = 0; i < partOfX.Length; i++ ) x [wFrom + i] = partOfX [i];
 }

 public static void Main( string[] args ) {
  if ( args.Length < 3 ) {
   Console.WriteLine( "Usage: HPL.exe n p q" );
   Console.WriteLine( "Where n - size of matrix A, p - height of process grid, q - width of process grid" );
   return;
  }
  int n = Int32.Parse( args [0] );
  int p = Int32.Parse( args [1] );
  int q = Int32.Parse( args [2] );
  double[,] a = new double [n, n];
  double[] b = new double [n];
//  int maxRandNum = 100;

  // Generate matrix A - Expected mean must be equal to zero
  Random rand = new Random();
  for ( int i = 0; i < n; i++ )
   for ( int j = 0; j < n; j++ ) a [i,j] = (i + j + 1) / (j + 1); //rand.NextDouble() * maxRandNum - maxRandNum/2;

  // Generate vector b - Expected mean must be equal to zero
  for ( int i = 0; i < n; i++ ) b [i] = 1 << i; //rand.NextDouble() * maxRandNum - maxRandNum/2;

  HPLAlgorithm hpl = new HPLAlgorithm();
  double[] x = hpl.Solve( a, b, n, p, q );
 }

 public static double[,] GetSubMatrix( double[,] a, int fromH, int toH, int fromW, int toW ) {
  int h = toH - fromH;
  int w = toW - fromW;
  double[,] b = new double[h, w];
  for ( int i = 0; i < h; i++ )
   for ( int j = 0; j < w; j++ ) b [i, j] = a [fromH+i, fromW+j];
  return b;
 }

 public static double[] GetSubVector( double[] b, int fromH, int toH ) {
  int h = toH - fromH;
  double[] c = new double[h];
  for ( int i = 0; i < h; i++ ) c [i] = b [fromH + i];
  return c;
 }

 public static void PrintMatrix (double[,] a, int h, int w)
 {
  for ( int i = 0; i < h; i++ ) {
   for ( int j = 0; j < w; j++ ) Console.Write( "\t{0:0.00}", a[i, j] );
   Console.WriteLine();
  }
  Console.WriteLine();
 }

 public static void PrintVector( double[] b ) {
  for ( int i = 0; i < b.Length; i++ ) Console.Write( "\t{0:0.00}", b[i] );
  Console.WriteLine();
 }
}