use strict;
use Time::HiRes qw(gettimeofday tv_interval sleep);
my $Start_time = [gettimeofday];
use Term::ANSIScreen qw/:color :cursor :screen/;
use Win32::Console::ANSI qw/ Cursor /;
use threads;
use threads::shared;
use constant FATAL => "\aFATAL ERROR: ";
#===== Copyright 2009, Webpraxis Consulting Ltd. - ALL RIGHTS RESERVED - Email: webpraxis@gmail.com ============================
# liebmann_threaded.pl: Threaded solution of Poisson's equation on a rectangular region using Liebmann's method.
#===============================================================================================================================
#           Usage : perl liebmann_threaded.pl NOBLOCKS SLOW
#       Arguments : NOBLOCKS = number of domain partitions.
#                   SLOW     = boolean flag for slowing down the iteration process
#     Input Files : None.
#    Output Files : None.
# Temporary Files : None.
#         Remarks : See http://www.webpraxis.ab.ca/poisson/liebmann_threaded.shtml for details.
#         History : v1.0.0 - January 20, 2009 - Original release.
#===============================================================================================================================
# 0) DEFINE THE PROBLEM PARTICULARS:
my $F				= sub { my( $x, $y ) = @_; return $x * exp $y; };					#anonymous sub for rhs of Poisson equation
my $G				= sub { my( $x, $y ) = @_; return $x * exp $y; };					#anonymous sub for boundary conditions
my $Xa				= 0.;																#domain boundary: lower endoint on the x-axis
my $Xb				= 2.;																#domain boundary: upper endoint on the x-axis
my $Yc				= 0.;																#domain boundary: lower endoint on the y-axis
my $Yd				= 1.;																#domain boundary: upper endoint on the y-axis
my $N				= 6;																#number of grid columns
my $M				= 5;																#number of grid rows
my $Tol				= 1.e-10;															#tolerance in successive estimates
my $MaxIterations	= 100;																#maximum number of iterations to attempt
#-------------------------------------------------------------------------------------------------------------------------------
# 1) INITIALIZE:
$| = 1; 																				#set STDOUT buffer to auto-flush
cls();																					#clear screen
print colored ['black on white'], "$0\n\n\n", color 'reset';							#display program name

my $NoBlocks	= shift || die FATAL, 'Number of partitions not specified';				#get number of domain blocks
my $Slow		= shift;																#get flag for slowing down the iterations

my $H			= ( $Xb - $Xa ) / $N;													#compute mesh size along x-axis
my $K			= ( $Yd - $Yc ) / $M;													#compute mesh size along y-axis
my $H2			= $H * $H;																#define computational constant
my $Lambda		= $H2 / ( $K * $K );													#define computational constant
my $Mu			= 2. * ( 1. + $Lambda );												#define computational constant

print " Problem Details:\n\n", 															#echo problem details
	  "\t           Grid Size: N = $N, M = $M => ", ($N-1) * ($M-1), " interior mesh points\n",
	  "\t    No of partitions: $NoBlocks\n",
	  "\t                  Xa: $Xa\n",
	  "\t                  Xb: $Xb\n",
	  "\t                   H: $H\n",
	  "\t                  Yc: $Yc\n",
	  "\t                  Yd: $Yd\n",
	  "\t                   K: $K\n",
	  "\t           Tolerance: $Tol\n",
	  "\tMax no of iterations: $MaxIterations\n",
	  "\t         Slow Motion: ", ($Slow) ? 'Yes' : 'No', "\n\n",
	  '-' x 80, "\n\n";
die FATAL, "N must exceed the number of partitions" unless $N > $NoBlocks;				#check for minimum grid size
die FATAL, "M must exceed 1" unless $M > 1;
#-------------------------------------------------------------------------------------------------------------------------------
# 2) DEFINE THE INTERIOR MESH POINTS, COMPUTE BOUNDARY VALUES & INITIALIZE SOLUTION ESTIMATE:
my @X :shared;																			#array for the interior mesh point x-values
my @Y :shared;																			#array for the interior mesh point y-values
my @W :shared;																			#array for the boundary values & solution estimates
   $W[$_] = &share([]) for ( 0..$N );													#extend as array of arrays by adding shared leaf nodes

$X[$_] = $Xa + $_ * $H for ( 1..$N - 1 );												#compute interior mesh point x-values
$Y[$_] = $Yc + $_ * $K for ( 1..$M - 1 );												#compute interior mesh point y-values

for ( 1..$N - 1 ) {																		#repeat for each interior horizontal node value
	$W[$_][$M] = &$G( $X[$_], $Yd );													# compute top boundary condition
	$W[$_][ 0] = &$G( $X[$_], $Yc );													# compute bottom boundary condition
}																						#until all values processed

for ( 1..$M - 1 ) {																		#repeat for each interior vertical node value
	$W[ 0][$_] = &$G( $Xa, $Y[$_] );													# compute left boundary condition
	$W[$N][$_] = &$G( $Xb, $Y[$_] );													# compute right boundary condition
}																						#until all values processed

for my $i ( 1..$N - 1 ) {																#repeat for each interior mesh point
	$W[$i][$_] = 0.	for ( 1..$M - 1 );													# init all estimates
}																						#until all interior mesh points processed
#-------------------------------------------------------------------------------------------------------------------------------
# 3) PARTITION DOMAIN INTO VERTICAL STRIPS; TRY TO BALANCE THE IMPENDING WORK LOAD:
my @Blocks;																				#array of hashes for the block node boundaries
{																						#begin naked block
	my $blockIdx	= 0;																# block index
	my @noVLines	= ( int( ($N - 1) / $NoBlocks ) ) x $NoBlocks;						# init number of vertical grid lines per block

	PARTITIONS: { 																		# begin domain partitioning
		$Blocks[0 ] = { TOP => $M - 1, BOTTOM => 1, LEFT => 1, RIGHT => $noVLines[0] };
		$Blocks[$_] = { TOP => $M - 1, BOTTOM => 1, LEFT => $Blocks[$_-1]{RIGHT} + 1, RIGHT => $Blocks[$_-1]{RIGHT} + $noVLines[$_] }
		 for( 1..$NoBlocks - 1 );

		++$noVLines[$blockIdx++], redo PARTITIONS										#  increase width of one block at the time
		 if $Blocks[$#Blocks]{RIGHT} < $N - 1;											#   until all interior nodes allocated
	}																					# end domain partitioning
}																						#end naked block
#-------------------------------------------------------------------------------------------------------------------------------
# 4) ITERATE THE BLOCKS:
my $CursorY;																			#cursor row coordinate
my @Spinners		= ( '-', '\\', '|', '/' );											#define symbols for spinner
my @threads;																			#array of threads
																						#variables shared between calling program and threads:
my $MaxDev			:shared;															# maximum deviation in successive estimates
my $NoIterations	:shared;															# counter for number of iterations used
																						#variables shared between threads only:
my $Halt			:shared;															# boolean flag for halting the threads
my $NoBlocksDone	:shared;															# counter for number of partitions processed

print " Process Details:\n\n";															#inform user of processing start
$CursorY = ( Cursor() )[1];  															#record cursor row position

$threads[$_] = threads->create( \&estimateBlock, $_ ) for ( 0..$#Blocks );				#launch a thread for each partition
$threads[$_]->join() for ( 0..$#Blocks );												#wait for threads to return

print locate( ++$CursorY + $NoBlocks, 1 ), '-' x 80, "\n\n";
die	FATAL, "Iteration limit reached (MaxDev = $MaxDev).\n" if $MaxDev > $Tol;			#halt if iteration limit reached
#-------------------------------------------------------------------------------------------------------------------------------
# 5) OUTPUT THE ESTIMATES, EXACT VALUES AND ERROR NORMS:
print	" $NoIterations Iterations Required:\n\n";										#report no of iterations
printf	"\t%3s %3s      %1s        %1s        %1s        %1s      %7s\n",				#output column headings
		'i','j','x','y','w','u','|u - w|';
print	"\t", '-' x 55, "\n";
for my $i ( 1..$N - 1 ) {																#repeat for each interior horizontal node value
	for my $j ( 1..$M - 1 ) {															# repeat for each interior vertical node value
		my $u = &$F( $X[$i],$Y[$j] );													#  compute exact solution
		my $e = abs( $u - $W[$i][$j] );													#  compute error norm
		printf	"\t%3s %3s   %7.5f  %7.5f  %7.5f  %7.5f  %4.2e\n",						#  output values
				$i,$j,$X[$i],$Y[$j],$W[$i][$j],$u,$e;
	}																					# until all vertical values processed
}																						#until all horizontal values processed
print "\n\tElapsed time: ", tv_interval( $Start_time ), "\n";							#report total elapsed time
exit;																					#end processing
#===== SUBROUTINES =============================================================================================================
#            Usage : &estimateBlock( %BLOCKIDX );
#          Purpose : Estimate the solution for a specified partition using Liebmann's method.
#        Arguments : %BLOCKIDX = partition index value
#        Externals : $CursorY, $H2, $Lambda, $Mu, $NoBlocks, $Slow, @Blocks, @Spinners
# Shared Externals : $Halt, $MaxDev, $NoBlocksDone, $NoIterations, @W, @X, @Y
#             Subs : $F
#          Remarks : None.
#          History : v1.0.0 - January 20, 2009 - Original release.

sub estimateBlock {																		#begin sub
	my $blockIdx	= shift;															# parameterize the argument
	my %block		= %{ $Blocks[ $blockIdx ] };										# parameterize the block boundaries
	my $maximum		= sub { my( $x, $y ) = @_; return ($x < $y ) ? $y : $x; };			# anonymous sub for computing the maximum of 2 numbers
	my $spinnerIdx	= 0;																# init spinner index
	my $spinner;																		# spinner

	do {																				# repeat
		my $maxBlockDev	= 0;															#  init maximum deviation for the block
		for my $i ( $block{LEFT}..$block{RIGHT} ) {										#  repeat for each block horizontal node value
			for my $j ( $block{BOTTOM}..$block{TOP} ) {									#   repeat for each block vertical node value
				$maxBlockDev	= &$maximum( $maxBlockDev, &estimateNode( $i, $j ) );	#    estimate solution at node & update maximum deviation
				$spinnerIdx 	= ++$spinnerIdx % 4;									#    update spinner index
				$spinner 		= "(" . $Spinners[$spinnerIdx] . ")";					#    update processing spinner
				print	locate( $CursorY + $blockIdx, 3 ), clline, "Partition #$blockIdx: ",
						colored ['bold green'], "Iterating " , $spinner;
				sleep( 0.2 ) if $Slow;													#    sleep if slowdown requested
			}																			#   until all vertical values processed
		}																				#  until all horizontal values processed

		lock( $NoBlocksDone );															#  lock counter for no of blocks processed
		++$NoBlocksDone;																#  increment counter
		$MaxDev	= &$maximum( $MaxDev, $maxBlockDev );									#  update maximum deviation across all blocks
		if( $NoBlocksDone == $NoBlocks ) {												#  if this was the last block to be processed
			++$NoIterations;															#   update iteration counter
			print	locate( $CursorY + $blockIdx, 3 ), clline, "Partition #$blockIdx: ",
					colored ['bold yellow'], "Checking after $NoIterations iteration(s); MaxDev = $MaxDev";
			sleep( 1.4 ) if $Slow;														#   sleep if slowdown requested
			$Halt			= ($MaxDev <= $Tol) || ($NoIterations == $MaxIterations);	#   update stopping criterion
			$NoBlocksDone	= $MaxDev = 0 unless $Halt;									#   reset control variables if won't be stopping
			cond_broadcast( $NoBlocksDone );											#   signal other threads to resume
		} else {																		#  else
			print	locate( $CursorY + $blockIdx, 3 ), clline, "Partition #$blockIdx: ",
					colored ['bold red'], "Waiting...";
			cond_wait( $NoBlocksDone );													#   wait for remaining threads to finish estimating
		}																				#  end if-else
	} until $Halt;																		# until stopping criterion in effect
	print	locate( $CursorY + $blockIdx, 3 ), clline, "Partition #$blockIdx: ",
			colored ['bold red'], "Stopped";

	#            Usage : &estimateNode( I, J );
	#          Purpose : Estimates the solution at a specified grid node. Returns the deviation.
	#        Arguments : I = horizontal node value
	#                    J = vertical node value
	#        Externals : $H2, $Lambda, $Mu
	# Shared Externals : @W, @X, @Y
	#             Subs : $F
	#          History : v1.0.0 - January 19, 2009 - Original release.

	sub estimateNode {																	# begin sub
		my $i = shift;																	#  parameterize horizontal node value
		my $j = shift;																	#  parameterize vertical node value

		my $oldW	= $W[$i][$j];														#  save previous estimate at current node
		$W[$i][$j]	= ( - $H2 * &$F( $X[$i],$Y[$j] ) + $W[$i+1][$j] + $W[$i-1][$j]		#  compute new estimate
							+ $Lambda * ( $W[$i][$j+1] + $W[$i][$j-1] )
					  ) / $Mu;
		return abs( $W[$i][$j] - $oldW );												#  return the deviation
	}																					# end sub estimateNode
}																						#end sub estimateBlock
#===== Copyright 2009, Webpraxis Consulting Ltd. - ALL RIGHTS RESERVED - Email: webpraxis@gmail.com ============================
# end of liebmann_threaded.pl