#include <stdio.h>
#include <shmem.h>
#ifdef _OPENMP
#include <omp.h>
#endif
	static float **dp;
	void initval( float *p[], float *rhs[], int nrows, int ncols,
		      float pborder, float pinside, float rhsborder, float rhsinside )
{
	int i, j;
/***/
	for ( i = 0; i < nrows; i++ )
	 {
	  for ( j = 0; j < ncols; j++ )
	   {
	    if ( (i == 0) || (j == 0) || 
	         (i == (nrows-1)) || (j == (ncols-1)) )
	     {
	      p[i][j] = pborder;
	      rhs[i][j] = rhsborder;
	     }	 
	    else
	     {
	      p[i][j] = pinside;
	      rhs[i][j] = rhsinside;
	     } 
	   }
	 } 
}
	void allocate( int nrows, int ncols, float **cop, float **corhs )
{
	int i;
/***/
	dp = emalloc( nrows*sizeof(*dp) );
	if ( !dp )
	 {
	  printf( "Cannot allocate\n" );
	 }
	for ( i = 0; i < nrows; i++ )
	 {
	  dp[i] = emalloc( ncols*sizeof(**dp) );
	 if ( !dp[i] )
	  {
	   printf( "Cannot allocate\n" );
	  }
	 }
} 
	void upload( float *p[], float *rhs[], int nrows, int ncols,
		      float *cop, float *corhs )
{
} 
	void download( float *p[], float *rhs[], int nrows, int ncols,
		      float *cop, float *corhs )
{
} 
	void poissn( float *p[], float *rhs[], float **pcop, int icop, float *corhs, 
	             int nrows, int ncols, float delx, float dely, float addcoeff[] )
{
	int i, j, k;
	float rdx2, rdy2, beta2, result;
/***/
	rdx2 = 1./delx/delx;
	rdy2 = 1./dely/dely;
	beta2 = 1.0/(2.0*(rdx2+rdy2));
#pragma omp parallel for schedule(static) num_threads( 7 )
        for ( i = 1; i < (nrows); i++ )
	 {
	  for ( j = 1; j < (ncols-1); j++ )
	   {
	    dp[i][j]  = (beta2*(
		    (p[i+1][j]+p[i-1][j])*rdx2
		  + (p[i][j+1]+p[i][j-1])*rdy2
		  - rhs[i][j])) - p[i][j];    
	   }
	 }
#pragma omp parallel for schedule(static) num_threads( 7 )
        for ( i = 1; i < (nrows-1); i++ )
	 {
	  for ( j = 1; j < (ncols-1); j++ )
	   {
	    p[i][j]  += dp[i][j];
	   }
	 }
}