LÖVE

Code: Select all

x = { r=2, i=-1}

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <mpi.h>
#include "plot.h"


#define THRESHOLD 1000.0
#define MAXITER 1500
#define SCALE 30000.0


double complex
    global_lower = -0.512  - 1.024*I,
    global_upper =  1.024  + 1.024*I;

/*
 * Estimate the divergence rate of our function for a given complex point
 */
double
test_point ( double complex c )
{
    double complex z_new, z, z_prev;
    double difference;

    z = z_prev = c;

    // Iterate until divergence or max. number of iterations
    for ( int iter=0; iter<MAXITER; iter++ )
    {
        // Update the value by the formula
        z_new = z*z + 0.56667 - (z_prev * 0.5);

        // Buffer the changes
        z_prev = z;
        z = z_new;

        // Test for "fast enough" divergence - if this is found, we can end
        // the iteration already, assuming that the sequence shortly will
        // disappear into infinity
        difference = cabs(z) - cabs(z_prev);
        if ( difference > THRESHOLD )
            break;
    }

    // Scale the difference value to make a nicer plot
    difference /= SCALE;
    return (difference > 1.0) ? 1.0 : difference;
}


/*
 * This function should evaluate the points in the rectangle bounded by
 * the lower/upper variables, discretized into a grid of w by h points.
 * When this is called, the 'data' pointer should already be appropriately
 * allocated to hold these w by h points.
 *
 * Note that while this array is one-dimensional, there is a convenience
 * macro AT(x,y) in 'plot.h' to address it as a 2D array (writing
 * data[AT(x,y)]).
 */
void
compute_region (
    double complex lower, double complex upper, int w, int h, double *data
)
{
   double complex delta = upper - lower;
   double real = creal(delta)/(double)w;
   double imag = cimag(delta)/(double)h;
   
   for(int i=0;i<w;i++)
   {
   	for(int j=0;j<h;j++)
   	{
   		double complex c = real * i + imag * j * I;
   		data[i*ROWS+j] = test_point( lower + c);
   	}
   }
}


int
main ( int argc, char **argv )
{
    double
        *img, *local_img,         // Global and local image parts
        t_start, t_finish, t_elapsed, max_elapsed;  // Timing variables

    complex double
        local_lower, local_upper, delta;    // Bounds and step size of domain

    int
        rank, size, blocksize;  // Used to compute problem partitioning

    // Initialize MPI, get world rank and size
    MPI_Init ( &argc, &argv );
    MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
    MPI_Comm_size ( MPI_COMM_WORLD, &size );

    // Rank 0 hosts the complete image
    if ( rank == 0 )
        img = calloc ( ROWS * COLS, sizeof(double) );

	// YOUR CODE HERE	
	
	// Assign local variables.
	delta = (creal(global_upper) - creal(global_lower))/(double)size + 0 * I;
	local_lower = global_lower + delta * (double)rank;
	local_upper = creal(local_lower+delta) + cimag(global_upper) * I;

	// Allocate memory.
	int memsize = (ROWS * COLS)/size;
	local_img = calloc(memsize, sizeof(double));
	
	// Compute region
	t_start = MPI_Wtime();
	compute_region(local_lower, local_upper, COLS/size, COLS, local_img);
	t_finish = MPI_Wtime();
	
	t_elapsed = t_finish - t_start;	

	// Gather at 0
	MPI_Gather(local_img, memsize, MPI_DOUBLE, img + memsize * rank, memsize, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	// Reduce
	MPI_Reduce(&t_elapsed, &max_elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

    // Rank 0 reports greatest time spent in computation, writes final image
    if ( rank == 0 )
    {
        printf ( "%lf [s] spent on computation\n", max_elapsed );
        write_image ( img );
        free ( img );
    }
    
    // Free the local image.
    free(local_img);
    
    MPI_Finalize ();
    exit ( EXIT_SUCCESS );
}