#include "Util/XPlatform.h"
#include "Util/MultiStreamIO.h"
#include "Util/Socket.h"
#include "Util/Half/half.h"

#define LOW_MEMORY_HACK 0

#define NEW_MULTIPLY 1

#if NEW_MULTIPLY
template< bool AddAverage >
Vector< double > Multiply( const SparseMatrix< double >& A , const Vector< double >& x )
{
	if( AddAverage )
	{
		Vector< double > y = A*x;
		double average = 0;
		for( size_t i=0 ; i<x.Dimensions() ; i++ ) average +=x [i];
		average /= x.Dimensions();
		for( size_t i=0 ; i<x.Dimensions() ; i++ ) y[i] += average;
		return y;
	}
	else return A*x;
}
template< bool AddAverage >
int MySolveConjugateGradient( const SparseMatrix< double >& A , const Vector< double >& b , const int& iters , Vector< double >& x )
{
	double eps=1e-16;
	int dim = (int)( b.Dimensions() );
	Vector< double > r ,d , Ad;

	double delta_new = 0 , delta_0;

	////////////////////////
	// d = r = b - M * x
	// \delta_new = ||r||^2
	d = r = b - Multiply< AddAverage >( A , x );
	for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i];
	////////////////////////

	delta_0 = delta_new;
	if( delta_new<eps )
	{
		fprintf( stderr , "[WARNING] Initial residual too low: %g < %g\n" , delta_new , eps );
		return 0;
	}
	int ii;
	for( ii=0 ; ii<iters && delta_new>eps*delta_0 ; ii++ )
	{
		////////////////////////////////////
		// \alpha = ||r||^2 / (d^t * M * d)
		Ad = Multiply< AddAverage >( A , d );
		double dDotMd = 0;
		for( int i=0 ; i<dim ; i++ ) dDotMd += Ad[i] * d[i];
		double alpha = double( delta_new / dDotMd );
		double delta_old = delta_new;
		////////////////////////////////////

		delta_new = 0;

		if( (ii%50)==(50-1) )
		{
			////////////////////////////////
			// x = x + d * alpha
			// r = b - M * ( x + d * alpha )
			//   = r - M * d * alpha
			// \delta_new = ||r||^2
			for( int i=0 ; i<dim ; i++ ) x[i] += d[i] * alpha;
			r.SetZero();
			r = b - Multiply< AddAverage >( A , x );
			for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i];
			////////////////////////////////
		}
		else
		{
			////////////////////////////////
			// x = x + d * alpha
			// r = r - M * d * alpha
			// \delta_new = ||r||^2
			r -= Ad * alpha;
			for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i] , x[i] += d[i] * alpha;
			////////////////////////////////
		}
		////////////////////////////////
		// beta = ||r||^2 / ||r_old||^2
		// d = r + d * beta
		double beta = delta_new / delta_old;
		for( int i=0 ; i<dim ; i++ ) d[i] = r[i] + d[i] * beta;
		////////////////////////////////
	}
	return ii;
}
#else // !NEW_MULTIPLY
template< int Type >
Vector< double > Multiply( const SparseMatrix< double >& A , const Vector< double >& x );

template< int Type >
Vector< double > Multiply( const SparseMatrix< double >& A , const Vector< double >& x )
{
	return A*x;
}
template<> inline Vector< double > Multiply< ZERO_DERIVATIVE >( const SparseMatrix<double>& A , const Vector<double>& x )
{
	Vector<double> y = A*x;
	double average=0;
	for(size_t i=0;i<x.Dimensions();i++)	average+=x[i];
	average/=x.Dimensions();
	for(size_t i=0;i<x.Dimensions();i++)	y[i]+=average;
	return y;
}
template<> inline Vector< double > Multiply< PERIODIC_BOUNDARY >( const SparseMatrix<double>& A , const Vector<double>& x )
{
	Vector<double> y = A*x;
	double average=0;
	for(size_t i=0;i<x.Dimensions();i++)	average+=x[i];
	average/=x.Dimensions();
	for(size_t i=0;i<x.Dimensions();i++)	y[i]+=average;
	return y;
}
template< class Type >
Vector< double > Multiply( const SparseMatrix< double >& A , const Vector< double >& x ) { return A*x; }
template< int Type >
int MySolveConjugateGradient( const SparseMatrix< double >& A,const Vector< double >& b , const int& iters , Vector< double >& x )
{
	double eps=1e-16;
	int dim = b.Dimensions();
	Vector< double > r ,d , Ad;

	double delta_new = 0 , delta_0;

	////////////////////////
	// d = r = b - M * x
	// \delta_new = ||r||^2
	d = r = b - A*x;
	for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i];
	////////////////////////

	delta_0 = delta_new;
	if( delta_new<eps )
	{
		fprintf( stderr , "[WARNING] Initial residual too low: %g < %g\n" , delta_new , eps );
		return 0;
	}
	int ii;
	for( ii=0 ; ii<iters && delta_new>eps*delta_0 ; ii++ )
	{
		////////////////////////////////////
		// \alpha = ||r||^2 / (d^t * M * d)
		Ad = A*d;
		double dDotMd = 0;
		for( int i=0 ; i<dim ; i++ ) dDotMd += Ad[i] * d[i];
		double alpha = double( delta_new / dDotMd );
		double delta_old = delta_new;
		////////////////////////////////////

		delta_new = 0;

		if( (ii%50)==(50-1) )
		{
			////////////////////////////////
			// x = x + d * alpha
			// r = b - M * ( x + d * alpha )
			//   = r - M * d * alpha
			// \delta_new = ||r||^2
			for( int i=0 ; i<dim ; i++ ) x[i] += d[i] * alpha;
			r.SetZero();
			r = b - A*x;
			for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i];
			////////////////////////////////
		}
		else
		{
			////////////////////////////////
			// x = x + d * alpha
			// r = r - M * d * alpha
			// \delta_new = ||r||^2
			r -= A * d * alpha;
			for( int i=0 ; i<dim ; i++ ) delta_new += r[i] * r[i] , x[i] += d[i] * alpha;
			////////////////////////////////
		}
		////////////////////////////////
		// beta = ||r||^2 / ||r_old||^2
		// d = r + d * beta
		double beta = delta_new / delta_old;
		for( int i=0 ; i<dim ; i++ ) d[i] = r[i] + d[i] * beta;
		////////////////////////////////
	}
	return ii;
}
#endif // NEW_MULTIPLY

void SetPaddedSize( int width , int height , int minRes , int& paddedWidth , int& paddedHeight , int& widthBlockSize, int& heightBlockSize )
{
	// Compute the size of the domain
	int domainW = FiniteElements1D< float , Type , Degree >::DomainSize( width  );
	int domainH = FiniteElements1D< float , Type , Degree >::DomainSize( height );
	if( domainW<minRes ) fprintf( stderr , "[ERROR] Width is smaller than minimum-resolution: %d<%d\n" , domainW , minRes ) , exit( 0 );
	if( domainH<minRes ) fprintf( stderr , "[ERROR] Height is smaller than minimum-resolution: %d<%d\n" , domainH , minRes ) , exit( 0 );
	// Want domainW <= _paddedWidth = a * 2^b for the smallest a \in [ minRes , 2*minRes )
	// b is as small as possible the property that domainW <= (2*minRes-1)^2^b
	// -- by assumption (2*minRes-1)/2 < domainW
	// => the smallest possible value of b is 0.
#if 1
	long long aW = (2*minRes-1) , aH = (2*minRes-1);
	int bW=0 , bH=0;
	// First find the smallest value of b
	while( domainW>(aW<<bW) ) bW++ ; while( domainH>(aH<<bH) ) bH++;
	widthBlockSize = 1<<bW , heightBlockSize = 1<<bH;
	// Now find the smallest value of a
	while( domainW<=( (aW-1)<<bW ) ) aW-- ; while( domainH<=( (aH-1)<<bH ) ) aH--;
	long long _paddedWidth = (aW<<bW) , _paddedHeight = (aH<<bH);
#else
	long long _paddedWidth = minRes*2-1 , _paddedHeight = minRes*2-1;
	widthBlockSize = heightBlockSize = 1;
	while( domainW>_paddedWidth  ) _paddedWidth  *= 2 , widthBlockSize  *= 2;
	while( domainH>_paddedHeight ) _paddedHeight *= 2 , heightBlockSize *= 2;
	// Now we have b such that (2*minRes-1)*2^{b-1} < domainW <= (2*minRes-1)*2^b

	while( _paddedWidth - widthBlockSize>=domainW ) _paddedWidth  -= widthBlockSize;
	while( _paddedHeight-heightBlockSize>=domainH ) _paddedHeight -= heightBlockSize;
#endif
	paddedWidth  = (int)( FiniteElements1D<float,Type,Degree>::Dimension( (int)( _paddedWidth  ) ) );
	paddedHeight = (int)( FiniteElements1D<float,Type,Degree>::Dimension( (int)( _paddedHeight ) ) );
}
void SetPaddedSize( int width , int height , int minRes , int& paddedWidth , int& paddedHeight )
{
	int widthBlockSize, heightBlockSize;
	SetPaddedSize( width , height , minRes , paddedWidth , paddedHeight , widthBlockSize , heightBlockSize );
}
bool IsDownSamplable( int width , int height , int iters , int& newWidth , int& newHeight )
{
	if( (width&3) || (height&3) || (width<20) )	return false;
	newWidth  = width>>1;
	newHeight = height>>1;
	return true;
}
bool IsDownSamplable(int width,int height,int& newWidth,int& newHeight)
{
	if( (width&3) || (height&3) || (width<20) )	return false;
	newWidth  = width>>1;
	newHeight = height>>1;
	return true;
}
bool IsDownSamplable(int dCount,ClientSocket* cSockets,int sockCount)
{
	for(int i=0;i<sockCount;i++)
		if ( (cSockets[i].clientData.start>>dCount)&3 || ((cSockets[i].clientData.end-cSockets[i].clientData.start)>>dCount)<20)	return false;
	return true;
}

template<class Real>
void SetRectangularLaplacianMatrix( SparseMatrix< Real > & lap , int width , int height , double iWeight=0 )
{
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::FullDotProductStencil dotStencil1,d2DotStencil1,dotStencil2,d2DotStencil2;
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,dotStencil1,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,d2DotStencil1,1,1,false);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,dotStencil2,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,d2DotStencil2,1,1,false);
	SetRectangularLaplacianMatrix( dotStencil1 , d2DotStencil1 , dotStencil2 , d2DotStencil2 , lap , width , height , iWeight );
}
template<class Real>
void SetCylindricalLaplacianMatrix( SparseMatrix< Real > & lap , int width , int height , double iWeight=0 )
{
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::FullDotProductStencil dotStencil1,d2DotStencil1,dotStencil2,d2DotStencil2;
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,dotStencil1,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,d2DotStencil1,1,1,false);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,dotStencil2,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,d2DotStencil2,1,1,false);
	SetCylindricalLaplacianMatrix( dotStencil1 , d2DotStencil1 , dotStencil2 , d2DotStencil2 , lap , width , height , iWeight );
}
template<class Real>
void SetSphericalLaplacianMatrix(SparseMatrix<Real>& lap,int width,int height)
{
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::FullDotProductStencil dotStencil1,d2DotStencil1,dotStencil2,d2DotStencil2;
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,dotStencil1,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,d2DotStencil1,1,1,false);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,dotStencil2,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,d2DotStencil2,1,1,false);
	SetSphericalLaplacianMatrix( dotStencil1 , d2DotStencil1 , dotStencil2 , d2DotStencil2 , lap , width , height );
}
template<class Real>
void SetRectangularLaplacianMatrix( DotProductStencil& dotMajor , DotProductStencil& d2DotMajor , DotProductStencil& dotMinor , DotProductStencil& d2DotMinor ,
								    SparseMatrix< Real > & lap , int width , int height , double iWeight=0 )
{
	lap.Resize( width*height );

	int ii , jj;
	int xStart , xEnd , yStart , yEnd;
	for(int y=0;y<height;y++)
	{
		if		(y<Degree)			jj = y							, yStart = -y		, yEnd = Degree;
		else if	(y>height-1-Degree)	jj = 2*Degree+(y-(height-1))	, yStart = -Degree	, yEnd = height-1-y;
		else						jj = Degree						, yStart = -Degree	, yEnd = Degree;

		for(int x=0;x<width;x++)
		{
			if		(x<Degree)			ii = x						, xStart = -x		, xEnd = Degree;
			else if	(x>width-1-Degree)	ii = 2*Degree+(x-(width-1))	, xStart = -Degree	, xEnd = width-1-x;
			else						ii = Degree					, xStart = -Degree	, xEnd = Degree;

			int idx = x+y*width;
			lap.SetGroupSize( idx , (yEnd-yStart+1)*(xEnd-xStart+1) );
			int _i = 0;
			for(int yy=yStart;yy<=yEnd;yy++)
				for(int xx=xStart;xx<=xEnd;xx++)
				{
#if 1
					lap[idx][_i  ].N=(x+xx)+(y+yy)*width;
					lap[idx][_i++].Value=
						dotMajor.caseTable[ii].values[xx+Degree]*d2DotMinor.caseTable[jj].values[yy+Degree]+
						d2DotMajor.caseTable[ii].values[xx+Degree]*dotMinor.caseTable[jj].values[yy+Degree]+
						(dotMajor.caseTable[ii].values[xx+Degree]*dotMinor.caseTable[jj].values[yy+Degree])*iWeight;
#else
					lap.m_ppElements[idx][_i  ].N=(x+xx)+(y+yy)*width;
					lap.m_ppElements[idx][_i++].Value=
						dotMajor.caseTable[ii].values[xx+Degree]*d2DotMinor.caseTable[jj].values[yy+Degree]+
						d2DotMajor.caseTable[ii].values[xx+Degree]*dotMinor.caseTable[jj].values[yy+Degree]+
						(dotMajor.caseTable[ii].values[xx+Degree]*dotMinor.caseTable[jj].values[yy+Degree])*iWeight;
#endif
				}
		}
	}
}
template<class Real>
void SetCylindricalLaplacianMatrix( DotProductStencil& dotMajor , DotProductStencil& d2DotMajor , DotProductStencil& dotMinor , DotProductStencil& d2DotMinor ,
								    SparseMatrix< Real > & lap , int width , int height , double iWeight=0 )
{
	lap.Resize( width*height );

	int jj;
	int yStart , yEnd;
	for( int y=0 ; y<height ; y++ )
	{
		if		(y<Degree)			jj = y							, yStart = -y		, yEnd = Degree;
		else if	(y>height-1-Degree)	jj = 2*Degree+(y-(height-1))	, yStart = -Degree	, yEnd = height-1-y;
		else						jj = Degree						, yStart = -Degree	, yEnd = Degree;

		for( int x=0 ; x<width ; x++ )
		{
			int idx = x+y*width;
			lap.SetGroupSize( idx , (yEnd-yStart+1)*(2*Degree+1) );
			int _i = 0;
			for(int yy=yStart;yy<=yEnd;yy++)
				for( int i=-Degree ; i<=Degree ; i++ )
				{
					int xx=(x+i+width)%width;
#if 1
					lap[idx][_i  ].N=xx+(y+yy)*width;
					lap[idx][_i++].Value=
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[jj].values[yy+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[jj].values[yy+Degree]+
						(dotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[jj].values[yy+Degree])*iWeight;
#else
					lap.m_ppElements[idx][_i  ].N=(x+xx)+(y+yy)*width;
					lap.m_ppElements[idx][_i++].Value=
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[jj].values[yy+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[jj].values[yy+Degree]+
						(dotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[jj].values[yy+Degree])*iWeight;
#endif
				}
		}
	}
}
template<class Real>
void SetSphericalLaplacianMatrix(DotProductStencil& dotMajor,DotProductStencil& d2DotMajor,DotProductStencil& dotMinor,DotProductStencil& d2DotMinor,
								 SparseMatrix<Real>& lap,int width,int height)
{
	lap.Resize(width*height);

	for(int x=0;x<width;x++)
		for(int y=0;y<height;y++)
		{
			int idx1 = x + y*width;
			lap.SetGroupSize( idx1 , (2*Degree+1)*(2*Degree+1) );
			for(int i=-Degree;i<=Degree;i++)
				for(int j=-Degree;j<=Degree;j++)
				{
					int idx2 = (i+Degree)+(j+Degree)*(2*Degree+1);
					int xx=(x+i+width)%width;
					int yy= y+j;
					if(yy<0)
					{
						yy = -yy-1;
						if(xx<width/2)	xx = width/2-1-xx;
						else			xx = width/2+width-1-xx;
					}
					else if(yy>=height)
					{
						yy = height-1-(yy-height);
						xx = width-1-xx;
					}
#if 1
					lap[idx1][idx2].N = xx+yy*width;
					lap[idx1][idx2].Value = 
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[Degree].values[j+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree];
#else
					lap.m_ppElements[idx1][idx2].N = xx+yy*width;
					lap.m_ppElements[idx1][idx2].Value = 
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[Degree].values[j+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree];
#endif
				}
		}
}
template< class Real >
void UpSampleSphericalSolution( const Vector< Real >&low , int width , int height , Vector< Real >& high )
{
	int width2 , height2;
	FiniteElements1D< double , Type , Degree >::FullProlongationStencil wStencil , hStencil;
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( width  , wStencil , width2  );
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( height , hStencil , height2 );
	high.Resize( width2 * height2 );
	for( int i=0 ; i<width ; i++ )
	{
		int iStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( i );
		for( int j=0 ; j<height ; j++ )
		{
			int jStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( j );
			for( int x=0 ; x<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; x++ )
			{
				for( int y=0 ; y<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; y++ )
				{
					int xx = ( iStart + x + width2 ) % width2;
					int yy= jStart+y;
					if( yy<0 )
					{
						yy = -yy-1;
						if( xx<width2/2 )	xx = width2/2          - 1 - xx;
						else				xx = width2/2 + width2 - 1 - xx;
					}
					else if( yy>=height2 )
					{
						yy = height2-1-(yy-height2);
						xx = width2-1-xx;
					}
					high[ yy*width2 + xx ] += Real( low[j*width+i] * wStencil.caseTable[Degree].values[x] * hStencil.caseTable[Degree].values[y] );
				}
			}
		}
	}
}
template< class Real >
void UpSampleCylindricalSolution( const Vector< Real >&low , int width , int height , Vector< Real >& high )
{
	int width2 , height2;
	FiniteElements1D< double , Type , Degree >::FullProlongationStencil wStencil , hStencil;
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( width  , wStencil , width2 );
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( height , hStencil , height2 );
	high.Resize( width2 * height2 );
	for( int i=0 ; i<width ; i++ )
	{
		int iStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( i );
		for( int j=0 ; j<height ; j++ )
		{
			int jj;
			if( j<Degree )				jj = j;
			else if( j>=height-Degree)	jj = 2*Degree+(j-(height-1));
			else						jj = Degree;
			int jStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( j );
			for( int x=0 ; x<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; x++ )
			{
				int xx = ( iStart + x + width2 ) % width2;
				for( int y=0 ; y<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; y++ )
					if( jStart+y<0 || jStart+y>=height2 ) continue;
					else 
						high[ (jStart+y)*width2 + xx ] += Real( low[j*width+i] * wStencil.caseTable[Degree].values[x] * hStencil.caseTable[jj].values[y] );
			}
		}
	}
}
template< class Real >
void UpSampleRectangularSolution( const Vector< Real >&low , int width ,int height , Vector< Real >& high )
{
	int width2 , height2;
	FiniteElements1D< double , Type , Degree >::FullProlongationStencil wStencil , hStencil;
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( width  , wStencil , width2 );
	FiniteElements1D< double , Type , Degree >::ProlongationStencil( height , hStencil , height2 );
	high.Resize( width2 * height2 );
	for( int i=0 ; i<width ; i++ )
	{
		int iStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( i );
		int ii;
		if( i<Degree )				ii = i;
		else if( i>=width-Degree)	ii = 2*Degree+(i-(width-1));
		else						ii = Degree;
		for( int j=0 ; j<height ; j++ )
		{
			int jj;
			if( j<Degree )				jj = j;
			else if( j>=height-Degree)	jj = 2*Degree+(j-(height-1));
			else						jj = Degree;
			int jStart = FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Start( j );
			for( int x=0 ; x<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; x++ )
				if( iStart+x<0 || iStart+x>=width2 ) continue;
				else
					for( int y=0 ; y<FiniteElements1D< double , Type , Degree >::FullProlongationStencil::ProlongationStencil::Size ; y++ )
						if( jStart+y<0 || jStart+y>=height2 ) continue;
						else 
							high[ (jStart+y)*width2 + (iStart+x) ] += Real( low[j*width+i] * wStencil.caseTable[ii].values[x] * hStencil.caseTable[jj].values[y] );
		}
	}
}

template<class Real>
void SetSphericalLaplacianMatrix(SparseMatrix<Real>& lap,int width,int height,double iWeight)
{
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::FullDotProductStencil dotStencil1,d2DotStencil1,dotStencil2,d2DotStencil2;
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,dotStencil1,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(width,d2DotStencil1,1,1,false);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,dotStencil2,0,0);
	FiniteElements1D<double,Type,Degree>::DotProduct<Type,Degree>::DotProductStencil(height,d2DotStencil2,1,1,false);
	SetSphericalLaplacianMatrix( dotStencil1 , d2DotStencil1 , dotStencil2 , d2DotStencil2 , lap , width , height , iWeight );
}
template<class Real>
void SetSphericalLaplacianMatrix(DotProductStencil& dotMajor,DotProductStencil& d2DotMajor,DotProductStencil& dotMinor,DotProductStencil& d2DotMinor,
								 SparseMatrix<Real>& lap,int width,int height,double iWeight)
{
	lap.Resize(width*height);

	for(int x=0;x<width;x++)
		for(int y=0;y<height;y++)
		{
			int idx1 = x + y*width;
			lap.SetGroupSize( idx1 , (2*Degree+1)*(2*Degree+1) );
			for(int i=-Degree;i<=Degree;i++)
				for(int j=-Degree;j<=Degree;j++)
				{
					int idx2 = (i+Degree)+(j+Degree)*(2*Degree+1);
					int xx=(x+i+width)%width;
					int yy= y+j;
					if(yy<0)
					{
						yy = -yy-1;
						if(xx<width/2)	xx = width/2-1-xx;
						else			xx = width/2+width-1-xx;
					}
					else if(yy>=height)
					{
						yy = height-1-(yy-height);
						xx = width-1-xx;
					}

#if 1
					lap[idx1][idx2].N = xx+yy*width;
					lap[idx1][idx2].Value = 
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[Degree].values[j+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree]+
						(dotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree])*iWeight;
#else
					lap.m_ppElements[idx1][idx2].N = xx+yy*width;
					lap.m_ppElements[idx1][idx2].Value = 
						dotMajor.caseTable[Degree].values[i+Degree]*d2DotMinor.caseTable[Degree].values[j+Degree]+
						d2DotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree]+
						(dotMajor.caseTable[Degree].values[i+Degree]*dotMinor.caseTable[Degree].values[j+Degree])*iWeight;
#endif
				}
		}
}
///////////////////////////
// SphericalSynchronizer //
///////////////////////////
SphericalSynchronizer::SphericalSynchronizer(void)
{
	_cCount = 0;
	_clientSockets = NULL;
}
void SphericalSynchronizer::init( int width , const int* widths , int cCount , int rPasses , int pPasses , int vCycles )
{
	_width=width;
	_cCount=cCount;
	_rPasses=rPasses;
	_pPasses=pPasses;
	_vCycles=vCycles;

	if( _clientSockets ) delete[] _clientSockets;

	_clientSockets = new ClientSocket[cCount];
	if(!_clientSockets)
	{
		fprintf(stderr,"Failed to allocate SphericalSynchronizer::_clientSockets\n");
		exit(0);
	}
	for( int i=0 ; i<cCount ; i++ )
	{
		_clientSockets[i].listener = GetListenSocket(_clientSockets[i].port);
		if(_clientSockets[i].listener == _INVALID_ACCEPTOR_SOCKET_ ) fprintf(stderr,"Failed to set server in SphericalSynchronizer\n") , exit(0);
		if( i ) _clientSockets[i].clientData.start = _clientSockets[i-1].clientData.end;
		else    _clientSockets[i].clientData.start = 0;
		_clientSockets[i].clientData.end = _clientSockets[i].clientData.start+widths[i];
	}
}
SphericalSynchronizer::~SphericalSynchronizer(void)
{
	delete[] _clientSockets;
	_clientSockets=NULL;
}
int SphericalSynchronizer::clients(void)		const	{ return _cCount; }
int SphericalSynchronizer::port(int cIndex)	const	{ return _clientSockets[cIndex].port; }
template< class DataType >
void SphericalSynchronizer::Run( void )
{
	for(int i=0;i<_cCount;i++)
	{
		_clientSockets[i].socket = AcceptSocket( _clientSockets[i].listener );
		if ( _clientSockets[i].socket == _INVALID_SOCKET_ ) exit(0);
	}

	Pointer( DataType ) buffer1;
	Pointer( DataType ) buffer2;
	buffer1 = AllocArray< DataType >( _width , "SphericalSynchronizer::Run (buffer1)" );
	buffer2 = AllocArray< DataType >( _width , "SphericalSynchronizer::Run (buffer2)" );
	for ( int v = 0 ; v < _vCycles ; v++ )
	{
		//////////////////////////////////////////
		// First synchronize on the restriction //
		//////////////////////////////////////////

		// Synchronize the head of the stream
		for(int k=0;k<_rPasses;k++)
		{
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				ReceiveOnSocket( _clientSockets[i].socket , buffer1 , sizeof(DataType)*w , "Failed to get restriction head from client" );
				for( int x=0 ; x<w ; x++ ) buffer2[x+_clientSockets[i].clientData.start] = buffer1[x];
			}
			for(int x=0;x<_width/2;x++)
			{
				int xx = _width/2-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int x=_width/2;x<_width;x++)
			{
				int xx = _width/2+_width-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				for( int x=0 ; x<w ; x++ ) buffer2[x] = buffer1[x+_clientSockets[i].clientData.start];
				SendOnSocket( _clientSockets[i].socket , buffer2 , sizeof(DataType)*w , "Failed to send restriction head to client" ) ;
			}
		}
		// Synchronize the tail of the stream
		for(int k=0;k<_rPasses;k++)
		{
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				ReceiveOnSocket( _clientSockets[i].socket , buffer1 , sizeof(DataType)*w , "Failed to get restriction tail from client" );
				for( int x=0 ; x<w ; x++ ) buffer2[x+_clientSockets[i].clientData.start] = buffer1[x];
			}
			for( int x=0 ; x<_width ; x++ )
			{
				int xx = _width-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				for( int x=0 ; x<w ; x++ ) buffer2[x] = buffer1[x+_clientSockets[i].clientData.start];
				SendOnSocket( _clientSockets[i].socket , buffer2 , sizeof(DataType)*w , "Failed to send restriction tail to client" );
			}
		}

		/////////////////////////////////////////
		// Now synchronize on the prolongation //
		/////////////////////////////////////////

		// Synchronize the head of the stream
		for(int k=0;k<_pPasses;k++)
		{
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				ReceiveOnSocket( _clientSockets[i].socket , buffer1 , sizeof(DataType)*w , "Failed to get prolongation head from client" );
				for( int x=0 ; x<w ; x++ ) buffer2[x+_clientSockets[i].clientData.start] = buffer1[x];
			}
			for(int x=0;x<_width/2;x++)
			{
				int xx = _width/2-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int x=_width/2;x<_width;x++)
			{
				int xx = _width/2+_width-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				for( int x=0 ; x<w ; x++ ) buffer2[x] = buffer1[x+_clientSockets[i].clientData.start];
				SendOnSocket( _clientSockets[i].socket , buffer2 , sizeof(DataType)*w , "Failed to send prolongation head to client" );
			}
		}
		// Synchronize the tail of the stream
		for(int k=0;k<_pPasses;k++)
		{
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				ReceiveOnSocket( _clientSockets[i].socket , buffer1 , sizeof(DataType)*w , "Failed to get prolongation tail from client" );
				for( int x=0 ; x<w ; x++ ) buffer2[x+_clientSockets[i].clientData.start] = buffer1[x];
			}
			for(int x=0;x<_width;x++)
			{
				int xx = _width-1-x;
				buffer1[xx] = buffer2[x];
			}
			for(int i=0;i<_cCount;i++)
			{
				int w = _clientSockets[i].clientData.end-_clientSockets[i].clientData.start;
				for( int x=0 ; x<w ; x++ ) buffer2[x] = buffer1[x+_clientSockets[i].clientData.start];
				SendOnSocket( _clientSockets[i].socket , buffer2 , sizeof(DataType)*w , "Failed to send prolongation tail to client" );
			}
		}
	}
	FreeArray( buffer1 );
	FreeArray( buffer2 );
}
template< class DataType >
THREAD_FUNCTION_OUT SphericalSynchronizer::RunThread( THREAD_FUNCTION_IN vparams )
{
	SphericalSynchronizer* synchronizer = (SphericalSynchronizer*)vparams;
	synchronizer->Run< DataType >();
	return 0;
}
/////////////////////////////
// SocketedMultigridServer //
/////////////////////////////
template< int PixelChannels , int LabelChannels , class SyncType >
SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::SocketedMultigridServer(void)
{
	_synchronizerHandles = NULL;
	_icSynchronizers = NULL;
	_icSyncSockets = NULL;
}
template< int PixelChannels , int LabelChannels , class SyncType >
SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::~SocketedMultigridServer(void)
{
	if( _globalData.periodicType==SPHERICAL_PERIODIC )
	{
		WaitOnThreads( _synchronizerHandles , _oocLevels+1+_icLevels , 10000 , "SocketedMultigridServer::~SocketedMultigridServer" );
		for( int i=0 ; i<_oocLevels+1+_icLevels ; i++ ) CloseThreadHandle( _synchronizerHandles[i] );
		delete[] _synchronizerHandles;
		_synchronizerHandles = NULL;
		if( _icSyncSockets )
		{
			for( int i=0;i<_icLevels ; i++ ) CloseSocket( _icSyncSockets[i] );
			delete[] _icSyncSockets;
		}
		_icSyncSockets = NULL;
	}
	if( _icSynchronizers )	delete[] _icSynchronizers;
	for( size_t i=0 ; i<_oocSynchronizers.size() ; i++ ) if( _oocSynchronizers[i] ) delete[] _oocSynchronizers[i];
	_oocSynchronizers.clear();
	_icSynchronizers = NULL;
}
template< int PixelChannels , int LabelChannels , class SyncType >
template< class LabelType >
bool SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::SetUp( char* address , int port , int clientCount , int iters , int inCoreRes , int minMGRes , int vCycles, int minBandSize ,
											   int tileWidth , int tileHeight , const char* tileExt , 
											   bool gammaCorrection , int quality , int lanes , bool verbose , int periodicType , double iWeight , double gScale ,
											   bool removeAverageGradient , int unknownType , bool showProgress , bool noCG , bool shortSync , int* clipDimensions )
{
	// Create a listening Socket for connecting to server
	Socket listenSocket = _INVALID_SOCKET_;
	listenSocket = GetListenSocket( port );
	if( listenSocket == _INVALID_SOCKET_ ) return false;
	if( address )
	{
		printf( "Server Address: %s:%d\n", address , port ) , fflush( stdout );
		return SetUp( address , listenSocket , clientCount , iters , inCoreRes , minMGRes , vCycles , minBandSize , tileWidth , tileHeight , tileExt , gammaCorrection , quality , lanes , verbose , periodicType , iWeight , gScale , removeAverageGradient , unknownType , showProgress , noCG , shortSync , clipDimensions );
	}
	else
	{
		char hostAddress[512];
		GetHostAddress( hostAddress );
		printf( "Host Address: %s:%d\n", hostAddress , port ) , fflush( stdout );
		return SetUp( hostAddress , listenSocket , clientCount , iters , inCoreRes , minMGRes , vCycles , minBandSize , tileWidth , tileHeight , tileExt , gammaCorrection , quality , lanes , verbose , periodicType , iWeight , gScale , removeAverageGradient , unknownType , showProgress , noCG , shortSync , clipDimensions );
	}
}
template< int PixelChannels , int LabelChannels , class SyncType >
template< class LabelType >
bool SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::SetUp( char* address , AcceptorSocket listenSocket , int clientCount , int iters , int inCoreRes , int minMGRes , int vCycles , int minBandSize ,
														   int tileWidth , int tileHeight , const char* tileExt , 
											   bool gammaCorrection , int quality , int lanes , bool verbose , int periodicType , double iWeight , double gScale ,
											   bool removeAverageGradient , int unknownType , bool showProgress , bool noCG , bool shortSync , int* clipDimensions )
{
	ClientSocket* clientSockets = new ClientSocket[clientCount];
	_noCG		= noCG;
	_inCoreRes	= inCoreRes;
	_minMGRes	= minMGRes;

	_globalData.periodicType= periodicType;
	_globalData.verbose		= verbose;
	_globalData.iters		= iters;
	_globalData.vCycles		= vCycles;
	_globalData.iWeight		= iWeight;
	_globalData.gScale		= gScale;
	_globalData.tileWidth	= tileWidth;
	_globalData.tileHeight	= tileHeight;
	_globalData.quality		= quality;
	_globalData.lanes		= lanes;
	_globalData.unknownType	= unknownType;
	_globalData.showProgress= showProgress;
	_globalData.shortSync	= shortSync;
	if( tileExt ) strcpy( _globalData.tileExt , tileExt );
	else          _globalData.tileExt[0] = 0;

	// Establish a connection to each of the processes
	for( int i=0 ; i<clientCount ; i++ )
	{
		clientSockets[i].socket = AcceptSocket( listenSocket );
		printf( "Connected to thread: %d /%d     \r" , i+1 , clientCount );
	}
	printf( "\n" ) , fflush( stdout );

	// Get the client info
	int width = 0 , height;
	for( int i=0 ; i<clientCount ; i++ )
	{
		ReceiveOnSocket( clientSockets[i].socket , GetPointer( clientSockets[i].clientData ) , sizeof( ClientData ) , "Failed to get client data from client" );
		if( !i ) height = clientSockets[i].clientData.height;
		else if( height!=clientSockets[i].clientData.height ) fprintf( stderr , "[ERROR] Band heights differ: %d != %d\n" , height , clientSockets[i].clientData.height ) , exit(0);
		width += clientSockets[i].clientData.width;
	}
	{
		std::vector< std::pair< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > > > averageMap;
		// Receive average gradient values
		std::map< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > > _averageMap;

		for( int i=0 ; i<clientCount ; i++ )
		{
			int count;
			ReceiveOnSocket( clientSockets[i].socket , GetPointer( count ) , sizeof( count ) , "Failed to get average gradient count from client" );
			if( count )
			{
				std::vector< std::pair< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > > > __averageMap( count );
				long long size = sizeof( std::pair< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > > ) * count;
				ReceiveOnSocket( clientSockets[i].socket , GetPointer( __averageMap[0] , count ) , (int)(size) ,  "Failed to get average gradients from client" );
				for( int j=0 ; j<__averageMap.size() ; j++ )
				{
					typename std::map< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > >::iterator iter = _averageMap.find( __averageMap[j].first );
					if( iter==_averageMap.end() ) _averageMap[ __averageMap[j].first ] = __averageMap[j].second;
					else iter->second += __averageMap[j].second;
				}
			}
		}
		for( typename std::map< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > >::iterator iter=_averageMap.begin() ; iter!=_averageMap.end() ; iter++ )
			averageMap.push_back( std::pair< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > >( iter->first , iter->second ) );

		for( int i=0 ; i<averageMap.size() ; i++ ) averageMap[i].second.average();

		// Send back average gradient values
		for( int i=0 ; i<clientCount ; i++ )
		{
			int count = (int)( averageMap.size() );
			SendOnSocket( clientSockets[i].socket , GetPointer( count ) , int (sizeof( count ) ) , "Failed to send average gradient count to client" );
			if( count )
			{
				long long size = sizeof( std::pair< LabelData< LabelType , LabelChannels > , GradientAverage< PixelChannels > > ) * count;
				SendOnSocket( clientSockets[i].socket , GetPointer( averageMap[0] , count ) , (int)(size) ,  "Failed to send average gradients to client" );
			}
		}
	}

	// Sort the clients by index order
	qsort( clientSockets , clientCount , sizeof(ClientSocket) , ClientSocket::Sort );

	// Store the clipping width and then compute the padded dimensions
	// Only compute the padded dimensions if we don't already know them.
	if( clipDimensions )
	{
		_globalData.cWidth	= clipDimensions[0];
		_globalData.cHeight	= clipDimensions[1];
		_globalData.width	=  width;
		_globalData.height	= height;
	}
	else
	{
		_globalData.cWidth  = width;
		_globalData.cHeight = height;

#if 0
#pragma message( "[WARNING] For the life of me, why do I have to make this change" )
		_globalData.width = width  , _globalData.height = height;
#else
		if     ( periodicType==         NO_PERIODIC ) _globalData.width = width+1, _globalData.height = height+1;
		else if( periodicType==CYLINDRICAL_PERIODIC ) _globalData.width = width  , _globalData.height = height+1;
		else                                          _globalData.width = width  , _globalData.height = height  ;
#endif

		int w = _globalData.width , h = _globalData.height , wBS, hBS;
		SetPaddedSize( w , h , _minMGRes, _globalData.width , _globalData.height , wBS , hBS );
		if( periodicType==CYLINDRICAL_PERIODIC && _globalData.width!=w ) fprintf( stderr , "Cylindrical width was not preserved: (%d , %d) -> (%d , %d)\t(%d , %d)\n" , w , h , _globalData.width , _globalData.height , wBS , hBS ) , exit( 0 );
		if( periodicType==  SPHERICAL_PERIODIC && ( _globalData.width!=w || _globalData.height!=h ) ) fprintf( stderr , "Spherical dimensions were not preserved: (%d , %d) -> (%d , %d)\t(%d , %d)\n" , w , h , _globalData.width , _globalData.height , wBS , hBS ) , exit( 0 );
	}
	// Send the general solver information
	for( int i=0 ; i<clientCount ; i++ ) SendOnSocket( clientSockets[i].socket , GetPointer( _globalData ) , sizeof( _globalData ) , "Failed to send global data to client" );

	// Set the local information for the bands
	int start=0;
	for( int i=0 ; i<clientCount ; i++ )
	{
		int w = clientSockets[i].clientData.width;
		clientSockets[i].clientData.start = start;
		clientSockets[i].clientData.end = clientSockets[i].clientData.cEnd = start+w;
		if( i==clientCount-1 ) clientSockets[i].clientData.end  = _globalData.width;
		start = clientSockets[i].clientData.end;
		SendOnSocket( clientSockets[i].socket , GetPointer( clientSockets[i].clientData ) , sizeof( ClientData ) , "Failed to send client data to client" );
	}

	// Figure out how many different out-of-core levels are needed
	_oocLevels = 1;
	while( (long long)( width>>_oocLevels ) * (long long)( height>>_oocLevels ) >= _inCoreRes*_inCoreRes ) _oocLevels++;
	if( !ProcessPartition<int>::IsValidBandSize( width>>_oocLevels , height>>_oocLevels , iters , minBandSize ) )
	{
		fprintf( stderr , "In-core resolution is too small: %d x %d\n" , width>>_oocLevels , height>>_oocLevels );
		return false;
	}
	// WARNING!!! It's possible that _oocLevels = 0 if the image fits in-core

	ProcessPartition< int > rootPartition;
	rootPartition.resize( clientCount );
	for( int i=0 ; i<clientCount ; i++ )
	{
		rootPartition[i].width = clientSockets[i].clientData.end - clientSockets[i].clientData.start;
		if( !ProcessPartition< int >::IsValidBandSize( rootPartition[i].width , _globalData.height , _globalData.iters , minBandSize ) )
		{
			fprintf( stderr , "Client[%d] has invalid dimension: %d x %d\n" , i , rootPartition[i].width , _globalData.height );
			return false;
		}
	}
	// WARNING!!! The "repeat" flag should be turned off if processors don't double-solve the same level

	if( !processPartition.Initialize( rootPartition , _globalData.height , _globalData.iters , _oocLevels-1 , false , minBandSize ) )
	{
		fprintf( stderr , "Failed to create process partition\n" );
		return false;
	}
	for( int i=0 ; i<clientCount ; i++ )
	{
		int idx = i;
		int tCount = 1;
		for( int j=1 ; j<processPartition.size() && !(idx&1) && ( (idx>>1)<processPartition[j].size() ); j++ ) tCount++ , idx>>=1;
		SendOnSocket( clientSockets[i].socket , GetPointer( tCount ) , sizeof( tCount ) , "failed to send tCount to client" );
		Pointer( ProcessData ) pd = AllocArray< ProcessData >( tCount , "SocketedMultigridServer::SetUp (pd)" );
		pd[0].index = 0 ;
		pd[0].maxIndex = processPartition.size()-1;
		pd[0].offset = i;
		pd[0].maxOffset = clientCount-1;
		pd[0].startDepth = processPartition[0].startDepth;
		pd[0].endDepth = processPartition[0].endDepth;
		processPartition[0].setBounds( i , pd[0].start , pd[0].stop );
		pd[0].width = _globalData.width;
		pd[0].height = _globalData.height;
		pd[0].children = (int)( processPartition[0][i].children.size() );

		idx = i;
		for( int j=1 ; j<tCount ; j++ )
		{
			idx >>= 1;
			pd[j].index  = j;
			pd[j].maxIndex = processPartition.size()-1;
			pd[j].offset = idx;
			pd[j].maxOffset = processPartition[j].size()-1;
			pd[j].startDepth = processPartition[j].startDepth;
			pd[j].endDepth = processPartition[j].endDepth;
			processPartition[j].setBounds( idx , pd[j].start , pd[j].stop );
			pd[j].width = _globalData.width >> processPartition[j].startDepth;
			pd[j].height = _globalData.height >> processPartition[j].startDepth;
			pd[j].children = (int)( processPartition[j][idx].children.size() );
		}
		SendOnSocket( clientSockets[i].socket , pd , sizeof(ProcessData) * tCount , "Failed to send process data to client" );
		FreeArray( pd );
	}
	delete[] clientSockets;

	// At this point we transition our communications so that we are talking to the threads, not processes
	for( int i=0 ; i<processPartition.size() ; i++ )
		for( int j=0 ; j<processPartition[i].size() ; j++ )
		{
			// Establish a connection to the process
			int index , offset;
			Socket sock = AcceptSocket( listenSocket );
			ReceiveOnSocket( sock , GetPointer( index  ) , sizeof( index  ) , "Failed to get client index" );
			ReceiveOnSocket( sock , GetPointer( offset ) , sizeof( offset ) , "Failed to get client offset" );
			processPartition[index][offset].data.sock = sock;
			if( index )
			{
				ReceiveOnSocket ( sock , GetPointer( processPartition[index][offset].data.childAddressX ) , sizeof( IPData ) , "Failed to get child address X" );
				ReceiveOnSocket ( sock , GetPointer( processPartition[index][offset].data.childAddressB ) , sizeof( IPData ) , "Failed to get child address B" );
			}
			if( periodicType!=NO_PERIODIC || offset )
			{
				ReceiveOnSocket ( sock , GetPointer( processPartition[index][offset].data.leftAddress ) , sizeof( IPData ) , "Failed to get left address from client" );
			}
		}

	// Create a listening socket to establish the connection between the in-core and out-of-core solver
	int listenerCount = 2;
	ThreadHandle listenerHandles[3];

	_inCoreConnectInfo.init( processPartition[processPartition.size()-1].size() );

	IPData inCoreAddressX , inCoreAddressB , inCoreAddressLeft;

	{
		IOServer::SystemLock lock;
		GetHostEndpointAddress( &inCoreAddressX.address );
		GetHostEndpointAddress( &inCoreAddressB.address );
		GetHostEndpointAddress( &inCoreAddressLeft.address );
	}
	_inCoreConnectInfo.listenChildX = GetListenSocket( inCoreAddressX.port );
	_inCoreConnectInfo.listenChildB = GetListenSocket( inCoreAddressB.port );
	if ( _inCoreConnectInfo.listenChildX == _INVALID_ACCEPTOR_SOCKET_ ) return false;
	if ( _inCoreConnectInfo.listenChildB == _INVALID_ACCEPTOR_SOCKET_ ) return false;
	listenerHandles[0] = SpawnListenerSocket( _inCoreConnectInfo.listenChildX , _inCoreConnectInfo.childX , processPartition[processPartition.size()-1].size() );
	listenerHandles[1] = SpawnListenerSocket( _inCoreConnectInfo.listenChildB , _inCoreConnectInfo.childB , processPartition[processPartition.size()-1].size() );
	if( !TestThreadHandle( listenerHandles[0] ) ) return false;
	if( !TestThreadHandle( listenerHandles[1] ) ) return false;

	Socket leftSocket = _INVALID_SOCKET_;
	if( periodicType!=NO_PERIODIC )
	{
		_inCoreConnectInfo.listenLeft = GetListenSocket( inCoreAddressLeft.port );
		if ( _inCoreConnectInfo.listenLeft == _INVALID_ACCEPTOR_SOCKET_ ) return false;
		listenerCount++;
		listenerHandles[2] = SpawnListenerSocket( _inCoreConnectInfo.listenLeft , &leftSocket , 1 );
		if( !TestThreadHandle( listenerHandles[2] ) ) return false;
	}

	// Now send everybody information about who to connect to for their data
	for( int i=0 ; i<processPartition.size() ; i++ )
		for( int j=0 ; j<processPartition[i].size() ; j++ )
		{
			if( i )
				for( size_t k=0 ; k<processPartition[i][j].children.size() ; k++ )
				{
					Socket& sock = processPartition[i-1][processPartition[i][j].children[k]].data.sock;
					SendOnSocket ( sock , GetPointer( processPartition[i][j].data.childAddressX ) , sizeof( IPData ) , "Failed to send child address X to client" );
					SendOnSocket ( sock , GetPointer( processPartition[i][j].data.childAddressB ) , sizeof( IPData ) , "Failed to send child address B to client" );
				}
		}
	{
		int i = processPartition.size()-1;
		for( int j=0; j<processPartition[i].size() ; j++ )
		{
			Socket& sock = processPartition[i][j].data.sock;
			SendOnSocket ( sock , GetPointer( inCoreAddressX ) , sizeof( IPData ) , "Failed to send address X to client" );
			SendOnSocket ( sock , GetPointer( inCoreAddressB ) , sizeof( IPData ) , "Failed to send address B to client" );
		}
	}
	for( int i=0 ; i<processPartition.size() ; i++ )
		for( int j=0 ; j<processPartition[i].size() ; j++ )
			if( periodicType!=NO_PERIODIC || j )
			{
				Socket& sock = processPartition[i][(j-1+processPartition[i].size())%processPartition[i].size()].data.sock;
				SendOnSocket ( sock , GetPointer( processPartition[i][j].data.leftAddress ) , sizeof( IPData ), "Failed to send left address to client" );
			}
	if( periodicType!=NO_PERIODIC )
		_inCoreConnectInfo.right = new SocketStream( inCoreAddressLeft.address , inCoreAddressLeft.port , 5 , false );
	WaitOnThreads( listenerHandles , listenerCount , 10000 , "SocketedMultigridServer::SetUp" );
	if( periodicType!=NO_PERIODIC ) _inCoreConnectInfo.left = new SocketStream( leftSocket );

	// Figure out how many different in-core resolutions are needed
	_icLevels = 1;
	int w = _globalData.width >> _oocLevels , h = _globalData.height >> _oocLevels;
	while( ProcessPartition<int>::IsValidBandSize( w>>_icLevels , h>>_icLevels , _globalData.iters , 0 ) && (w>>_icLevels)>=_minMGRes && (h>>_icLevels)>=_minMGRes )
		_icLevels++;

	// Now for the spherical synchronizers
	if( _globalData.periodicType==SPHERICAL_PERIODIC )
	{
		_icSynchronizers = new SphericalSynchronizer[ _icLevels ];
		if( !_icSynchronizers )
		{
			fprintf(stderr,"Failed to allocate in-core synchronizers\n");
			return false;
		}
		_oocSynchronizers.resize( processPartition.size() );
		for( int i=0 ; i<processPartition.size() ; i++ )
		{
			if( !i ) _oocSynchronizers[i] = new SphericalSynchronizer[ (processPartition[i].endDepth-processPartition[i].startDepth+1)+1 ];
			else	 _oocSynchronizers[i] = new SphericalSynchronizer[ (processPartition[i].endDepth-processPartition[i].startDepth+1)   ];
		}

		int *oocWidths , icWidths;
		oocWidths = new int[clientCount];

		for( int i=0 ; i<processPartition.size() ; i++ )
		{
			if(!i)
			{
				for( int j=0 ; j<processPartition[i].size() ; j++ ) oocWidths[j] = processPartition[i][j].width;
				imageSynchronizer.init( _globalData.width , oocWidths , clientCount , Degree , 0 , 1 );
			}
			for( int d=processPartition[i].startDepth ; d<=processPartition[i].endDepth ; d++ )
			{
				int depth = d-processPartition[i].startDepth;
				for( int j=0 ; j<processPartition[i].size() ; j++ ) oocWidths[j] = processPartition[i][j].width >> depth;
				if( _globalData.verbose )
					_oocSynchronizers[i][depth].init( _globalData.width>>d , oocWidths , processPartition[i].size() , (_globalData.iters*Degree+Degree)+Degree , (_globalData.iters*Degree+Degree)+Degree , _globalData.vCycles );
				else
					_oocSynchronizers[i][depth].init( _globalData.width>>d , oocWidths , processPartition[i].size() , (_globalData.iters*Degree+Degree)+Degree , (_globalData.iters*Degree+Degree)+     0 , _globalData.vCycles );
			}
		}
		delete[] oocWidths;
		for( int d=0 ; d<_icLevels ; d++ )
		{
			int depth = d+_oocLevels;
			icWidths = _globalData.width >> depth;
			if( _globalData.verbose )
				_icSynchronizers[d].init( icWidths , &icWidths , 1 , (_globalData.iters*Degree+Degree)+Degree , (_globalData.iters*Degree+Degree)+Degree , _globalData.vCycles );
			else
				_icSynchronizers[d].init( icWidths , &icWidths , 1 , (_globalData.iters*Degree+Degree)+Degree , (_globalData.iters*Degree+Degree)+     0 , _globalData.vCycles );
		}
		_synchronizerHandles = new ThreadHandle[ _oocLevels + 1 + _icLevels ];

		int idx = 0;
		{
			// [WARNING] The spherical synchronizer should only be enabled depending on whether or not there are
			// labels, not the interpolation weight
			if( _globalData.iWeight )
				_synchronizerHandles[idx] = RunThread( SphericalSynchronizer::RunThread< SyncType > , &imageSynchronizer );
			else
				_synchronizerHandles[idx] = RunThread( SphericalSynchronizer::RunThread< ImageData< PixelChannels , LabelChannels , SyncType , uint16_t > > , &imageSynchronizer );
			if( !TestThreadHandle( _synchronizerHandles[idx++] ) ) fprintf( stderr , "Failed to create sync thread\n" ) , exit(0);
		}
		for( int i=0 ; i<processPartition.size() ; i++ )
		{
			int depths = processPartition[i].endDepth - processPartition[i].startDepth + 1;
			for( int d=0 ; d<depths ; d++ )
			{
				_synchronizerHandles[idx] = RunThread( SphericalSynchronizer::RunThread< SyncType > , &_oocSynchronizers[i][d] );
				if( !TestThreadHandle( _synchronizerHandles[idx++] ) ) fprintf( stderr , "Failed to create sync thread\n" ) , exit(0);
			}
		}
		for(int i=0;i<_icLevels;i++)
		{
			_synchronizerHandles[idx] = RunThread( SphericalSynchronizer::RunThread< float > , &_icSynchronizers[i] );
			if( !TestThreadHandle( _synchronizerHandles[idx++] ) ) fprintf( stderr , "Failed to create sync thread\n" ) , exit(0);
		}

		int* icPorts = new int[_icLevels];
		for( int d=0 ; d<_icLevels ; d++ )	icPorts[d] = _icSynchronizers[d].port(0);

		for( int i=0 ; i<processPartition.size() ; i++ )
		{
			int depths = processPartition[i].endDepth - processPartition[i].startDepth + 1;
			Pointer( int ) oocPorts = AllocArray< int >( depths , "SocketedMultigridServer::SetUp (oocPorts)" );
			for( int j=0 ; j< processPartition[i].size() ; j++ )
			{
				if( !i )
				{
					int imagePort = imageSynchronizer.port(j);
					SendOnSocket( processPartition[i][j].data.sock , GetPointer( imagePort ) , sizeof( int ) , "Failed to send image port to client" );
				}
				for( int d=0 ; d<depths ; d++ ) oocPorts[d] = _oocSynchronizers[i][d].port(j);
				SendOnSocket( processPartition[i][j].data.sock , oocPorts , sizeof(int)*depths , "Failed to send out of core ports to client" );
			}
			FreeArray( oocPorts );
		}
		_icSyncSockets = new Socket[_icLevels];
		for( int i=0 ; i<_icLevels ; i++ )
		{
			_icSyncSockets[i] = GetConnectSocket( address , icPorts[i] );
			if( _icSyncSockets[i] == _INVALID_SOCKET_ )	return false;
		}
		delete[] icPorts;
	}
	CloseAcceptorSocket( listenSocket );
	return true;
}
template< int PixelChannels , int LabelChannels , class SyncType >
void SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::Run( void )
{
#if MISHA_DENORMAL_CONTROL
//	unsigned int control_word_x87  , control_word_sse2;
//	__control87_2( _DN_FLUSH , _MCW_DN , &control_word_x87 , &control_word_sse2 );
//	_controlfp( _DN_FLUSH , _MCW_DN );
	_MM_SET_FLUSH_ZERO_MODE( _MM_FLUSH_ZERO_ON );
#endif // MISHA_DENORMAL_CONTROL
	StreamingGrid *inB , *outX;
	printf( "Running server\n" ) , fflush( stdout );
	Pointer( SolverInfo< PixelChannels > )  solverInfo = AllocArray< SolverInfo< PixelChannels > >( _oocLevels , "SocketedMultigridServer::Run (solverInfo)" );
	Pointer( SolverInfo< PixelChannels > ) tSolverInfo = AllocArray< SolverInfo< PixelChannels > >( _oocLevels , "SocketedMultigridServer::Run (tSolverInfo)" );
	double zeroAverage[PixelChannels];
	for(int c=0;c<PixelChannels;c++)	zeroAverage[c]=0;
	Vector<float> globalB,globalX;
	int gWidth  = _globalData.width  >> _oocLevels;
	int gHeight = _globalData.height >> _oocLevels;

	globalB.Resize( gWidth * gHeight * PixelChannels );
	globalX.Resize( gWidth * gHeight * PixelChannels );

//	if(_globalData.verbose)	MultigridSolver<float,ZERO_DERIVATIVE,2,Channels>::verbose = MultigridSolver<float,ZERO_DERIVATIVE,2,Channels>::FULL_VERBOSE;

	memset(_average, 0 , sizeof(_average) );
	inB  = GetMultiSocketBackedGrid< PixelChannels , SyncType > ( _inCoreConnectInfo.childB , _inCoreConnectInfo.children , gHeight    , false , true );
	outX = GetMultiSocketBackedGrid< PixelChannels , SyncType > ( _inCoreConnectInfo.childX , _inCoreConnectInfo.children , gHeight<<1 , false , true );

	for( int ii=0 ; ii<_globalData.vCycles ; ii++ )
	{
		double t;

		t=Time();
		memset( solverInfo , 0 , sizeof( SolverInfo< PixelChannels > ) * _oocLevels );
		inB->reset( true , 1 );

		for( int j=0 ; j<gHeight ; j++ )
		{
			Pointer( SyncType ) row = ( Pointer( SyncType ) )(*inB)[j];
//			memcpy( &globalB[gWidth*j*Channels] , row , sizeof(float)*Channels*gWidth );
			for( int i=0 ; i<PixelChannels*gWidth ; i++ ) globalB[gWidth*j*PixelChannels + i ] = float( row[i] );
			inB->advance();
		}
		for( int i=0 ; i<processPartition.size() ; i++ )
			for( int j=0 ; j<processPartition[i].size() ; j++ )
			{
				ReceiveOnSocket( processPartition[i][j].data.sock , tSolverInfo + processPartition[i].startDepth , sizeof(SolverInfo<PixelChannels>) * (processPartition[i].endDepth-processPartition[i].startDepth+1) , "Failed to get process data from process" );
				for( int d=processPartition[i].startDepth ; d<=processPartition[i].endDepth ; d++ )
				{
//printf( "%d %d %f %f %f\n" , i , d , tSolverInfo[d].bSquareNorm , tSolverInfo[d].rSquareNorm , tSolverInfo[d].xSquareNorm );
					solverInfo[d].bSquareNorm += tSolverInfo[d].bSquareNorm;
					solverInfo[d].rSquareNorm += tSolverInfo[d].rSquareNorm;
					solverInfo[d].xSquareNorm += tSolverInfo[d].xSquareNorm;
					for( int c=0 ; c<PixelChannels ; c++ ) solverInfo[d].solutionSum[c] += tSolverInfo[d].solutionSum[c];
				}
			}
		for( int j=0 ; j<processPartition[0].size() ; j++ )
		{
			double avg[PixelChannels];
			ReceiveOnSocket ( processPartition[0][j].data.sock , GetPointer( avg ) , sizeof( avg ) , "Failed to get average from process" );
			for( int c=0 ; c<PixelChannels ; c++ )	_average[c]+=avg[c];
		}
		if(ii==0)	for( int c=0 ; c<PixelChannels ; c++ )	_average[c] /= _globalData.cWidth	,	_average[c] /= _globalData.cHeight;
		printf("Client Restriction:    %f\n",Time()-t) , fflush( stdout );
		if(_globalData.verbose)
			for( int i=0 ; i<processPartition.size() ; i++ )
				for( int d=processPartition[i].startDepth ; d<=processPartition[i].endDepth ; d++ )
				{
					int dd = processPartition[i].endDepth - ( d - processPartition[i].startDepth );
					printf("\tError[%d x %d] %g -> %g\t(%g)\n" , _globalData.width >> d , _globalData.height >> d , sqrt(solverInfo[dd  ].bSquareNorm) , sqrt(solverInfo[dd  ].rSquareNorm) , sqrt(solverInfo[dd  ].xSquareNorm) ) , fflush( stdout );
				}

		// Get the base solution
		DotProductStencil dotMajor , d2DotMajor , dotMinor ,d2DotMinor;
		SetDownSampledStencil( _globalData.width  , _oocLevels , dotMajor , d2DotMajor );
		SetDownSampledStencil( _globalData.height , _oocLevels , dotMinor , d2DotMinor );
		if( ii==_globalData.vCycles-1)	SolveInCore( dotMajor , d2DotMajor , dotMinor ,d2DotMinor , globalB , globalX , _average);
		else							SolveInCore( dotMajor , d2DotMajor , dotMinor ,d2DotMinor , globalB , globalX , zeroAverage);

		t=Time();
		int gWidth2 = gWidth<<1 , gHeight2 = gHeight<<1;
		outX->reset( false , 1 );
		// Get the up-sampled solution
		float *upX = new float[ PixelChannels * gWidth2 * gHeight2 ];

		for( int c=0 ; c<PixelChannels ; c++ )
		{
			Vector< float > low , high;
			low.Resize( gWidth * gHeight );
			for( int j=0 ; j<gHeight  ; j++ ) memcpy( &low[ j*gWidth ] , &globalX[ j*gWidth*PixelChannels + gWidth*c ] ,  sizeof(float) * gWidth );
			if     ( _globalData.periodicType==  SPHERICAL_PERIODIC ) UpSampleSphericalSolution  ( low , gWidth , gHeight , high );
			else if( _globalData.periodicType==CYLINDRICAL_PERIODIC ) UpSampleCylindricalSolution( low , gWidth , gHeight , high );
			else                                                      UpSampleRectangularSolution( low , gWidth , gHeight , high );
			for( int j=0 ; j<gHeight2 ; j++ ) memcpy( &upX[ j*gWidth2*PixelChannels + c*gWidth2 ] , &high[ j*gWidth2 ] ,  sizeof(float) * gWidth2 );
		}

		// Send the up-sampled solution to the socket
		for( int j=0 ; j<gHeight2 ; j++ )
		{
			Pointer( SyncType ) row = ( Pointer( SyncType ) )(*outX)[j];
//			memcpy( row , &upX[gWidth2*j*Channels] , sizeof(float)*Channels*gWidth2 );
			for( int i=0 ; i<PixelChannels*gWidth2 ; i++ ) row[i] = SyncType( upX[gWidth2*j*PixelChannels+i] );
			outX->advance();
		}
		delete upX;

		memset( solverInfo , 0 , sizeof(SolverInfo< PixelChannels >) * _oocLevels );
		for( int i=0 ; i<processPartition.size() ; i++ )
			for( int j=0 ; j<processPartition[i].size() ; j++ )
			{
				ReceiveOnSocket ( processPartition[i][j].data.sock , tSolverInfo + processPartition[i].startDepth , sizeof(SolverInfo<PixelChannels>) * (processPartition[i].endDepth-processPartition[i].startDepth+1) , "Failed to get process data from process" );
				for( int d=processPartition[i].startDepth ; d<=processPartition[i].endDepth ; d++ )
				{
					solverInfo[d].bSquareNorm += tSolverInfo[d].bSquareNorm;
					solverInfo[d].rSquareNorm += tSolverInfo[d].rSquareNorm;
					solverInfo[d].xSquareNorm += tSolverInfo[d].xSquareNorm;
					for( int c=0 ; c<PixelChannels ; c++ ) solverInfo[d].solutionSum[c] += tSolverInfo[d].solutionSum[c];
				}
			}
		printf("Client Prolongation:    %f\n",Time()-t) , fflush( stdout );
		if( _globalData.verbose )
			for( int i=processPartition.size()-1 ; i>=0 ; i-- )
				for( int d=processPartition[i].endDepth ; d>=processPartition[i].startDepth ; d-- )
				{
					int dd = processPartition[i].startDepth - ( d - processPartition[i].endDepth );
					printf("\tError[%d x %d] %g -> %g\n" , _globalData.width >> d , _globalData.height >> d , sqrt(solverInfo[dd  ].bSquareNorm) , sqrt(solverInfo[dd  ].rSquareNorm)) , fflush( stdout );
				}
	}
	FreeArray( solverInfo  );
	FreeArray( tSolverInfo );
	delete inB;
	delete outX;
}
template< int PixelChannels , int LabelChannels , class SyncType >
void SocketedMultigridServer< PixelChannels , LabelChannels , SyncType >::SolveInCore(DotProductStencil& dotMajor,DotProductStencil& d2DotMajor,
																 DotProductStencil& dotMinor,DotProductStencil& d2DotMinor,
																 Vector<float>& in,Vector<float>& out,double average[PixelChannels])
{
	double t;
	int w = _globalData.width >> _oocLevels , h = _globalData.height >> _oocLevels;
	SocketedMultiGridStreamingSolver< PixelChannels , float , float >* solvers;
	DotProductStencil lowDotMajor , lowD2DotMajor , lowDotMinor , lowD2DotMinor;

	SetDownSampledStencil( dotMajor , d2DotMajor , _globalData.width  >> _oocLevels , _icLevels-1 , lowDotMajor , lowD2DotMajor );
	SetDownSampledStencil( dotMinor , d2DotMinor , _globalData.height >> _oocLevels , _icLevels-1 , lowDotMinor , lowD2DotMinor );
	Vector<float> lowB,lowX;
	StreamingGrid *B,*X;
	solvers=new SocketedMultiGridStreamingSolver< PixelChannels , float , float >[_icLevels];

	for( int i=1 ; i<_icLevels   ; i++ )	solvers[i].parent=&solvers[i-1];
	for( int i=0 ; i<_icLevels-1 ; i++ )	solvers[i].rChild =solvers[i].pChild =&solvers[i+1];
	for( int i=0 ; i<_icLevels   ; i++ )
	{
		solvers[i].bSquareNorm=solvers[i].rSquareNorm=solvers[i].xSquareNorm=0;
		for( int c=0 ; c<PixelChannels ; c++ ) solvers[i].solutionSum[c]=0;
		solvers[i].setResidual = true;
	}
	// BADNESS!!! For spherical domains, the server socket has not been set here.
#if LOW_MEMORY_HACK // WARNING!!!: Big hack over here: To avoid memory overflow we are using an ooc solver here
	MultiStreamIOServer server;
	if( _globalData.iWeight )	solvers[_icLevels-1].Initialize( dotMajor , d2DotMajor , dotMinor , d2DotMinor , _globalData.iWeight , 0 , w , w , h , _globalData.iters , _inCoreConnectInfo.left , _icSyncSockets , _inCoreConnectInfo.right , true , _globalData.spherical , &server );
	else						solvers[_icLevels-1].Initialize( dotMajor , d2DotMajor , dotMinor , d2DotMinor ,                       0 , w , w , h , _globalData.iters , _inCoreConnectInfo.left , _icSyncSockets , _inCoreConnectInfo.right , true , _globalData.spherical , &server );
#else // !LOW_MEMORY_HACK
	if( _globalData.iWeight )	solvers[_icLevels-1].Initialize( dotMajor , d2DotMajor , dotMinor , d2DotMinor , _globalData.iWeight , 0 , w , w , h , _globalData.iters , _inCoreConnectInfo.left , _icSyncSockets , _inCoreConnectInfo.right , false , _globalData.periodicType , NULL );
	else						solvers[_icLevels-1].Initialize( dotMajor , d2DotMajor , dotMinor , d2DotMinor ,                       0 , w , w , h , _globalData.iters , _inCoreConnectInfo.left , _icSyncSockets , _inCoreConnectInfo.right , false , _globalData.periodicType , NULL );
#endif // LOW_MEMORY_HACK

	lowX.Resize( solvers[0].major*solvers[0].minor*PixelChannels );
	lowB.Resize( solvers[0].major*solvers[0].minor*PixelChannels );
	B = new MemoryBackedGrid( ( Pointer( byte ) )GetPointer( lowB[0] , (int)( lowB.Dimensions() ) ) , solvers[0].major * PixelChannels * sizeof(float) , solvers[0].minor , true );
	X = new MemoryBackedGrid( ( Pointer( byte ) )GetPointer( lowX[0] , (int)( lowX.Dimensions() ) ) , solvers[0].major * PixelChannels * sizeof(float) , solvers[0].minor , true );
	solvers[_icLevels-1].inB = new MemoryBackedGrid( ( Pointer( byte ) )GetPointer(  in[0] , (int)(  in.Dimensions() ) ) , w*PixelChannels*sizeof(float) , h , true );
	solvers[_icLevels-1].inX = new MemoryBackedGrid( ( Pointer( byte ) )GetPointer( out[0] , (int)( out.Dimensions() ) ) , w*PixelChannels*sizeof(float) , h , true );

	solvers[0].outX = X;
	solvers[0].outB = B;

	// Data for the interleaved streaming multigrid
	t=Time();
	// Initialize
	solvers[_icLevels-1].InitRestriction();
	solvers[_icLevels-1].SetRestriction();
#if LOW_MEMORY_HACK
	server.Reset();
	server.StartIO();
#endif // LOW_MEMORY_HACK
	solvers[_icLevels-1].SolveRestriction();
#if LOW_MEMORY_HACK
	server.Reset();
#endif // LOW_MEMORY_HACK

	// Clean up
	delete solvers[_icLevels-1].inB;
	solvers[_icLevels-1].inB=NULL;
	delete solvers[_icLevels-1].inX;
	solvers[_icLevels-1].inX=NULL;

	printf( "Server Restriction:        %f\n" , Time()-t ) , fflush( stdout );
	if(_globalData.verbose)
		for( int i=_icLevels-1 ; i>=0 ; i-- )
			printf("\tError[%d x %d] %g -> %g\n",solvers[i].major,solvers[i].minor,sqrt(solvers[i].bSquareNorm),sqrt(solvers[i].rSquareNorm)) , fflush( stdout );

	solvers[_icLevels-1].UnSetRestriction();
	// Get the base solution
	SparseMatrix<double> lap;
	t=Time();
	if     ( _globalData.periodicType==SPHERICAL_PERIODIC )
		if( _globalData.iWeight==0 )	SetSphericalLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor);
		else							SetSphericalLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor,_globalData.iWeight);
	else if( _globalData.periodicType==CYLINDRICAL_PERIODIC )
		if( _globalData.iWeight==0 )	SetCylindricalLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor);
		else							SetCylindricalLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor,_globalData.iWeight);
	else
		if( _globalData.iWeight==0)	SetRectangularLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor);
		else						SetRectangularLaplacianMatrix(lowDotMajor,lowD2DotMajor,lowDotMinor,lowD2DotMinor,lap,solvers[0].major,solvers[0].minor,_globalData.iWeight);
	{
		Vector<double> myLowX,myLowB;
		myLowB.Resize(solvers[0].major*solvers[0].minor);
		myLowX.Resize(solvers[0].major*solvers[0].minor);
		for(int c=0;c<PixelChannels;c++)
		{
			for(int i=0;i<solvers[0].major;i++)
				for(int j=0;j<solvers[0].minor;j++)
				{
					myLowX[i+j*solvers[0].major] = lowX[i+j*solvers[0].major*PixelChannels+c*solvers[0].major];
					myLowB[i+j*solvers[0].major] = lowB[i+j*solvers[0].major*PixelChannels+c*solvers[0].major];
				}
			myLowB = myLowB - lap * myLowX;
			myLowX.SetZero();
		if(!_noCG)
#if NEW_MULTIPLY
			if( _globalData.iWeight==0 ) MySolveConjugateGradient< true  >( lap , myLowB , 20*(int)( sqrt(lap.groups+1.0) ) , myLowX );
			else						 MySolveConjugateGradient< false >( lap , myLowB , 20*(int)( sqrt(lap.groups+1.0) ) , myLowX );
#else // !NEW_MULTIPLY
			if( _globalData.iWeight==0 ) MySolveConjugateGradient< Type >( lap , myLowB , 20*(int)( sqrt(lap.groups+1.0) ) , myLowX );
			else						 MySolveConjugateGradient< Type >( lap , myLowB , 20*(int)( sqrt(lap.groups+1.0) ) , myLowX );
#endif // NEW_MULTIPLY

			for(int i=0;i<solvers[0].major;i++)
				for(int j=0;j<solvers[0].minor;j++)
				{
					lowX[i+j*solvers[0].major*PixelChannels+c*solvers[0].major] = float( myLowX[i+j*solvers[0].major] );
					myLowX[i+j*solvers[0].major] = 0.0;
				}
		}
	}
	printf( "Server Conjugate Gradient: %f\n" , Time()-t ) , fflush( stdout );

	for( int i=0 ; i<_icLevels ; i++ )
	{
		solvers[i].bSquareNorm = solvers[i].rSquareNorm = solvers[i].xSquareNorm=0;
		for( int c=0 ; c<PixelChannels ; c++ ) solvers[i].solutionSum[c] = 0;
		solvers[i].setResidual = _globalData.verbose;
	}
	// Solve the prolongation
	t=Time();
	solvers[0].inX = X;
	// Set the child dependencies
	solvers[_icLevels-1].outX = new MemoryBackedGrid( ( Pointer( byte ) )GetPointer( out[0] , (int)( out.Dimensions() ) ) , w*PixelChannels*sizeof(float) , h , true );
	solvers[_icLevels-1].InitProlongation();
	solvers[0].SetProlongation();
	// Solve
#if LOW_MEMORY_HACK
	server.StartIO();
#endif // LOW_MEMORY_HACK
	solvers[0].SolveProlongation();
#if LOW_MEMORY_HACK
	server.Reset();
#endif // LOW_MEMORY_HACK
	printf( "Server Prolongation:       %f\n" , Time()-t ) , fflush( stdout );
	if(_globalData.verbose)
	{
		for(int i=0;i<_icLevels;i++)
		{
			printf("\tError[%d x %d] %g -> %g\n",solvers[i].major,solvers[i].minor,sqrt(solvers[i].bSquareNorm),sqrt(solvers[i].rSquareNorm)) , fflush( stdout );
		}
	}
	solvers[0].UnSetProlongation();
	delete solvers[_icLevels-1].outX;
	delete[] solvers;
	delete B;
	delete X;
	B=X=NULL;

	int ww=( (long long)(_globalData.cWidth) *w+_globalData.width -1)/_globalData.width;
	int hh=( (long long)(_globalData.cHeight)*h+_globalData.height-1)/_globalData.height;
	for(int c=0;c<PixelChannels;c++)
	{
		double newAverage=0;
		int pCount=0;

		for(int j=0;j<hh;j++)
		{
			float* o=&out[j*w*PixelChannels+c*w];
			for(int i=0;i<ww;i++)
			{
				newAverage+=o[i];
				pCount++;
			}
		}

		newAverage/=pCount;
		newAverage=average[c]-newAverage;

		if( _globalData.iWeight==0 )
			for(int j=0;j<h;j++)
			{
				float* o=&out[j*w*PixelChannels+c*w];
				for( int i=0 ; i<w ; i++ ) o[i] += float( newAverage );
			}
	}
}

//////////////////////////////////
// SocketedSuperMultigridServer //
//////////////////////////////////
bool SocketedSuperMultigridServer::SetUp( char* prefix , int port , int clientCount , int iters , int inCoreRes , int minMGRes , int vCycles , int minBandSize ,
												    int tileWidth , int tileHeight , const char* tileExt , 
												    bool gammaCorrection , int quality , int lanes , bool verbose , int periodicType , double iWeight , double gScale ,
													bool removeAverageGradient , int unknownType , bool showProgress ,
#if DEBUG_FLAGS
													DebugFlags debugFlags ,
#endif // DEBUG_FLAGS
												    bool noCG , bool shortSync )
{
#if DEBUG_FLAGS
	memcpy( &(this->debugFlags ) , &debugFlags , sizeof(debugFlags) );
#endif // DEBUG_FLAGS
	_noCG = noCG;
	_port = port;
	_clientCount = clientCount;
	_inCoreRes = inCoreRes;
	_minMGRes = minMGRes;
	_minBandSize = minBandSize;
	if( _minBandSize > _inCoreRes ) _inCoreRes = _minBandSize;
	_removeAverageGradient = removeAverageGradient;
	_globalData.iters = iters;
	_globalData.vCycles = vCycles;
	_globalData.tileWidth = tileWidth;
	_globalData.tileHeight = tileHeight;
	if( tileExt ) strcpy( _globalData.tileExt , tileExt );
	else _globalData.tileExt[0] = 0;
	_globalData.gammaCorrection = gammaCorrection;
	_globalData.quality = quality;
	_globalData.lanes = lanes;
	_globalData.verbose = verbose;
	_globalData.periodicType = periodicType;
	_globalData.iWeight = iWeight;
	_globalData.gScale = gScale;
	_globalData.unknownType = unknownType;
	_globalData.showProgress = showProgress;
	_globalData.shortSync = shortSync;

	GetHostAddress( _address , prefix );

	return true;
}
template< int PixelChannels , int LabelChannels , class SyncType , class LabelType >
bool SocketedSuperMultigridServer::Run( void )
{
	int subClientCount = 0;

	// Create a listening SOCKET for connecting to server
	AcceptorSocket listenSocket = _INVALID_ACCEPTOR_SOCKET_;
	listenSocket = GetListenSocket( _port );
	if ( listenSocket == _INVALID_ACCEPTOR_SOCKET_ ) return false;
	printf( "Server Address: %s:%d\n", _address , _port ) , fflush( stdout );

	ClientSocket* clientSockets = new ClientSocket[_clientCount];
	int coreCount = 0;

	// Establish a connection to the clients
	for( int i=0 ; i<_clientCount ; i++ )
	{
		clientSockets[i].socket  = AcceptSocket( listenSocket );
		printf( "Connected to process: %d /%d     \r" , i+1 , _clientCount );
	}
	printf( "\n" ) , fflush( stdout );
	// Get the client info
	int width = 0 , height;
	for( int i=0 ; i<_clientCount ; i++ )
	{
		ReceiveOnSocket( clientSockets[i].socket , GetPointer( clientSockets[i].clientData ) , sizeof( ClientData ) , "Failed to get client data from client" );
		if(!i) height = clientSockets[i].clientData.height;
		else if( height!=clientSockets[i].clientData.height ) fprintf( stderr , "Band heights differ: %d != %d\n" , height , clientSockets[i].clientData.height ) , exit(0);
		width += clientSockets[i].clientData.width;
		subClientCount += clientSockets[i].clientData.subClients;
	}
	// Sort the clients by index order
	qsort( clientSockets , _clientCount , sizeof(ClientSocket) , ClientSocket::Sort );

	// Store the clipping width and then compute the padded dimensions
	_globalData.cWidth  = width;
	_globalData.cHeight = height;

#if 0
#pragma message( "[WARNING] For the life of me, why do I have to make this change" )
	_globalData.width = width , _globalData.height = height;
#else
	if     ( _globalData.periodicType==         NO_PERIODIC ) _globalData.width = width+1, _globalData.height = height+1;
	else if( _globalData.periodicType==CYLINDRICAL_PERIODIC ) _globalData.width = width  , _globalData.height = height+1;
	else                                                      _globalData.width = width  , _globalData.height = height;
#endif
	int w = _globalData.width , h = _globalData.height , wBS , hBS;
	SetPaddedSize( w , h , _minMGRes , _globalData.width , _globalData.height , wBS , hBS );
	if( _globalData.periodicType==CYLINDRICAL_PERIODIC && _globalData.width!=w ) fprintf( stderr , "Cylindrical width was not preserved: (%d , %d) -> (%d , %d)\t(%d , %d)\n" , w , h , _globalData.width , _globalData.height , wBS , hBS ) , exit( 0 );
	if( _globalData.periodicType==  SPHERICAL_PERIODIC && ( _globalData.width!=w || _globalData.height!=h ) ) fprintf( stderr , "Spherical dimensions were not preserved: (%d , %d) -> (%d , %d)\t(%d , %d)\n" , w , h , _globalData.width , _globalData.height , wBS , hBS ) , exit( 0 );
	int oocLevels = 1;

	
#if 1 // QUESTION: Shouldn't this computed relative to the globalData.width/height?
	while( (long long)( _globalData.width >> oocLevels ) * (long long)( _globalData.height >> oocLevels ) >= _inCoreRes*_inCoreRes ) oocLevels++;
	if( !ProcessPartition<int>::IsValidBandSize( _globalData.width>>oocLevels , _globalData.height>>oocLevels , _globalData.iters , _minBandSize ) )
	{
		fprintf( stderr , "In-core resolution is too small: %d x %d\n" , _globalData.width>>oocLevels , _globalData.height>>oocLevels );
		return false;
	}
#else
	while( (long long)( width >> oocLevels ) * (long long)( height >> oocLevels ) >= _inCoreRes*_inCoreRes ) oocLevels++;
	if( !ProcessPartition<int>::IsValidBandSize( width>>oocLevels , height>>oocLevels , _globalData.iters , _minBandSize ) )
	{
		fprintf( stderr , "In-core resolution is too small: %d x %d\n" , width>>oocLevels , height>>oocLevels );
		return false;
	}
#endif
	int blockSize = 4<<oocLevels;
#if 1	// Since the coarsest level is a multiple of four, and in the next finer one pixels have a width of two, we only need to
		// double that to ensure multipler of four.
	blockSize >>= 1;
#endif

	for( int i=0 ; i<_clientCount-1 ; i++ )
		if( clientSockets[i].clientData.width % blockSize )
		{
			fprintf( stderr , "Non-paddable client[%d] width not multiple of block size: %d mod %d != 0 \n" , i , clientSockets[i].clientData.width , blockSize );
			exit(0);
		}
	printf( "  Image Dimension: %d x %d -> %d x %d\n" , width , height , _globalData.width , _globalData.height ) , fflush( stdout );
	printf( "In-Core Dimension: %d x %d\n" , _globalData.width>>oocLevels , _globalData.height>>oocLevels ) , fflush( stdout );
	printf( "       Block Size: %d\n" , blockSize ) , fflush( stdout );

	// Set the local information for the bands
	int start = 0 , index = 0;
	for( int i=0 ; i<_clientCount ; i++ )
	{
		int w = clientSockets[i].clientData.width;
		clientSockets[i].clientData.index = index;
		index += clientSockets[i].clientData.subClients;
		clientSockets[i].clientData.start = start;
		clientSockets[i].clientData.end = clientSockets[i].clientData.cEnd = start+w;
		if( i==_clientCount-1 ) clientSockets[i].clientData.end  = _globalData.width;
		start = clientSockets[i].clientData.end;
	
		SendOnSocket( clientSockets[i].socket , GetPointer( _globalData ) , sizeof( _globalData ), "Failed to send global data to client" );
		SendOnSocket( clientSockets[i].socket , GetPointer( clientSockets[i].clientData ) , sizeof( ClientData ) ,"Failed to send client data to client" );
		SendOnSocket( clientSockets[i].socket , GetPointer( blockSize ), sizeof( blockSize ) , "Failed to send block size to client" );
		int getGradientAverage = _removeAverageGradient ? 1 : 0;
		SendOnSocket( clientSockets[i].socket , GetPointer( getGradientAverage ) , sizeof( getGradientAverage ) , "Failed to send gradient average status to client" );
		int useGrayImage = PixelChannels==3 ? 0 : 1;
		SendOnSocket( clientSockets[i].socket , GetPointer( useGrayImage ) , sizeof( useGrayImage ) , "Failed to send use gray status to client" );
#if DEBUG_FLAGS
		SendOnSocket( clientSockets[i].socket , GetPointer( debugFlags ) , sizeof( debugFlags ) , "Failed to send debug flags to client" );
#endif // DEBUG_FLAGS
	}
	delete[] clientSockets;
	int clipDimensions[2];
	clipDimensions[0] = _globalData.cWidth;
	clipDimensions[1] = _globalData.cHeight;
	SocketedMultigridServer< PixelChannels , LabelChannels , SyncType > server;
	if( !server.template SetUp< LabelType >
		(
		_address , listenSocket , subClientCount , _globalData.iters , _inCoreRes , _minMGRes , _globalData.vCycles , _minBandSize ,
		_globalData.tileWidth , _globalData.tileHeight , _globalData.tileExt ,
		_globalData.gammaCorrection , _globalData.quality , _globalData.lanes , _globalData.verbose , _globalData.periodicType , _globalData.iWeight , _globalData.gScale ,
		_removeAverageGradient , _globalData.unknownType , _globalData.showProgress ,
		_noCG , _globalData.shortSync , clipDimensions )
		) return false;
	server.Run();
	return true;
}