/*--------------------------------------------------------------------
  NAS Parallel Benchmarks 2.3 UPC version - BT source code tree
  2004 - GWU - HPCL -- check bt.c for complete information
  --------------------------------------------------------------------*/

#include <stdio.h>
#include <upc_relaxed.h>
#include "functions.h"
 
void y_solve_cell( int first, int last, int c )
{
  /*--------------------------------------------------------------------
    --------------------------------------------------------------------*/

  /*--------------------------------------------------------------------
    c performs guassian elimination on this cell
    c 
    c assumes that unpacking routines for non-first cells
    c preload C' and rhs' from previous cell.
    c 
    c assumed send happens outside this routine, but that
    c c'(JMAX) and rhs'(JMAX) will be sent to next cell
    --------------------------------------------------------------------*/
  int i, k, isize, jsize, ksize, jstart;
  int j;

  j = 0;
  jstart = 0;
  isize = cell_size[c][0] - endc[c][0] - 1;
  jsize = cell_size[c][1] - 1;
  ksize = cell_size[c][2] - endc[c][2] - 1;

  /*--------------------------------------------------------------------
    c outer most do loops - sweeping in j direction
    c-------------------------------------------------------------------*/
  if( first == 1 )
    {
      for( i=startc[c][0]; i<=isize; i++ )
	{
	  for( k=startc[c][2]; k<=ksize; k++ )
	    {
	      /*--------------------------------------------------------------------
		c multiply c(i,jstart,k) by b_inverse and copy back to c
		c multiply rhs(jstart) by b_inverse(jstart) and copy to rhs
		c-------------------------------------------------------------------*/
	      binvcrhs( &lhs_priv_d(c,i+1,jstart+1,k+1,BB,0,0),
			&lhs_priv_d(c,i+1,jstart+1,k+1,CC,0,0),
			&rhs_priv_d(c,i+1,jstart+1,k+1,0) );
	    }
	}
    }

  /*--------------------------------------------------------------------
    c begin inner most do loop
    c do all the elements of the cell unless last
    c-------------------------------------------------------------------*/
  for( i=startc[c][0]; i<=isize; i++ )
    {
      for( j=jstart+first; j<jsize-last; j++ )
	{
	  for( k=startc[c][2]; k<=ksize; k++ )
	    {
	      /*--------------------------------------------------------------------
		c rhs(j) -= A*rhs(j-1)
		c-------------------------------------------------------------------*/
	      matvec_sub( &lhs_priv_d(c,i+1,j+1,k+1,AA,0,0),
			  &rhs_priv_d(c,i+1,j-1+1,k+1,0),
			  &rhs_priv_d(c,i+1,j+1,k+1,0) );

	      /*--------------------------------------------------------------------
		c B(j) -= C(j-1)*A(j)
		c-------------------------------------------------------------------*/
	      matmul_sub( &lhs_priv_d(c,i+1,j+1,k+1,AA,0,0),
			  &lhs_priv_d(c,i+1,j-1+1,k+1,CC,0,0),
			  &lhs_priv_d(c,i+1,j+1,k+1,BB,0,0) );

	      /*--------------------------------------------------------------------
		c multiply c(i,j,k) by b_inverse and copy back to c
		c multiply rhs(i,j,k) by b_inverse(i,j,k) and copy to rhs
		c-------------------------------------------------------------------*/
	      binvcrhs( &lhs_priv_d(c,i+1,j+1,k+1,BB,0,0),
			&lhs_priv_d(c,i+1,j+1,k+1,CC,0,0),
			&rhs_priv_d(c,i+1,j+1,k+1,0) );
	    }
	}
    }

  if( first == 0 )
    upc_wait;

  for( i=startc[c][0]; i<=isize; i++ )
    {
      for( k=startc[c][2]; k<=ksize; k++ )
	{
	  /*--------------------------------------------------------------------
	    c rhs(j) -= A*rhs(j-1)
	    c-------------------------------------------------------------------*/
	  matvec_sub( &lhs_priv_d(c,i+1,j+1,k+1,AA,0,0),
		      &rhs_priv_d(c,i+1,j-1+1,k+1,0),
		      &rhs_priv_d(c,i+1,j+1,k+1,0) );

	  /*--------------------------------------------------------------------
	    c B(j) -= C(j-1)*A(j)
	    c-------------------------------------------------------------------*/
	  matmul_sub( &lhs_priv_d(c,i+1,j+1,k+1,AA,0,0),
		      &lhs_priv_d(c,i+1,j-1+1,k+1,CC,0,0),
		      &lhs_priv_d(c,i+1,j+1,k+1,BB,0,0) );

	  /*--------------------------------------------------------------------
	    c multiply c(i,j,k) by b_inverse and copy back to c
	    c multiply rhs(i,j,k) by b_inverse(i,j,k) and copy to rhs
	    c-------------------------------------------------------------------*/
	  binvcrhs( &lhs_priv_d(c,i+1,j+1,k+1,BB,0,0),
		    &lhs_priv_d(c,i+1,j+1,k+1,CC,0,0),
		    &rhs_priv_d(c,i+1,j+1,k+1,0) );
	}
    }
  j++;
  
  /*--------------------------------------------------------------------
    c Now finish up special cases for last cell
    c-------------------------------------------------------------------*/
  if( last == 1 )
    {
      for( i=startc[c][0]; i<=jsize; i++ )
	{
	  for( k=startc[c][2]; k<=ksize; k++ )
	    {  
	      /*--------------------------------------------------------------------
		c rhs(jsize) -= A*rhs(jsize-1)
		c-------------------------------------------------------------------*/
	      matvec_sub( &lhs_priv_d(c,i+1,jsize+1,k+1,AA,0,0),
			  &rhs_priv_d(c,i+1,jsize-1+1,k+1,0),
			  &rhs_priv_d(c,i+1,jsize+1,k+1,0) );

	      /*--------------------------------------------------------------------
		c B(jsize) -= C(jsize-1)*A(jsize)
		c-------------------------------------------------------------------*/
	      matmul_sub( &lhs_priv_d(c,i+1,jsize+1,k+1,AA,0,0),
			  &lhs_priv_d(c,i+1,jsize-1+1,k+1,CC,0,0),
			  &lhs_priv_d(c,i+1,jsize+1,k+1,BB,0,0) );

	      /*--------------------------------------------------------------------
		c multiply rhs() by b_inverse() and copy to rhs
		c-------------------------------------------------------------------*/
	      binvrhs( &lhs_priv_d(c,i+1,j+1,k+1,BB,0,0),
		       &rhs_priv_d(c,i+1,j+1,k+1,0) );
	    }
	}
    }
}