#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

// Simulates heat transmission through a plate represented by a matrix
typedef struct
{
	// Dimensions of the matrix
	size_t rows;
	size_t cols;
	// When temperature changes are lower than epsilon, simulation stops
	double epsilon;
	// Matrix of temperatures representing the plate
	double** matrix;
	// An auxiliary matrix to store the state of the plate after an update
	double** mirror;
	// The starting time of the simulation
	struct timespec start_time;
	// Current generation of the simulation
	size_t generation;
	// Pointer to the matrix with the current state of the simulation
	double** current;
	// Pointer to the matrix that will hold the next generation
	double** next;
} heat_t;

int read_matrix(heat_t* heat, const char* filename);
void simulate(heat_t* heat);
double update_generation(heat_t* heat);
double update_cell(heat_t* heat, size_t row, size_t col);
double elapsed(const struct timespec* start_time);
void** create_matrix(const size_t rows, const size_t columns, const size_t cell_size);
void** destroy_matrix(void** matrix, const size_t rows);

int main(int argc, char* argv[])
{
	if ( argc != 3 )
		return fprintf(stderr, "usage: heat matrix.bin epsilon\n"), 1;

	heat_t heat;
	memset(&heat, 0, sizeof(heat_t));

	heat.epsilon = strtod(argv[2], NULL);
	if ( heat.epsilon <= 0.0 )
		return fprintf(stderr, "error: invalid epsilon %lg\n", heat.epsilon), 2;

	int error = read_matrix(&heat, argv[1]);
	if ( error == EXIT_SUCCESS )
	{
		clock_gettime(CLOCK_MONOTONIC, &heat.start_time);
		simulate(&heat);

		destroy_matrix((void**)heat.matrix, heat.rows);
		destroy_matrix((void**)heat.mirror, heat.rows);
	}

	return error;
}

int read_matrix(heat_t* heat, const char* filename)
{
	int error = EXIT_SUCCESS;

	FILE* input = fopen(filename, "rb");
	if ( input == NULL )
		return fprintf(stderr, "error opening %s\n", filename), 11;

	if ( fread(&heat->rows, sizeof(size_t), 1, input) > 0 &&
		 fread(&heat->cols, sizeof(size_t), 1, input) > 0 &&
		 heat->rows >= 2 && heat->cols >= 2 )
	{
		heat->matrix = (double**) create_matrix(heat->rows, heat->cols, sizeof(double));
		heat->mirror = (double**) create_matrix(heat->rows, heat->cols, sizeof(double));

		if ( heat->matrix && heat->mirror )
		{
			for ( size_t row = 0; row < heat->rows; ++row )
			{
				if ( fread( heat->matrix[row], sizeof(double), heat->cols, input ) > 0 )
					memcpy( heat->mirror[row], heat->matrix[row], sizeof(double) * heat->cols );
				else
				{
					fprintf(stderr, "error reading %zu row\n", row + 1);
					error = 14;
					break;
				}
			}
		}
		else
		{
			fprintf(stderr, "error allocating %zu x %zu matrixes\n", heat->rows, heat->cols);
			error = 13;
		}
	}
	else
	{
		fprintf(stderr, "error reading matrix dimensions\n");
		error = 12;
	}

	if ( error )
	{
		heat->matrix = (double**) destroy_matrix((void**)heat->matrix, heat->rows);
		heat->mirror = (double**) destroy_matrix((void**)heat->mirror, heat->rows);
	}
	else
		printf("%zu x %zu matrix loaded (%0.3lf GiB)\n", heat->rows, heat->cols
			, 2 * 8 * heat->rows * heat->cols / 1073741824.0);

	fclose(input);
	return error;
}

void simulate(heat_t* heat)
{
	heat->current = heat->matrix;
	heat->next = heat->mirror;

	double delta = 0.0;
	while ( (delta = update_generation(heat)) > heat->epsilon )
		printf("%.9lfs: generation %zu delta=%0.6lf\n", elapsed(&heat->start_time), heat->generation, delta);
}

double update_generation(heat_t* heat)
{
	// The maximum temperature change detected in this generation
	double max_change = DBL_MIN;

	// Update all rows
	for ( size_t row = 1; row < heat->rows - 1; ++row )
	{
		// Upate all columns
		for ( size_t col = 1; col < heat->cols - 1; ++col )
		{
			// Update cell from current matrix to the next matrix
			double change = update_cell(heat, row, col);

			// If the cell change is greater than our maximum, update the maximum
			if ( change > max_change )
				max_change = change;
		}
	}

	// The next matrix is enterely updated, it is now the current matrix
	// The current matrix will be used for the next generation. So swap them
	double** temp = heat->current;
	heat->current = heat->next;
	heat->next = temp;

	// A new generation was completed
	++heat->generation;

	// Return the max temperature change detected in this generation
	return max_change;
}

double update_cell(heat_t* heat, size_t row, size_t col)
{
	// The new value is the average of the neighbors in the previous generation
	heat->next[row][col] = ( heat->current[row - 1][col]
		+ heat->current[row][col - 1]
		+ heat->current[row][col + 1]
		+ heat->current[row + 1][col] ) / 4.0;

	// Return the difference between the cell's new temperature and old one
	return fabs( heat->next[row][col] - heat->current[row][col] );
}

double elapsed(const struct timespec* start_time)
{
	struct timespec current_time;
	clock_gettime(CLOCK_MONOTONIC, &current_time);

	return current_time.tv_sec - start_time->tv_sec
		+ 1e-9 * (current_time.tv_nsec - start_time->tv_nsec);
}

void** create_matrix(const size_t rows, const size_t columns, const size_t cell_size)
{
	void** matrix = (void**) calloc(rows, sizeof(void*));
	if ( matrix == NULL )
		return NULL;

	for ( size_t row = 0; row < rows; ++row )
		if ( ( matrix[row] = calloc( columns, cell_size ) ) == NULL )
			return destroy_matrix(matrix, rows), NULL;

	return matrix;
}

void** destroy_matrix(void** matrix, const size_t rows)
{
	if ( matrix )
		for ( size_t row = 0; row < rows; ++row )
			free( matrix[row] );

	free( matrix );
	return NULL;
}