diff --git a/TP/Makefile b/TP/Makefile new file mode 100755 index 0000000000000000000000000000000000000000..ff1401dc21283232c42f6448b4dc95ce5ed48051 --- /dev/null +++ b/TP/Makefile @@ -0,0 +1,25 @@ +#Purge implicites rules +.SUFFIXES: +.SUFFIXES: .pdf .tex + +# Tools +TEX = pdflatex +BIBMAKE = bibtex +# name of Pdf File + + +tp2.pdf : tp2.tex + $(TEX) $< + $(TEX) $< +# $(TEX) $< + +# Clean the directory +clean::clean_aux + + +clean_pdf: + @echo Cleaning pdf File... + @rm -f tp2.pdf +clean_aux: + @echo Cleaning Auxiliary files... + @rm -f *.vrb *.out *.aux *.toc *.nav *.lof *.log *.snm diff --git a/TP/cmake_def.tex b/TP/cmake_def.tex new file mode 100644 index 0000000000000000000000000000000000000000..ebfd2e7f6f520aa03df2ceb8a3c4d8a3c528f243 --- /dev/null +++ b/TP/cmake_def.tex @@ -0,0 +1,11 @@ +\lstdefinelanguage{mycmake} + { + morekeywords={if,endif,foreach,endforeach,else,copy,file,install,add_executable,add_library, + target_link_libraries,endforeach,add_directories,project,while, endwhile, function, + enable_testing,add_test,enable_language, set, cmake_minimum_required, math, + find_package, find_program, find_library, string, list, elseif,include_directories, add_definitions, + set_tests_properties, message, include, add_subdirectory}, + sensitive=true, + morecomment=[l]{\#},% + morestring=[b]" + } diff --git a/TP/sources/session2/C/.gitignore b/TP/sources/session2/C/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..1287e62cb1e7b8c20f5ad33f413dba63b6f2cc25 --- /dev/null +++ b/TP/sources/session2/C/.gitignore @@ -0,0 +1,7 @@ +build +sol_* +*.o +heat +heat_par +heat_seq +heat.avi diff --git a/TP/sources/session2/C/CMakeLists.txt b/TP/sources/session2/C/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..247f59012cbebe314c82bd7547adbd0324a40a84 --- /dev/null +++ b/TP/sources/session2/C/CMakeLists.txt @@ -0,0 +1,9 @@ +cmake_minimum_required (VERSION 2.8) +project (PatcHeat) +enable_language(C) + +set(HEAT_SOURCES heat_seq.c heat.c mat_utils.c) +add_executable(patcHeatSeq ${HEAT_SOURCES} + heat.h mat_utils.h) + +target_link_libraries(patcHeatSeq m) diff --git a/TP/sources/session2/C/README b/TP/sources/session2/C/README new file mode 100644 index 0000000000000000000000000000000000000000..6ed03cc51f2dd5c3bc45257315a06fbbe08fa7ab --- /dev/null +++ b/TP/sources/session2/C/README @@ -0,0 +1,35 @@ +A heat propagation computation program + +This program can be compiled using CMake (minimum version 1.8). Simple +build can be done by 'mkdir build && cd build && cmake .. && make'. +Installing can then be performed with 'make install'. See http://cmake.org +for more details on CMake. + +To run it do "./heat nx ny iter_max save" where: + nx and ny are the number of discretisation points in X and Y + iter_max is the maximal number of iterations in temporal loop + save is a boolean flag (1 or 0) telling wether to save the state of + the matrix after each iteration + +Actually the program propagate on a (nx + 2) x (ny + 2) matrix the mean value +using a cross stencil at each iteration. So, for each x in [1..nx] and y +in [1..ny]: + +m_{i+1}(x,y) = wx * (m_i(x-1,y) + m_i(x+1,y)) + + wy * (m_i(x,y-1) + m_i(x,y+1)) + + d * m_i(x, y) +where wx = dt / hx ^ 2 and wy = dt / hy ^ 2 (hx and hy are the precision of the +derivation over x and y), d = 1 - 2 * wx - 2 * wy. + +So, if you call the software by "./heat 10 10 12 1" it will do 12 iterations +with a 10x10 matrix and produce 12 files "log_00001" to "log_00012" +containing the state of the matrix after each iteration. It should output +something like: +$ ./heat 10 10 12 1 +it = 0, t = 0.000e+00, err = 1.732e+00 +it = 10, t = 2.500e-02, err = 2.518e-01 + +To convert the resulting "log_XXXXX" files, you can use the heat_plot.py +script (simply run it with "python heat_plot.py") and it will generates +"XXXXX.png" files and a heat.avi video that you can watch. This python script +requires the FFMPEG library and the python library called matplotlib. diff --git a/TP/sources/session2/C/README.dox b/TP/sources/session2/C/README.dox new file mode 100644 index 0000000000000000000000000000000000000000..ef4f32f7c6debf3190e97572ea40d38510bb1b82 --- /dev/null +++ b/TP/sources/session2/C/README.dox @@ -0,0 +1,56 @@ +/*! +@mainpage A heat propagation computation program + +The heat propagation computation program computes the heat propagation equation +approximation using an iterative approach. This file decribes how the program +works. + +@section Installation + + +This program can be compiled using CMake (minimum version 1.8). Simple +build can be done by @code mkdir build && cd build && cmake .. && make @endcode. +Installing can then be performed with @code make install @endcode. See http://cmake.org +for more details on CMake. + +@section Usage + +To run it do @code ./heat nx ny iter_max save @endcode where: + @li @e nx and @e ny are the number of discretisation points in X and Y + @li @e iter_max is the maximal number of iterations in temporal loop + @li @e save is a boolean flag (1 or 0) telling wether to save the state of + the matrix after each iteration + +Actually the program propagates on a (@e nx + 2) x (@e ny + 2) matrix +the mean value using a cross stencil at each iteration. +So, + +\f$\forall x \in [1 \dots nx], y \in [1 \dots ny], \forall i \geq 0 + +m_{i+1}(x,y) = w_x (m_i(x-1,y) + m_i(x+1,y)) + + w_y (m_i(x,y-1) + m_i(x,y+1)) + + d m_i(x, y) + +\ where\ +w_x = \frac{dt}{h_x^2}, w_y = \frac{d_t}{hy^2}, d = 1 - 2 w_x - 2 w_y +\f$ + +where @e hx and @e hy are the precision of the derivation over @e x and @e y. + +So, if you call the software by @"./heat 10 10 12 1" it will do 12 iterations +with a 10x10 matrix and produce 12 files "log_00001" to "log_00012" +containing the state of the matrix after each iteration. It should output +something like: +@code +$ ./heat 10 10 12 1 +it = 0, t = 0.000e+00, err = 1.732e+00 +it = 10, t = 2.500e-02, err = 2.518e-01 +@endcode + +\section Conversion of the output to video + +To convert the resulting "log_XXXXX" files, you can use the heat_plot.py +script (simply run it with "python heat_plot.py") and it will generates +"XXXXX.png" files and a heat.avi video that you can watch. This python script +requires the FFMPEG library and the python library called matplotlib. +*/ diff --git a/TP/sources/session2/C/heat.c b/TP/sources/session2/C/heat.c new file mode 100644 index 0000000000000000000000000000000000000000..61cd5d5836aa06699dd2425fb8ecc8a8b1184f7d --- /dev/null +++ b/TP/sources/session2/C/heat.c @@ -0,0 +1,27 @@ +#include "heat.h" + +double +heat (double hx, double hy, double dt, + int size_x, int size_y, + const double *u_in, double *u_out) +{ + + int i, j; + double w_x, w_y, err, d; + + w_x = dt / (hx * hx); + w_y = dt / (hy * hy); + d = 1. - 2. * w_x - 2. * w_y; + + err = 0.; + for (i = 1; i < size_x - 1; ++i) + { + for (j = 1; j < size_y - 1; ++j) + { + u_out[i * size_y + j] = d * u_in[ i * size_y + j] + w_x * (u_in[(i - 1) * size_y + j] + u_in[(i + 1) * size_y + j]) + + w_y * (u_in[i * size_y + j - 1] + u_in[i * size_y + j + 1]); + err += SQR (u_out[i * size_y + j] - u_in[i * size_y + j]); + } + } + return err; +} diff --git a/TP/sources/session2/C/heat.h b/TP/sources/session2/C/heat.h new file mode 100644 index 0000000000000000000000000000000000000000..39c4a57cd12c1dd76312be6210258dce4b0f6c6e --- /dev/null +++ b/TP/sources/session2/C/heat.h @@ -0,0 +1,55 @@ +/** + * @file heat.h + * @author Inria SED Bordeaux + * @brief Heat computation iteration code + * + * @details This file declares the heat() function. + */ +#ifndef __HEAT_H +#define __HEAT_H + +/** + * Computes the square of @e X, e.g., 4 for @e X = 2, 9 for @e X = 3, etc. + * + * @param X the value to square + * @return the square of @e X, that is @e X * @e X + */ +#define SQR(X) ((X) * (X)) + +/** + * Computes the minimum between @e X and @e Y. E.g., + * <code>MIN(1,2) = MIN(2,1) = 1</code> + * + * @param X the first value to test + * @param Y the second value to test + * @return the minimum value between @e X and @e Y + */ +#define MIN(X,Y) (((X) <= (Y)) ? (X) : (Y)) + + +/** + * Do a single iteration of the heat equation iterative computation on + * a 2-D cartesian map. + * + * For each cell in @e u_in that is in <code>[1..size_x]x[1..size_y]</code>, + * it put in @e u_out the result of one iteration computation using a + * cross-stencil. + * The complete equation can be found in the @e README file of the project. + * + * @param hx precision of the derivation over x + * @param hy precision of the derivation over y + * @param dt precision of the derivation over time + * @param size_x the size of the cartesian map in x + * @param size_y the size of the cartesion map in y + * @param u_in the input map + * @param u_out the output map + * @return the square of the quadratic differences between @e u_in and @e u_out after + * the execution of the heat function + */ +double +heat (double hx, double hy, double dt, + int size_x, int size_y, + const double *u_in, double *u_out); + + +#endif diff --git a/TP/sources/session2/C/heat_par.c b/TP/sources/session2/C/heat_par.c new file mode 100644 index 0000000000000000000000000000000000000000..0acafadad8881f71f25ab5924eac0ce1f42c3908 --- /dev/null +++ b/TP/sources/session2/C/heat_par.c @@ -0,0 +1,305 @@ +/** + * @file heat_par.c + * @author Inria SED Bordeaux + * @brief Heat computation main procedure code for parallele computation + * + * @details This file declares the entry point (main() procedure) of the program. + * + */ +#include <stdio.h> +#include <math.h> +#include <stdlib.h> +#include <string.h> +#include <mpi.h> +#include "heat.h" +#include "mat_utils.h" + +/** + * Procedure which put ones on the boundaries of the global solution. + * + * If the current process touchs the boundary of the 2D cartesian topology + * we have to puts ones of the corresponding column or row. + * + * @param coo integer pair containing the coordinates in 2D cartesian topology + * @param nc_x number of processes in the x-dimension + * @param nc_y number of processes in the y-dimension + * @param size_x rows number of the local part of the solution + * @param size_y columns number of the local part of the solution + * @param u local part of the solution + * + * @internal + */ +static void +set_bounds (const int *coo, int nc_x, int nc_y, int size_x, int size_y, + double *u) +{ + + int i, j; + for (i = 0; i < size_x; ++i) + { + if (coo[1] == 0) + u[i * size_y + 0] = 1.; + if (coo[1] == (nc_y - 1)) + u[i * size_y + size_y - 1] = 1.; + } + + for (j = 0; j < size_y; ++j) + { + if (coo[0] == 0) + u[0 * size_y + j] = 1.; + if (coo[0] == (nc_x - 1)) + u[(size_x - 1) * size_y + j] = 1.; + } + +} + +/** + * Procedure for swapping the boundaries (2 columns and 2 rows) + * with the neighbour processes + * + * @param comm the 2D communicator + * @param col the column type + * @param neighbours the set of the local neighbours processes + * @param size_x rows number of the local part of the solution + * @param size_y columns number of the local part of the solution + * @param u local part of the solution + * + * @internal + */ +static void +ghosts_swap (MPI_Comm comm, MPI_Datatype col, + const int *neighbours, int size_x, + int size_y, double *u) +{ + + int N = 0, S = 1, E = 2, W = 3; + int s_tag, r_tag; + MPI_Status status; + + /* N --> S + N block last significant row goes to S block first ghost row */ + s_tag = 0; + r_tag = 0; + MPI_Sendrecv (&u[(size_x - 2) * size_y + 0], size_y, MPI_DOUBLE, neighbours[S], + s_tag, &u[+0 * size_y + 0], size_y, MPI_DOUBLE, neighbours[N], + r_tag, comm, &status); + /* S --> N + S block first significant row goes to N block last ghost row */ + s_tag = 1; + r_tag = 1; + MPI_Sendrecv (&u[1 * size_y + 0], size_y, MPI_DOUBLE, neighbours[N], s_tag, + &u[(size_x - 1) * size_y + 0], size_y, MPI_DOUBLE, neighbours[S], + r_tag, comm, &status); + + /* W --> E + W block last significant column goes to E block first ghost column */ + s_tag = 2; + r_tag = 2; + MPI_Sendrecv (&u[1 * size_y + size_y - 2], 1, col, neighbours[E], s_tag, + &u[1 * size_y + 0], 1, col, neighbours[W], r_tag, comm, &status); + + /* E --> W + E block first significant column goes to W block last ghost column */ + s_tag = 3; + r_tag = 3; + MPI_Sendrecv (&u[1 * size_y + 1], 1, col, neighbours[W], s_tag, + &u[1 * size_y + size_y - 1], 1, col, neighbours[E], r_tag, comm, + &status); + +} + + +/** + * A usage function + * + * This function prints out the normal usage of the program and exit the program + * with failure. + * + * @param argv the array of arguments passed to the main procedure + * @internal + */ +static void +usage(char *argv[]) +{ + fprintf(stderr, "Usage: mpirun -np (px*py) %s nx ny iter_max px py save\n", argv[0]); + fprintf(stderr, "\tnx number of discretisation points in X\n"); + fprintf(stderr, "\tny number of discretisation points in Y\n"); + fprintf(stderr, "\titer_max maximal number of iterations in temporal loop\n"); + fprintf(stderr, "\tpx X process number\n"); + fprintf(stderr, "\tpy Y process number\n"); + fprintf(stderr, "\tsave boolean flag (1 or 0) for recording states\n"); + exit(EXIT_FAILURE); +} + +/** + * Main procedure + * + * Of course, this is the entry point :P + * + * @param argc the number of program arguments + * @param argv the list of program arguments + * @return the error code of the program + * @internal + */ +int +main (int argc, char *argv[]) +{ + + int nx, ny, i, j, r, size_x, size_y, cell_x, cell_y, nc_x, nc_y; + + double hx, hy, dt, err_loc, err, iter_max, prec; + + double *u_in, *u_out, *vec_temp, *solution; + + int rank_w, size_w; + + int rank_2D; + MPI_Comm comm2D; + MPI_Datatype type_col; + int ndims = 2; + int save; + int dims[2], periods[2], coords[2], reorder; + + int N = 0, S = 1, E = 2, W = 3; + int neighbours[4]; + + + MPI_Init (&argc, &argv); + MPI_Comm_rank (MPI_COMM_WORLD, &rank_w); + MPI_Comm_size (MPI_COMM_WORLD, &size_w); + + + if (argc < 7) + { + usage(argv); + } + + + nx = atoi (argv[1]); + ny = atoi (argv[2]); + iter_max = atoi (argv[3]); + nc_x = atoi (argv[4]); + nc_y = atoi (argv[5]); + save = atoi (argv[6]); + + if (nc_x * nc_y != size_w) + { + printf + (" the total number of processus differs from the product px * py : %d x %d != %d \n", + nc_x, nc_y, size_w); + exit (-1); + } + dims[0] = nc_x; + dims[1] = nc_y; + periods[0] = 0; + periods[1] = 0; + reorder = 1; + + // creation of a 2D cartesian topology (nc_x x nc_y processus) + MPI_Cart_create (MPI_COMM_WORLD, ndims, dims, periods, reorder, &comm2D); + // fetch the rank of the topology + MPI_Comm_rank (comm2D, &rank_2D); + // fetch the coords of the current process + MPI_Cart_coords (comm2D, rank_2D, ndims, coords); + + // we identify neighbors in the cartesian topology + // neighbours[point] contains the rank in the 2-D communicator + // of the cell corresponding to point + MPI_Cart_shift (comm2D, 0, 1, &neighbours[N], &neighbours[S]); + MPI_Cart_shift (comm2D, 1, 1, &neighbours[W], &neighbours[E]); + + + cell_x = nx / nc_x; + cell_y = ny / nc_y; + + size_x = cell_x + 2; + size_y = cell_y + 2; + + + // creation of a MPI column type to transfer a column (non contiguous in C) + // from one cell to another + MPI_Type_vector (cell_x, 1, size_y, MPI_DOUBLE, &type_col); + MPI_Type_commit (&type_col); + + + u_in = (double *) calloc (size_x * size_y, sizeof (double)); + if (u_in == NULL) + { + printf("not enough memory!\n"); + exit(-1); + } + + u_out = (double *) calloc (size_x * size_y, sizeof (double)); + if (u_out == NULL) + { + printf("not enough memory!\n"); + exit(-1); + } + + set_bounds (coords, nc_x, nc_y, size_x, size_y, u_in); + set_bounds (coords, nc_x, nc_y, size_x, size_y, u_out); + + + hx = 1. / nx; + hy = 1. / ny; + dt = MIN (SQR (hx) / 4., SQR (hy) / 4.); + prec = 1e-4; + err = 1e10; + // temporal loop + for (i = 0; i < iter_max; ++i) + { + + err_loc = heat (hx, hy, dt, size_x, size_y, u_in, u_out); + + // retrieve local error to compute the global error + MPI_Allreduce (&err_loc, &err, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + err = sqrt (err); + if (rank_w == 0 && (i % 10 == 0)) + printf ("it = %d, t = %.3e, err = %.3e\n", i, i * dt, err); + memcpy (u_in, u_out, sizeof (double) * size_x * size_y); + + ghosts_swap (comm2D, type_col, neighbours, size_x, size_y, u_in); + + if (err <= prec) + break; + } + + vec_temp = (double *) calloc (nc_x * nc_y * size_x * size_y, sizeof (double)); + // gather the big solution on the process of rank 0 + MPI_Gather (u_in, size_x * size_y, MPI_DOUBLE, vec_temp, size_x * size_y, + MPI_DOUBLE, 0, MPI_COMM_WORLD); + + if (rank_w == 0) + { + solution = (double *) calloc (nx * ny, sizeof (double)); + // loop over coordinates to re-organize the big flat vector + for (r = 0; r < size_w; ++r) + { + int coo2[2]; + MPI_Cart_coords (comm2D, r, ndims, coo2); + // we copy the flattened submatrix with coords (i,j) + // in the global solution matrix + for (i = 1; i < size_x - 1; ++i) + { + for (j = 1; j < size_y - 1; ++j) + { + solution[(coo2[0] * (size_x - 2) + (i - 1)) * ny + + coo2[1] * (size_y - 2) + j - 1] = + vec_temp[(coo2[0] * nc_y + coo2[1]) * size_x * size_y + + i * size_y + j]; + } + } + if (save) + save_mat ("sol_para.txt", nx, ny, solution); + } + free (solution); + } + + free (vec_temp); + MPI_Type_free (&type_col); + free (u_in); + free (u_out); + + MPI_Finalize (); + return 0; +} diff --git a/TP/sources/session2/C/heat_plot.py b/TP/sources/session2/C/heat_plot.py new file mode 100644 index 0000000000000000000000000000000000000000..1b8cea9a7c3d7cad61145a596a4fb8621eff73bc --- /dev/null +++ b/TP/sources/session2/C/heat_plot.py @@ -0,0 +1,9 @@ +import os +import numpy as np +import matplotlib +import matplotlib.pyplot as plt +sols= np.array([ np.loadtxt(f) for f in os.listdir('.') if 'sol' in f]) +for i in range(1,sols.shape[0]): + plt.imsave('%05d.png' % i, sols[i,:,:],cmap = plt.cm.jet,vmin=0,vmax=1) +comm="ffmpeg -i %05d.png heat.avi" +os.system(comm) diff --git a/TP/sources/session2/C/heat_seq.c b/TP/sources/session2/C/heat_seq.c new file mode 100644 index 0000000000000000000000000000000000000000..65c62785c7201f413f6453c66bc2f6fdfbe1858d --- /dev/null +++ b/TP/sources/session2/C/heat_seq.c @@ -0,0 +1,186 @@ +/** + * @file heat_seq.c + * @author Inria SED Bordeaux + * @brief Heat computation main procedure code + * + * @details This file declares the entry point (main() procedure) of the program. + * + */ +#include <stdio.h> +#include <math.h> +#include <stdlib.h> +#include <string.h> +#include "heat.h" +#include "mat_utils.h" + +/** + * Allocates code with error handling + * + * Allocates a pointer @e var to @e s object of type @e type. This + * macros exit with a proper message on failure of the allocation. + * For example: + * @code + * int *buf; + * xalloc(buf, 10, int); + * @endcode + * allocates an array of 10 integer in buf. + * + * @param var the pointer to allocate + * @param s the number of elements to allocate + * @param type the type of elements to allocate + * + * @internal + */ +#define xalloc(var, s, type) do { \ + var = (type *) calloc(s, sizeof(type)); \ + if(var == NULL) { \ + perror("calloc"); \ + exit(EXIT_FAILURE); \ + } \ + } while(0) + +/** + * Set the boundaries of a 2-D map to 1.0 + * + * This function will set to 1.0 all the cells on the boundaries of a 2-D + * double maps @e u of size @e size_x in X and @e size_y in Y. @e u should + * be an array of size (@e size_x + 2) * (@e size_y + 2). + * + * For example, if @e u = [0, 0, 0, 0, 0, 0, 0, 0, 0] then after a call to + * set_bounds(1, 1, u), @e u will contains [1, 1, 1, 1, 0, 1, 1, 1, 1]. In + * the same way, if @e u = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + * then after a call to set_bounds(2, 2, u), @e u will contains + * [1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1]. That is in 2-D: + * @code + * 1 1 1 1 + * 1 0 0 1 + * 1 0 0 1 + * 1 1 1 1 + * @endcode + * + * @param size_x The size in X of the @e u map + * @param size_y The size in Y of the @e u map + * @param u The map to set boundaries to 1.0 + * + * @internal + */ +static void +set_bounds(int size_x, int size_y, double *u) +{ + int i, j; + + for (i = 0; i < size_x; ++i) + { + u[i * size_y + 0] = 1.; + u[i * size_y + size_y -1] = 1.; + } + + for (j = 0; j < size_y; ++j) + { + u[0 * size_y + j] = 1.; + u[(size_x - 1) * size_y + j] = 1.; + } +} + +/** + * Compute the heat propagation equation using iterative approach + * until convergence. + * + * This function will compute the heat propagation equation on the map + * @e u_in using the iterative approach outputing the result in @e u_out. + * @e dt, @e hx and @e hy give the derivation approximation steps over time, + * X and Y. @e iter_max limits the maximum number of iteration if the computation + * does not converge. + * + * @param hx the derivation approximation step in X + * @param hy the derivation approximation step in Y + * @param dt the derivation approximation step in time + * @param iter_max the maximum number of iteration to perform. + * @param save a boolean indicating if the results is to be saved in a file after + * each step + * @param size_x The size in X of the @e u_in and @e u_out maps + * @param size_y The size in Y of the @e u_in and @e u_out maps + * @param u_in the input map, it will be modified and invalidated by the + * call of this function + * @param u_out the output map, it will contains the result of the computation + * + * @internal + */ +static void +compute_heat_propagation(double hx, double hy, double dt, + int iter_max, int save, + int size_x, int size_y, + double *u_in, double *u_out) +{ + int i; + char name_fic[120]; + double err, prec; + + prec = 1e-7; + err = 1e10; + for (i = 0; i < iter_max; ++i) + { + err = heat (hx, hy, dt, size_x, size_y, u_in, u_out); + err = sqrt (err); + if (i % 10 == 0) + { + printf ("it = %d, t = %.3e, err = %.3e\n", i, i * dt, err); + } + if (save) + { + sprintf (name_fic,"sol_%05d", i); + save_mat (name_fic, size_x, size_y, u_in); + } + memcpy (u_in, u_out, sizeof (double) * size_x * size_y); + if (err <= prec) + break; + } +} + +/** + * Main procedure + * + * Of course, this is the entry point :P + * + * @param argc the number of program arguments + * @param argv the list of program arguments + * @return the error code of the program + * @internal + */ +int +main (int argc, char *argv[]) +{ + int nx, ny, size_x, size_y, save, iter_max; + double hx, hy, dt; + double *u_in, *u_out; + + if (argc < 5) + { + return EXIT_FAILURE; + } + + nx = atoi (argv[1]); + ny = atoi (argv[2]); + iter_max = atoi (argv[3]); + save = atoi(argv[4]); + + hx = 1. / nx; + hy = 1. / ny; + dt = MIN (SQR (hx) / 4., SQR (hy) / 4.); + + size_x = nx + 2; + size_y = ny + 2; + + xalloc (u_in, size_x * size_y, double); + xalloc (u_out, size_x * size_y, double); + + set_bounds (size_x, size_y, u_in); + set_bounds (size_x, size_y, u_out); + + compute_heat_propagation (hx, hy, dt, iter_max, save, size_x, size_y, + u_in, u_out); + + free (u_in); + + return EXIT_SUCCESS; +} diff --git a/TP/sources/session2/C/mat_utils.c b/TP/sources/session2/C/mat_utils.c new file mode 100644 index 0000000000000000000000000000000000000000..e0c670c82a8905a47bff676fea8438914795329f --- /dev/null +++ b/TP/sources/session2/C/mat_utils.c @@ -0,0 +1,51 @@ +#include <stdio.h> +#include "mat_utils.h" + +void +print_mat (int size_x, int size_y, const double *u) +{ + + int i, j; + printf ("\n["); + for (i = 0; i < size_x - 1; ++i) + { + for (j = 0; j < size_y - 1; ++j) + { + printf ("%f,", u[i * size_y + j]); + } + printf ("%f\n", u[i * size_y + size_y - 1]); + } + for (j = 0; j < size_y - 1; ++j) + { + printf ("%f,", u[(size_x - 1) * size_y + j]); + } + printf ("%f]\n", u[(size_x - 1) * size_y + size_y - 1]); +} + +void +save_mat (const char *filename, int size_x, int size_y, + const double *u) +{ + int i, j; + FILE *fid; + + if (filename == NULL) { + fprintf(stderr, "Error: input argument <filename> NULL\n"); + return; + } + fid = fopen (filename, "w"); + if (fid == NULL) { + fprintf(stderr, "Error: unable to open <%s>\n", filename); + return; + } + for (i = 0; i < size_y; ++i) + { + for (j = 0; j < size_y; ++j) + { + fprintf (fid, "%.15e ", u[i * size_y + j]); + } + fprintf (fid, "\n"); + } + fclose (fid); +} + diff --git a/TP/sources/session2/C/mat_utils.h b/TP/sources/session2/C/mat_utils.h new file mode 100644 index 0000000000000000000000000000000000000000..45a61eb0082814983d379600d106352f4c96a790 --- /dev/null +++ b/TP/sources/session2/C/mat_utils.h @@ -0,0 +1,28 @@ +#ifndef __MAT_UTILS_H +#define __MAT_UTILS_H + +void +print_mat (int size_x, int size_y, const double *u); + +/** + * A matrix saving function + * + * It saves in a file the matrix given in the parameter @e u. + * For instance if <code>u = [[1., 2.], [3., 4.]]</code>, then + * save_mat("toto",2,2,u) will output in the file toto: + * @code + * 1.000000000000000e+00 2.000000000000000e+00 + * 3.000000000000000e+00 4.000000000000000e+00 + * @endcode + * + * @param filename the file to output the result + * @param size_x the size in X of the matrix u + * @param size_y the size in Y of the matrix u + * @param u the matrix to print + */ +void +save_mat (const char *filename, int size_x, int size_y, + const double *u); + +#endif + diff --git a/TP/sources/session2/C/test.sh b/TP/sources/session2/C/test.sh new file mode 100755 index 0000000000000000000000000000000000000000..636adbc19f1a6b78e85985430c775cdf52b78558 --- /dev/null +++ b/TP/sources/session2/C/test.sh @@ -0,0 +1 @@ +echo $0 $* diff --git a/TP/sources/session2/Fortran/.gitignore b/TP/sources/session2/Fortran/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..1287e62cb1e7b8c20f5ad33f413dba63b6f2cc25 --- /dev/null +++ b/TP/sources/session2/Fortran/.gitignore @@ -0,0 +1,7 @@ +build +sol_* +*.o +heat +heat_par +heat_seq +heat.avi diff --git a/TP/sources/session2/Fortran/CMakeLists.txt b/TP/sources/session2/Fortran/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4ad0802bc32d5b5520858870b63ccd9f95a78b41 --- /dev/null +++ b/TP/sources/session2/Fortran/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required (VERSION 2.8) +project (PatcHeat) +enable_language(Fortran) + +# add the subdirectory +add_subdirectory(lib) +# add the executable +set(HEAT_SOURCES heat_seq.f90 mat_utils.f90) +add_executable(patcHeatSeq ${HEAT_SOURCES}) +target_link_libraries(patcHeatSeq Heat) + +find_package(MPI) +if (MPI_FOUND) + add_executable(parallel heat_par.f90 mat_utils.f90) + target_link_libraries(parallel ${MPI_Fortran_LIBRARIES} Heat) +else () +endif() # MPI_FOUND +enable_testing() +message(STATUS "bin = ${CMAKE_BINARY_DIR}") +add_test(seq10x10 ${CMAKE_BINARY_DIR}/patcHeatSeq 10 10 1000 0) +add_test(seq100x100 ${CMAKE_BINARY_DIR}/patcHeatSeq 100 100 1000 0) diff --git a/TP/sources/session2/Fortran/README b/TP/sources/session2/Fortran/README new file mode 100644 index 0000000000000000000000000000000000000000..6ed03cc51f2dd5c3bc45257315a06fbbe08fa7ab --- /dev/null +++ b/TP/sources/session2/Fortran/README @@ -0,0 +1,35 @@ +A heat propagation computation program + +This program can be compiled using CMake (minimum version 1.8). Simple +build can be done by 'mkdir build && cd build && cmake .. && make'. +Installing can then be performed with 'make install'. See http://cmake.org +for more details on CMake. + +To run it do "./heat nx ny iter_max save" where: + nx and ny are the number of discretisation points in X and Y + iter_max is the maximal number of iterations in temporal loop + save is a boolean flag (1 or 0) telling wether to save the state of + the matrix after each iteration + +Actually the program propagate on a (nx + 2) x (ny + 2) matrix the mean value +using a cross stencil at each iteration. So, for each x in [1..nx] and y +in [1..ny]: + +m_{i+1}(x,y) = wx * (m_i(x-1,y) + m_i(x+1,y)) + + wy * (m_i(x,y-1) + m_i(x,y+1)) + + d * m_i(x, y) +where wx = dt / hx ^ 2 and wy = dt / hy ^ 2 (hx and hy are the precision of the +derivation over x and y), d = 1 - 2 * wx - 2 * wy. + +So, if you call the software by "./heat 10 10 12 1" it will do 12 iterations +with a 10x10 matrix and produce 12 files "log_00001" to "log_00012" +containing the state of the matrix after each iteration. It should output +something like: +$ ./heat 10 10 12 1 +it = 0, t = 0.000e+00, err = 1.732e+00 +it = 10, t = 2.500e-02, err = 2.518e-01 + +To convert the resulting "log_XXXXX" files, you can use the heat_plot.py +script (simply run it with "python heat_plot.py") and it will generates +"XXXXX.png" files and a heat.avi video that you can watch. This python script +requires the FFMPEG library and the python library called matplotlib. diff --git a/TP/sources/session2/Fortran/README.dox b/TP/sources/session2/Fortran/README.dox new file mode 100644 index 0000000000000000000000000000000000000000..ef4f32f7c6debf3190e97572ea40d38510bb1b82 --- /dev/null +++ b/TP/sources/session2/Fortran/README.dox @@ -0,0 +1,56 @@ +/*! +@mainpage A heat propagation computation program + +The heat propagation computation program computes the heat propagation equation +approximation using an iterative approach. This file decribes how the program +works. + +@section Installation + + +This program can be compiled using CMake (minimum version 1.8). Simple +build can be done by @code mkdir build && cd build && cmake .. && make @endcode. +Installing can then be performed with @code make install @endcode. See http://cmake.org +for more details on CMake. + +@section Usage + +To run it do @code ./heat nx ny iter_max save @endcode where: + @li @e nx and @e ny are the number of discretisation points in X and Y + @li @e iter_max is the maximal number of iterations in temporal loop + @li @e save is a boolean flag (1 or 0) telling wether to save the state of + the matrix after each iteration + +Actually the program propagates on a (@e nx + 2) x (@e ny + 2) matrix +the mean value using a cross stencil at each iteration. +So, + +\f$\forall x \in [1 \dots nx], y \in [1 \dots ny], \forall i \geq 0 + +m_{i+1}(x,y) = w_x (m_i(x-1,y) + m_i(x+1,y)) + + w_y (m_i(x,y-1) + m_i(x,y+1)) + + d m_i(x, y) + +\ where\ +w_x = \frac{dt}{h_x^2}, w_y = \frac{d_t}{hy^2}, d = 1 - 2 w_x - 2 w_y +\f$ + +where @e hx and @e hy are the precision of the derivation over @e x and @e y. + +So, if you call the software by @"./heat 10 10 12 1" it will do 12 iterations +with a 10x10 matrix and produce 12 files "log_00001" to "log_00012" +containing the state of the matrix after each iteration. It should output +something like: +@code +$ ./heat 10 10 12 1 +it = 0, t = 0.000e+00, err = 1.732e+00 +it = 10, t = 2.500e-02, err = 2.518e-01 +@endcode + +\section Conversion of the output to video + +To convert the resulting "log_XXXXX" files, you can use the heat_plot.py +script (simply run it with "python heat_plot.py") and it will generates +"XXXXX.png" files and a heat.avi video that you can watch. This python script +requires the FFMPEG library and the python library called matplotlib. +*/ diff --git a/TP/sources/session2/Fortran/heat_par.f90 b/TP/sources/session2/Fortran/heat_par.f90 new file mode 100644 index 0000000000000000000000000000000000000000..0c648a272119c18254443800a845e277fd76d201 --- /dev/null +++ b/TP/sources/session2/Fortran/heat_par.f90 @@ -0,0 +1,232 @@ +!------------------------------------------------------------------------------- +!> @file heat_par.f90 +!> @author Inria SED Bordeaux +!> @brief Heat computation main procedure code for parallel computation. +!> +!> @details This file declares the entry point (main() procedure) of the program. +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +!> Main procedure +!> +!> Of course, this is the entry point :P +!> +!> @internal +!------------------------------------------------------------------------------- +program heat_par + implicit none + include 'mpif.h' + integer :: nx, ny, nc_x, nc_y, i, r, cell_x, cell_y, size_x, size_y, iter_max, ierr + integer, parameter :: ndims = 2 + integer :: rank_w, size_w, rank_2D, comm2D, type_row , alloc_status, narg + integer :: N=1, S = 2, E = 3, W = 4 + logical :: verbose, reorder = .true. + integer, dimension(4) :: neighbour + double precision hx, hy, dt, error, error_loc, prec + double precision, allocatable :: u_in(:,:), u_out(:,:), solution(:,:), vec_temp(:), u_vec(:) + character(len=10)::param + + integer, dimension(ndims) :: dims , coords, coo + logical, dimension(ndims) :: periods = .false. + + + call MPI_Init(ierr) + call MPI_Comm_rank(MPI_COMM_WORLD, rank_w, ierr) + call MPI_Comm_size(MPI_COMM_WORLD, size_w, ierr) + + verbose = (rank_w == 0) + + narg = command_argument_count() + if (narg < 5) then + call usage(); + end if + + call get_command_argument(1, param) + read(param,*) nx + + call get_command_argument(2, param) + read(param,*) ny + + call get_command_argument(3, param) + read(param,*) iter_max + + call get_command_argument(4, param) + read(param,*) nc_x + + call get_command_argument(5, param) + read(param,*) nc_y + + + if (nc_x * nc_y /= size_w) then + if (verbose) print *, 'the total number of processus differs from the product px * py ', & + & nc_x, ' x ' , nc_y , '/= ', size_w + call MPI_finalize(ierr) + stop + endif + + cell_x = nx / nc_x + cell_y = ny / nc_y + + size_x = cell_x + 2 + size_y = cell_y + 2 + + + hx = 1.d0/ nx + hy = 1.d0/ ny + dt = min(hx*hx/4.d0, hy*hy/4.d0) + + ! construction of the cartesion topology + dims(1) = nc_x + dims(2) = nc_y + + call MPI_Cart_create(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm2D, ierr) + call MPI_Comm_rank(comm2D, rank_2D, ierr) + + + ! Fetch the processus coordinates in the 2-D grid + call MPI_Cart_coords(comm2D, rank_2D, ndims, coords, ierr) + + ! Creation of a non-contiguous in memory column type + ! to address Fortran storage: no stride of size_x + call MPI_Type_vector(cell_y, 1, size_x, MPI_DOUBLE_PRECISION, type_row, ierr) + call MPI_Type_commit(type_row, ierr) + + ! fetching the neighbor ranks + call MPI_Cart_shift(comm2D, 0, 1, neighbour(N), neighbour(S), ierr) + call MPI_Cart_shift(comm2D, 1, 1, neighbour(W), neighbour(E), ierr) + + allocate(u_in(1:size_x,1:size_y), stat=alloc_status) + if (alloc_status /=0) stop "Not enough memory" + allocate(u_out(1:size_x,1:size_y), stat=alloc_status) + if (alloc_status /=0) stop "Not enough memory" + + call set_bounds(coords, nc_x, nc_y, size_x, size_y, u_in) + call set_bounds(coords, nc_x, nc_y, size_x, size_y, u_out) + + + prec = 1e-4 + error = 1e10 + do i=0, iter_max + + call heat(hx, hy, dt, size_x, size_y, u_in, u_out, error_loc) + call MPI_Allreduce(error_loc, error, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr) + error = sqrt(error) + + if (verbose .and. mod(i,10) == 0) print * , 'it =', i, 't = ', i*dt, 'err = ', error + + u_in = u_out + call ghosts_swap(comm2D, type_row, neighbour, size_x, size_y, u_in) + if (error <= prec) exit + + end do + + ! We gather the solution on process 0 + + if (verbose) then + allocate(vec_temp(1:nx * ny)) + end if + allocate(u_vec(1:cell_x * cell_y)) + u_vec = reshape(u_in(2:size_x - 1, 2:size_y - 1),(/cell_x * cell_y/)) + call MPI_Gather(u_vec , cell_x * cell_y , MPI_DOUBLE_PRECISION, & + vec_temp, cell_x * cell_y, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + if (verbose) then + allocate(solution(1:nx,1:ny)) + solution = reshape(vec_temp, (/nx,ny/)) + do r=1,size_w + call MPI_Cart_coords(comm2D, r-1, ndims, coo, ierr) + solution(coo(1) * cell_x + 1:(coo(1)+1) * cell_x, & + coo(2) * cell_y + 1:(coo(2)+1) * cell_y) = & + reshape(vec_temp(cell_x * cell_y * (r - 1) + 1: cell_x * cell_y * r),(/cell_x,cell_y/)) + end do + !do i=1,nx + ! print *, (real(solution(i,j)), j=1,ny) + !enddo + deallocate(vec_temp) + deallocate(solution) + end if + + deallocate(u_vec) + deallocate(u_in) + deallocate(u_out) + + call MPI_Type_free(type_row, ierr) + call MPI_Finalize(ierr) + +contains + + +!------------------------------------------------------------------------------- +!> A usage function +!> +!> This function prints out the normal usage of the program and exit the program +!> with failure. +!> +!> @internal +!------------------------------------------------------------------------------- +subroutine usage() + implicit none + character(len=255)::name + + call get_command_argument(0,name) + print *, 'Usage: mpirun -np (px*py) ', TRIM(name), ' nx ny iter_max px py' + print *, ' nx number of discretisation points in X' + print *, ' ny number of discretisation points in Y' + print *, ' iter_max maximal number of iterations in temporal loop' + print *, ' px X process number' + print *, ' py Y process number' + stop +end subroutine usage + + +subroutine ghosts_swap(comm, type_row, neighbour, size_x, size_y, u) + + integer, intent(in) :: size_x, size_y, comm, type_row + integer, dimension(4), intent(in) :: neighbour + double precision, dimension(1:size_x, 1:size_y), intent(inout) :: u + integer, parameter :: N = 1, S = 2, E = 3, W = 4 + integer :: ierr, s_tag, r_tag + integer, dimension(MPI_STATUS_SIZE) :: stat + + ! N --> S + ! N block last significant row goes to S block first ghost row + s_tag =0; r_tag = 0 + call MPI_Sendrecv(u(size_x - 1, 2), 1, type_row, neighbour(S), s_tag, & + & u(1, 2) , 1, type_row, neighbour(N), r_tag, comm, stat, ierr) + + ! S --> N + ! S block first significant row goes to N block last ghost row + s_tag =1; r_tag = 1 + call MPI_Sendrecv(u(2, 2), 1, type_row, neighbour(N), s_tag, & + & u(size_x , 2), 1, type_row, neighbour(S), r_tag, comm, stat, ierr) + + ! W --> E + ! W block last significant column goes to E block first ghost column + s_tag =2; r_tag = 2 + call MPI_Sendrecv(u(1, size_y - 1), size_x, MPI_DOUBLE_PRECISION, neighbour(E), s_tag,& + & u(1, 1) , size_x, MPI_DOUBLE_PRECISION, neighbour(W), r_tag, comm, stat, ierr) + + ! E --> W + ! E block first significant column goes to W block last ghost column + s_tag =3; r_tag = 3 + call MPI_Sendrecv(u(1, 2), size_x , MPI_DOUBLE_PRECISION, neighbour(W), s_tag, & + & u(1, size_y) , size_x, MPI_DOUBLE_PRECISION, neighbour(E), r_tag, comm, stat, ierr) + +end subroutine ghosts_swap + +subroutine set_bounds(coo, nc_x, nc_y, size_x, size_y, u) + + implicit none + integer, intent(in) :: size_x, size_y, nc_x, nc_y + integer, dimension(2), intent(in) :: coo + double precision, dimension(1:size_x, 1:size_y), intent(out) :: u + u = 0.d0 + if (coo(1) == 0) u(1, :) = 1.d0 + if (coo(1) == (nc_x - 1)) u(size_x, :) = 1.d0 + + if (coo(2) == 0) u(:, 1) = 1.d0 + if (coo(2) == (nc_y - 1)) u(:, size_y) = 1.d0 + +end subroutine set_bounds + +end program diff --git a/TP/sources/session2/Fortran/heat_plot.py b/TP/sources/session2/Fortran/heat_plot.py new file mode 100644 index 0000000000000000000000000000000000000000..1b8cea9a7c3d7cad61145a596a4fb8621eff73bc --- /dev/null +++ b/TP/sources/session2/Fortran/heat_plot.py @@ -0,0 +1,9 @@ +import os +import numpy as np +import matplotlib +import matplotlib.pyplot as plt +sols= np.array([ np.loadtxt(f) for f in os.listdir('.') if 'sol' in f]) +for i in range(1,sols.shape[0]): + plt.imsave('%05d.png' % i, sols[i,:,:],cmap = plt.cm.jet,vmin=0,vmax=1) +comm="ffmpeg -i %05d.png heat.avi" +os.system(comm) diff --git a/TP/sources/session2/Fortran/heat_seq.f90 b/TP/sources/session2/Fortran/heat_seq.f90 new file mode 100644 index 0000000000000000000000000000000000000000..0d69244516b8b37ad95932dc083e55afe17c2bac --- /dev/null +++ b/TP/sources/session2/Fortran/heat_seq.f90 @@ -0,0 +1,153 @@ +!------------------------------------------------------------------------------- +!> @file heat_seq.f90 +!> @author Inria SED Bordeaux +!> @brief Heat computation main procedure code +!> +!> @details This file declares the entry point (main() procedure) of the program. +!------------------------------------------------------------------------------- + + + +!------------------------------------------------------------------------------- +!> Main procedure +!> +!> Of course, this is the entry point :P +!> +!> @internal +!------------------------------------------------------------------------------- +program heat_seq + implicit none + integer :: nx, ny, size_x, size_y, iter_max, narg, save_flag + double precision :: hx, hy, dt + double precision, allocatable :: u_in(:,:), u_out(:,:) + character(len=10)::param + character(len=10)::name_save + + narg = command_argument_count() + if (narg < 4) then + stop + end if + + call get_command_argument(1,param) + read(param,*) nx + + call get_command_argument(2,param) + read(param,*) ny + + call get_command_argument(3,param) + read(param,*) iter_max + + call get_command_argument(4,param) + read(param,*) save_flag + + hx = 1./ nx; + hy = 1./ ny + dt = min(hx*hx/4., hy*hy/4.) + + size_x = nx + 2 + size_y = ny + 2 + + allocate(u_in(1:size_x,1:size_y)) + allocate(u_out(1:size_x,1:size_y)) + + call set_bounds(size_x, size_y, u_in) + call set_bounds(size_x, size_y, u_out) + + + call compute_heat_propagation(hx, hy, dt, iter_max, save_flag, size_x, size_y, & + u_in, u_out) + + deallocate(u_in) + +contains + +!------------------------------------------------------------------------------- +!> Set the boundaries of a 2-D map to 1.0 +!> +!> This function will set to 1.0 all the cells on the boundaries of a 2-D +!> double maps @e u of size @e size_x in X and @e size_y in Y. @e u should +!> be an array of size (@e size_x + 2) * (@e size_y + 2). +!> +!> For example, if @e u = [0, 0, 0, 0, 0, 0, 0, 0, 0] then after a call to +!> set_bounds(1, 1, u), @e u will contains [1, 1, 1, 1, 0, 1, 1, 1, 1]. In +!> the same way, if @e u = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +!> then after a call to set_bounds(2, 2, u), @e u will contains +!> [1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1]. That is in 2-D: +!> @code +!> 1 1 1 1 +!> 1 0 0 1 +!> 1 0 0 1 +!> 1 1 1 1 +!> @endcode +!> +!> @param size_x The size in X of the @e u map +!> @param size_y The size in Y of the @e u map +!> @param u The map to set boundaries to 1.0 +!> +!> @internal +!------------------------------------------------------------------------------- +subroutine set_bounds(size_x, size_y, u) + implicit none + integer, intent(in) :: size_x, size_y + double precision, dimension(1:size_x, 1:size_y), intent(inout) :: u + + u(1, :) = 1.d0 + u(size_x, :) = 1.d0 + + u(:, 1) = 1.d0 + u(:, size_y) = 1.d0 + +end subroutine set_bounds + + +!------------------------------------------------------------------------------- +!> Compute the heat propagation equation using iterative approach +!> until convergence. +!> +!> This function will compute the heat propagation equation on the map +!> @e u_in using the iterative approach outputing the result in @e u_out. +!> @e dt, @e hx and @e hy give the derivation approximation steps over time, +!> X and Y. @e iter_max limits the maximum number of iteration if the computation +!> does not converge. +!> +!> @param hx the derivation approximation step in X +!> @param hy the derivation approximation step in Y +!> @param dt the derivation approximation step in time +!> @param iter_max the maximum number of iteration to perform. +!> @param save a boolean indicating if the results is to be saved in a file after +!> each step +!> @param size_x The size in X of the @e u_in and @e u_out maps +!> @param size_y The size in Y of the @e u_in and @e u_out maps +!> @param u_in the input map, it will be modified and invalidated by the +!> call of this function +!> @param u_out the output map, it will contains the result of the computation +!> +!> @internal +!------------------------------------------------------------------------------- +subroutine compute_heat_propagation(hx, hy, dt, iter_max, save_flag, & + size_x, size_y, u_in, u_out) + + implicit none + double precision :: hx, hy, dt, error, prec + integer :: size_x, size_y, i, save_flag, iter_max + double precision, dimension(1:size_x, 1:size_y), intent(inout) :: u_in, u_out + + prec = 1e-4 + error = 1e10 + do i=0, iter_max + call heat(hx, hy, dt, size_x, size_y, u_in, u_out, error) + error = sqrt(error) + if (mod(i,10) == 0) then + print * , 'it =', i, 't = ', i*dt, 'err = ', error + endif + if (save_flag /= 0) then + write(name_save, '("sol_"i5.5)') i + call save_mat(name_save, size_x, size_y, u_out) + endif + u_in = u_out + if (error <= prec) exit + end do + +end subroutine compute_heat_propagation + +end program heat_seq diff --git a/TP/sources/session2/Fortran/lib/CMakeLists.txt b/TP/sources/session2/Fortran/lib/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..3cf42dfacc9cd971c42b062701a099ec1f81b713 --- /dev/null +++ b/TP/sources/session2/Fortran/lib/CMakeLists.txt @@ -0,0 +1 @@ +add_library(Heat SHARED heat.f90) diff --git a/TP/sources/session2/Fortran/lib/heat.f90 b/TP/sources/session2/Fortran/lib/heat.f90 new file mode 100644 index 0000000000000000000000000000000000000000..6754117f14c1efc5d4325e5b306079141d8b04dd --- /dev/null +++ b/TP/sources/session2/Fortran/lib/heat.f90 @@ -0,0 +1,51 @@ +!------------------------------------------------------------------------------- +!> @file heat.f90 +!> @author Inria SED Bordeaux +!> @brief Heat computation iteration code +!> +!> @details This file declares the heat() function. +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +!> Do a single iteration of the heat equation iterative computation on +!> a 2-D cartesian map. +!> +!> For each cell in @e u_in that is in <code>[1..size_x]x[1..size_y]</code>, +!> it put in @e u_out the result of one iteration computation using a +!> cross-stencil. +!> The complete equation can be found in the @e README file of the project. +!> +!> @param hx precision of the derivation over x +!> @param hy precision of the derivation over y +!> @param dt precision of the derivation over time +!> @param size_x the size of the cartesian map in x +!> @param size_y the size of the cartesion map in y +!> @param u_in the input map +!> @param u_out the output map +!> @return the square of the quadratic differences between @e u_in and @e u_out after +!> the execution of the heat function +!------------------------------------------------------------------------------- +subroutine heat(hx, hy, dt, size_x, size_y, u_in, u_out, error) + + implicit none + double precision hx, hy, dt + integer size_x, size_y, i, j + double precision, dimension(1:size_x, 1:size_y) :: u_in, u_out + double precision w_x, w_y, error, d + + w_x = dt / (hx * hx) + w_y = dt / (hy * hy) + d = 1.d0 - 2.d0 * w_x - 2.d0 * w_y + + error = 0.d0 + do i=2, size_x - 1 + do j = 2, size_y - 1 + u_out(i,j) = u_in(i,j) * d + & + (u_in(i - 1, j) + u_in(i + 1, j)) * w_x + & + (u_in(i, j - 1) + u_in(i, j + 1)) * w_y + error = error + (u_out(i,j) - u_in(i, j))**2 + end do + end do + +end subroutine heat + diff --git a/TP/sources/session2/Fortran/mat_utils.f90 b/TP/sources/session2/Fortran/mat_utils.f90 new file mode 100644 index 0000000000000000000000000000000000000000..83106a6cc85fca8035d98477339fc43da41d7cd8 --- /dev/null +++ b/TP/sources/session2/Fortran/mat_utils.f90 @@ -0,0 +1,49 @@ +subroutine print_mat(size_x, size_y, u) + implicit none + integer, intent(in) :: size_x, size_y + integer :: i,j + double precision, dimension(1:size_x, 1:size_y) :: u + write(6, *) '['; + do i=1, size_x + do j=1, size_y + write(6, '(e15.5)', advance='no') u(i,j) + end do + print * + end do + write(6, *) ']'; + print * + close(12) +end subroutine print_mat + +!------------------------------------------------------------------------------- +!> A matrix saving function +!> +!> It saves in a file the matrix given in the parameter @e u. +!> For instance if <code>u = [[1., 2.], [3., 4.]]</code>, then +!> save_mat("toto",2,2,u) will output in the file toto: +!> @code +!> 1.000000000000000e+00 2.000000000000000e+00 +!> 3.000000000000000e+00 4.000000000000000e+00 +!> @endcode +!> +!> @param filename the file to output the result +!> @param size_x the size in X of the matrix u +!> @param size_y the size in Y of the matrix u +!> @param u the matrix to print +!------------------------------------------------------------------------------- +subroutine save_mat(filename, size_x, size_y, u) + implicit none + character(len=*), intent(in) :: filename + integer, intent(in) :: size_x, size_y + integer :: i,j + double precision, dimension(1:size_x, 1:size_y) :: u + + open(unit=12,file=trim(filename),form='formatted') + do i=1, size_y + do j=1, size_y + write(12, '(e15.5)', advance='no') u(i,j) + end do + write(12, '(/)') + end do + close(12) +end subroutine save_mat diff --git a/TP/tp2.tex b/TP/tp2.tex new file mode 100644 index 0000000000000000000000000000000000000000..d1c92dca504a304e0330c3b94b9d5757109a1e41 --- /dev/null +++ b/TP/tp2.tex @@ -0,0 +1,204 @@ +% +% Session 2 of the Software Development Tools training session -- +% Hands-on session +% (10-11 Dec. 2012) +% +% Compilation: rubber -I "$INRIA_ROOT" --pdf tp1 +% +\documentclass{beamer} + +\usepackage{helvet} + +\usepackage[utf8]{inputenc} +\usepackage{hyperref,xspace,multicol} +\usepackage[absolute,overlay]{textpos} +\usepackage{tikz} +\usepackage{listings} +\usepackage{ulem} +\usepackage{relsize} +\usetikzlibrary{arrows,shapes,shadows} + +% Remember the position of every picture. +\tikzstyle{every picture}+=[remember picture] + +% \inriaswitchcolors COLOR +% +% Where COLOR is one of red, blue, orange, darkblue, violet, +% pastelgreen, grey, or green. +\newcommand{\inriaswitchcolors}[1]{% +\pgfaliasimage{figfootline}{figfootline-#1}% !!! +\pgfaliasimage{figbackground}{figbackground-#1}% !!! +\pgfaliasimage{figbackground}{figbackground-#1}% !!! +} + +\AtBeginSection[]{ + \begin{frame}[plain] + \partpage + \end{frame} +} +\newcommand{\sect}[2]{% +\inriaswitchcolors{#2}% +\section{#1}% +} +\newcommand{\subsect}[1]{% +\begin{frame} + \center\huge\alert{\textbf{#1}} +\end{frame} +} + +\title{Hands on CMake} +\author{Marc Fuentes} +\institute{SED Bordeaux Sud-Ouest} +\date{26 Novembre 2018} +\definecolor{dkgreen}{rgb}{0,0.6,0} +\definecolor{gray}{rgb}{0.5,0.5,0.5} +\definecolor{mauve}{rgb}{0.58,0,0.82} +\include{cmake_def} +\lstset{ % + language=mycmake, + framerule=0pt, + basicstyle=\relsize{-2}\ttfamily, % the size of the fonts that are used for the code + showspaces=false, % show spaces adding particular underscores + showstringspaces=false, % underline spaces within strings + showtabs=false, % show tabs within strings adding particular underscores + %frame=single, % adds a frame around the code + rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. commens (green here)) + breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace + keywordstyle=\color{blue}, % keyword style + commentstyle=\color{dkgreen}, % comment style + stringstyle=\color{mauve} +} +\begin{document} + +\begin{frame}[plain] + \titlepage +\end{frame} + +\inriaswitchcolors{grey} + +\begin{frame} + \frametitle{Outline of the session} + \begin{itemize}[<+->] + \item basic \texttt{CMakeLists.txt} + \item add a personal library + \item system inspection + \item add some tests + \item parametrize installation + \item make a source package + \end{itemize} + \begin{center} + \end{center} +\end{frame} + +\begin{frame}[fragile] + \frametitle{Basic {\tt CMakeLists.txt}} + \begin{itemize} + \item open \texttt{CMakeLists.txt} + \lstinputlisting{sources/session2/Fortran/CMakeLists.txt} + \item \texttt{mkdir build ; cd build ; cmake .. } (try with \texttt{ccmake ..} or \texttt{cmake-gui ..}) + \item try to change the build type (\texttt{Release}, \texttt{RelWithDebInfo}) \\ + \texttt{cmake -DCMAKE\_BUILD\_TYPE=Debug ..} + \item run {\tt Make}: \texttt{make} \\ + \texttt{make VERBOSE=1} if you need some info + \end{itemize} +\end{frame} + +\begin{frame}[fragile] + \frametitle{Add a personal library} + \begin{itemize}[<+->] + \item we want to put \texttt{heat.f90} (\texttt{heat.c}) in a library: + \item move \texttt{heat.f90} (\texttt{heat.c}) to a subdirectory \texttt{lib} and remove it from \lstinline[basicstyle=\relsize{0}]!add_executable! + \item use \lstinline[basicstyle=\relsize{0}]!add_subdirectory(lib)! + \item for C: move \texttt{heat.h} in a subdirectory \texttt{include} and add with\lstinline[basicstyle=\relsize{0}]!include_directories(include)! + \item create in \texttt{lib} subdirectory a \texttt{CMakeLists.txt} containing + \begin{itemize} + \item \lstinline[basicstyle=\relsize{0}]!add_library(Heat SHARED heat.f90)! for Fortran + \item \lstinline[basicstyle=\relsize{0}]!add_library(Heat SHARED heat.c)! for C + \end{itemize} + \item in the main \texttt{CMakeLists.txt}: \\ \lstinline[basicstyle=\relsize{0}]!target_link_libraries(patcHeatSeq Heat)! + \item \textcolor{gray}{optional: try to switch between \texttt{SHARED} and \texttt{STATIC} in} \lstinline[basicstyle=\relsize{0}]!add_library! + \end{itemize} +\end{frame} + +\begin{frame}[fragile] + \frametitle{System inspection} + \begin{itemize}[<+->] + \item we want to use \texttt{MPI} for the parallel version + \item look at \texttt{/usr/share/cmake/Modules/FindMPI.cmake} + \item use \lstinline[basicstyle=\relsize{0}]!find_package(MPI)! + \item add a conditional block + \begin{lstlisting} +if(MPI_FOUND) + include_directories(${MPI_INCLUDE_PATH}) + set(HEAT_PAR_SOURCES heat_par.c mat_utils.c) + add_executable(patcHeatPar ${HEAT_PAR_SOURCES} mat_utils.h) + target_link_libraries(patcHeatPar ${MPI_LIBRARIES} m Heat) +else (MPI_FOUND) + message(WARNING "MPI not found; only the sequential version will be built") +endif(MPI_FOUND) +\end{lstlisting} + \end{itemize} +\end{frame} + +\begin{frame}[fragile] + \frametitle{Add some tests} + \begin{itemize}[<+->] + \item first you need \lstinline[basicstyle=\relsize{0}]!enable_testing()! + \item you could add tests with \lstinline[basicstyle=\relsize{0}]!add_test(name_test command arg1 arg2 ...)! + \item you can filter the output with \\ + \lstinline[basicstyle=\relsize{-2}]!set_tests_properties(name_test PROPERTIES PASS_REGULAR_EXPRESSION "regexp")! + \item try the following tests: + \begin{itemize} + \item \lstinline!add_test(patcHeatSeqUsage ./patcHeatSeq )! + \\ \lstinline!set_tests_properties(patcHeatSeqUsage PROPERTIES PASS_REGULAR_EXPRESSION "Usage*")! + \item \lstinline! add_test(patcHeatSeq10Err ./patcHeatSeq 10 10 1 0)! + \\ \lstinline!set_tests_properties(patcHeatSeq10Err PROPERTIES PASS_REGULAR_EXPRESSION "1.732*")! + \item \lstinline!add_test(patcHeatPar4 ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4! + \\ \lstinline! ./patcHeatPar 10 10 200 2 2 0)! + \end{itemize} + \item run \texttt{make test} or try \texttt{ctest -VV } + \end{itemize} +\end{frame} + +\begin{frame}[fragile] + \frametitle{Parametrize installation} + \begin{itemize}[<+->] + \item edit the rpath config + \begin{lstlisting} +# set the rpath config +set(CMAKE_SKIP_BUILD_RPATH FALSE) +set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE) +set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") + \end{lstlisting} + \item install binaries in their directories: + \begin{itemize} + \item \texttt{patcHeatSeq} $\rightarrow$ \texttt{bin} \\ + \lstinline[basicstyle=\relsize{0}]!install(TARGETS patcHeatSeq DESTINATION bin)! + \item \texttt{patcHeatPar} $\rightarrow$ \texttt{bin} \\ + \lstinline[basicstyle=\relsize{0}]!install(TARGETS patcHeatPar DESTINATION bin)! + \item \texttt{libHeat} $\rightarrow$ \texttt{lib} \\ + \lstinline[basicstyle=\relsize{0}]!install(TARGETS Heat DESTINATION lib)! + \end{itemize} + \item edit \texttt{CMAKE\_INSTALL\_PREFIX}: \texttt{make edit\_cache} + \item do the install: \texttt {make install} + \end{itemize} +\end{frame} + +\begin{frame}[fragile] + \frametitle{Make a source package} + \begin{itemize}[<+->] + \item avoid inclusion of build directory in package: \\ + \lstinline[basicstyle=\relsize{0}]!set(CPACK_SOURCE_IGNORE_FILES "/build/" "/.*\\\\.swp$")! + \item set some \texttt{CPACK} variables: + \begin{itemize} + \item \lstinline!set(CPACK_SOURCE_GENERATOR "TGZ")! + \item \lstinline!set(CPACK_PACKAGE_NAME "mon_super_package_qui_dechire")! + \end{itemize} + \item add \lstinline[basicstyle=\relsize{0}]!include(CPack)! + \item run \texttt{make package\_source} + \item try to install the tarball in another directory + \item \alert{caveat}: binaries are often non portable, esp. on free OSes + \end{itemize} +\end{frame} + +\end{document}