first draft

.gitignore

+*.txt
+/main

cart.c

+#include <mpi.h> 
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "project.h"
+
+void Optimize_Cart_Cluster(int dim0_size, int dim1_size, MPI_Comm comm, int rank, int *pro_per_dim, int *cell_per_pro) {
+	int num_p;
+	float float_pro_per_dim0, float_pro_per_dim1;
+
+	if(MPI_Comm_size(comm, &num_p)) {
+		fprintf(stderr, "Cannot fetch size of cluster\n");		
+		exit(1);
+	}
+
+	if(rank == 0) {
+		printf("number of processes: %d\n", num_p);
+	}
+
+	float_pro_per_dim1 = sqrt((float)(num_p*dim0_size)/dim1_size);
+	float_pro_per_dim0 = num_p/float_pro_per_dim1;
+	float_pro_per_dim1 = floor(float_pro_per_dim1);
+	float_pro_per_dim0 = floor(float_pro_per_dim0);		// try ceil
+	pro_per_dim[0] = (int)float_pro_per_dim0;
+	pro_per_dim[1] = (int)float_pro_per_dim1;
+
+	cell_per_pro[0] = ceil(dim0_size/(float)pro_per_dim[1]);
+	cell_per_pro[1] = ceil(dim1_size/(float)pro_per_dim[0]);
+
+	if(rank==0) {
+		printf("dim0: %d dim1: %d\n", pro_per_dim[0], pro_per_dim[1]);
+		printf("size per pro: %dx%d\n", cell_per_pro[0], cell_per_pro[1]);	
+	}
+}
+
+MPI_Comm Create_MPI_Cart_Cluster(MPI_Comm comm, int rank, int *pro_per_dim) {
+	MPI_Comm cart_comm;
+	int periods[] = {0,0};	// edges are not connected
+	
+	if(MPI_Cart_create(comm, 2, pro_per_dim, periods, 0, &cart_comm)) {
+		fprintf(stderr, "Cannot create topology\n");		
+		exit(1);
+	}
+	if(cart_comm == MPI_COMM_NULL) {
+		printf("process %d not in use. exiting...\n", rank);
+		MPI_Finalize();
+		exit(0);
+	}
+
+	return cart_comm;
+}
\ No newline at end of file

jacobi.c

+#define min(X,Y) (((X) < (Y)) ? (X) : (Y))
+#define max(X,Y) (((X) > (Y)) ? (X) : (Y))
+
+#include <omp.h>
+#include <mpi.h> 
+#include <stdio.h>
+#include <stdlib.h>
+#include "project.h"
+
+void Alloc_Partial_Field(int *matrix_size, double ***partial_field, double ***partial_field_clipboard) {
+	*partial_field = New_Matrix(matrix_size[0], matrix_size[1]);
+	if (partial_field == NULL) {
+		fprintf(stderr, "Can't allocate partial_field\n");
+		exit(1);
+	}
+	*partial_field_clipboard = New_Matrix(matrix_size[0], matrix_size[1]);
+	if (partial_field_clipboard == NULL) {
+		fprintf(stderr, "Can't allocate partial_field_clipboard\n");
+		exit(1);
+	}
+
+	Init_Matrix(*partial_field, matrix_size[0], matrix_size[1], 0);
+	Init_Matrix(*partial_field_clipboard, matrix_size[0], matrix_size[1], 0);
+}
+
+void Init_Jacobi(int dim0_size, int dim1_size, int alpha, double *delta_t, double *hx, double *hy, double *hx_square, double *hy_square) {
+	*hx = 1.0/(double)dim0_size;
+	*hy = 1.0/(double)dim1_size;
+	*hx_square = *hx * *hx;
+	*hy_square = *hy * *hy;
+
+	double max_delta_t = 0.25*((min(*hx,*hy))*(min(*hx,*hy)))/alpha;  /* minimaler Wert für Konvergenz */
+	if (*delta_t > max_delta_t) { 
+		*delta_t = max_delta_t;
+		printf ("Info: delta_t set to %.10lf.\n", *delta_t);
+	}
+}
+
+void Init_Edges(int dim0_size, int dim1_size, int *matrix_size, int neighbor_dim0_left, int neighbor_dim0_right, int neighbor_dim1_left, int neighbor_dim1_right, double **partial_field, double **partial_field_clipboard, t_process_info pi) {
+	if(neighbor_dim1_left == MPI_PROC_NULL) {
+		for(int i = pi.start_m; i <= pi.end_m; i++) {
+			partial_field[i - pi.start_m + 1][1] = (double)i / (dim0_size-1);
+			partial_field_clipboard[i - pi.start_m + 1][1] = (double)i / (dim0_size-1);
+		}
+	}
+	if(neighbor_dim1_right == MPI_PROC_NULL) {
+		for(int i = pi.start_m; i <= pi.end_m; i++) {
+			partial_field[i - pi.start_m + 1][matrix_size[1]-2] = 1 - (double)i / (dim0_size-1);
+			partial_field_clipboard[i - pi.start_m + 1][matrix_size[1]-2] = 1 - (double)i / (dim0_size-1);
+		}
+	}
+	if(neighbor_dim0_left == MPI_PROC_NULL) {
+		for(int i = pi.start_n; i <= pi.end_n; i++) {
+			partial_field[1][i - pi.start_n + 1] = (double)i / (dim1_size-1);
+			partial_field_clipboard[1][i - pi.start_n + 1] = (double)i / (dim1_size-1);
+		}
+	}
+	if(neighbor_dim0_right == MPI_PROC_NULL) {
+		for(int i = pi.start_n; i <= pi.end_n; i++) {
+			partial_field[matrix_size[0] - 2][i - pi.start_n + 1] = 1 - (double)i / (dim1_size-1);
+			partial_field_clipboard[matrix_size[0] - 2][i - pi.start_n + 1] = 1 - (double)i / (dim1_size-1);
+		}
+	}
+}
+
+int Jacobi_Iterate(int neighbor_dim0_left, int neighbor_dim0_right, int neighbor_dim1_left, int neighbor_dim1_right, double alpha, double delta_t, double eps, double hx_square, double hy_square, t_process_info pi, double ***partial_field, double ***partial_field_clipboard) {
+	double delta_a;
+	double **swap;
+	double maxdiff = 0;
+	double inner_maxdiff;
+	for(
+		int i = (neighbor_dim0_left == MPI_PROC_NULL) ? 2 : 1; // catch edges
+		i < pi.end_m - pi.start_m + ((neighbor_dim0_right == MPI_PROC_NULL) ? 1 : 2); 
+		i++
+	) {
+		inner_maxdiff = 0;
+		#pragma omp parallel for private(delta_a) reduction(max: inner_maxdiff) schedule(dynamic, 512) default(shared) num_threads(48)
+		for(
+			int j = (neighbor_dim1_left == MPI_PROC_NULL) ? 2 : 1; // catch edges
+			j < pi.end_n - pi.start_n + ((neighbor_dim1_right == MPI_PROC_NULL) ? 1 : 2); 
+			j++
+		) {
+			//printf("%d working on %d, %d\n", omp_get_thread_num(), i, j);
+			delta_a = alpha * 
+				    ( ((*partial_field)[i+1][j] + (*partial_field)[i-1][j] - 2.0 * (*partial_field)[i][j]) / (hy_square)
+					 +((*partial_field)[i][j-1] + (*partial_field)[i][j+1] - 2.0 * (*partial_field)[i][j]) / (hx_square) );
+			delta_a = delta_a * delta_t;
+			(*partial_field_clipboard)[i][j] = (*partial_field)[i][j] + delta_a;
+			
+			if(delta_a > inner_maxdiff)
+				inner_maxdiff = delta_a;
+		}
+		if(inner_maxdiff > maxdiff)
+			maxdiff = inner_maxdiff;
+	}
+	swap = *partial_field_clipboard;
+	*partial_field_clipboard = *partial_field;
+	*partial_field = swap;
+
+	return maxdiff <= eps;
+}
\ No newline at end of file

main.c

+
+#define min(X,Y) (((X) < (Y)) ? (X) : (Y))
+#define max(X,Y) (((X) > (Y)) ? (X) : (Y))
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <mpi.h> 
+#include <sys/time.h>
+#include "project.h"
+
+int main(int argc, char **argv)
+{
+	double eps = 0.001;
+	double delta_t = 0.000001;
+	double alpha = 1;
+	int dim0_size, dim1_size;
+	double **root_field;				// complete field owned by root
+	double **partial_field;				// partial field where process works on
+	double **partial_field_clipboard;	// copy of partial field
+	t_process_info pi;
+	t_process_info *infos;
+	int pro_per_dim[2];
+	int cell_per_pro[2];
+	MPI_Comm cart_comm;
+	int matrix_size[2];
+	int neighbor_dim0_left, neighbor_dim0_right, neighbor_dim1_left, neighbor_dim1_right;
+	double hx, hy, hx_square, hy_square;
+	double *dim1_own_edge_values;
+	double *dim1_neighbor_egde_values;
+	MPI_Request sync_requests[9];	// 2 for each edge + 1 completion
+	int iterations = 0;
+	double time_start, time_iterate, time_end;
+
+	MPI_Init(&argc, &argv);
+
+	Process_Args(argc, argv, &dim0_size, &dim1_size, &eps, &delta_t);
+
+	int rank, cart_cluster_size;
+	if(MPI_Comm_rank(MPI_COMM_WORLD, &rank)) {
+		fprintf(stderr, "Cannot fetch rank\n");
+		exit(1);
+	}
+
+	if(rank == 0) {
+		time_start = MPI_Wtime();
+	}
+
+	if(rank == 0) {
+		root_field = New_Matrix(dim0_size, dim1_size);
+		if (root_field == NULL) {
+			fprintf(stderr, "rank %s: Can't allocate root_field !\n", rank);
+			exit(1);
+		}
+	}
+
+	// optimize cart cluster
+	Optimize_Cart_Cluster(dim0_size, dim1_size, MPI_COMM_WORLD, rank, pro_per_dim, cell_per_pro);
+
+	cart_comm = Create_MPI_Cart_Cluster(MPI_COMM_WORLD, rank, pro_per_dim);
+
+	pi = Calculate_Process_Info(cart_comm, rank, dim0_size, dim1_size, cell_per_pro);
+
+	matrix_size[0] = pi.end_m - pi.start_m + 3;
+	matrix_size[1] = pi.end_n - pi.start_n + 3;
+
+	if(MPI_Comm_size(cart_comm, &cart_cluster_size)) {
+		fprintf(stderr, "Cannot fetch size of cart\n");		
+		exit(1);
+	}
+
+	infos = Gather_Process_Info(&pi, rank, cart_cluster_size, cart_comm);
+
+	if(rank == 0) {
+		for(int i = 0; i < cart_cluster_size; i++) {
+			Print_Process_Info(infos[i]);
+		}
+	}
+
+	Alloc_Partial_Field(matrix_size, &partial_field, &partial_field_clipboard);
+
+	Init_Neighbor_Comm(cart_comm, sync_requests, matrix_size, &neighbor_dim0_left, &neighbor_dim0_right, &neighbor_dim1_left, &neighbor_dim1_right, &dim1_own_edge_values, &dim1_neighbor_egde_values);
+
+	Init_Jacobi(dim0_size, dim1_size, alpha, &delta_t, &hx, &hy, &hx_square, &hy_square);
+
+	Init_Edges(dim0_size, dim1_size, matrix_size, neighbor_dim0_left, neighbor_dim0_right, neighbor_dim1_left, neighbor_dim1_right, partial_field, partial_field_clipboard, pi);
+
+	int *completions = malloc(sizeof(int) * cart_cluster_size);
+
+	/*
+	*
+	* START ITERATION
+	*
+	*/
+	if(rank == 0) {
+		time_iterate = MPI_Wtime();
+	}
+
+	while (1) {			// iterate until break;
+		iterations++;
+
+		int completion = Jacobi_Iterate(neighbor_dim0_left, neighbor_dim0_right, neighbor_dim1_left, neighbor_dim1_right, alpha, delta_t, eps, hx_square, hy_square, pi, &partial_field, &partial_field_clipboard);
+
+		/*
+		if(iterations == 2) {
+			MPI_Finalize();
+			return 0;
+		}
+		*/
+
+		int all_completed = Sync(cart_comm, completion, completions, cart_cluster_size, matrix_size, sync_requests, neighbor_dim0_left, neighbor_dim0_right, neighbor_dim1_left, neighbor_dim1_right, partial_field, dim1_own_edge_values, dim1_neighbor_egde_values);
+
+		if(all_completed) {
+			printf("rank: %d: break after %d iterations\n", rank, iterations);
+			break;
+		}
+	}
+	if(rank == 0) {
+		time_end = MPI_Wtime();
+	}
+
+	/*
+	*
+	* END ITERATION
+	*
+	*/
+
+	if(rank == 0) {
+		printf("init: %.10lf\n", time_iterate- time_start);
+		printf("iterate %.10lf\n", time_end - time_iterate);
+		printf("total: %.10lf\n", time_end - time_start);
+	}
+
+	Send_To_Root(cart_comm, rank, dim0_size, dim1_size, cart_cluster_size, matrix_size, infos, partial_field, root_field);
+
+	MPI_Finalize();
+
+	return 0;
+}
\ No newline at end of file

makefile

+all: main
+
+main: matrix.c proj.c args.c pid0.c pi.c cart.c mpi_util.c
+	mpicc matrix.c main.c args.c pid0.c pi.c cart.c mpi_util.c jacobi.c -std=c99 -lm -o main
+
+ViewMatrix.class: ViewMatrix.java
+	javac ViewMatrix.java
+
+clean:
+	rm -f *.o *.c~ heat Matrix.txt
+	rm -f ViewMatrix.class
+ 
\ No newline at end of file

pi.c

+#include <mpi.h> 
+#include <stdio.h>
+#include <stdlib.h>
+#include "project.h"
+
+t_process_info Calculate_Process_Info(MPI_Comm cart_comm, int rank, int dim0_size, int dim1_size, int *cell_per_pro) {
+	t_process_info pi;
+	int coord[2];
+	if(MPI_Cart_coords(cart_comm, rank, 2, coord)) {
+		fprintf(stderr, "Cannot get coordinates\n");		
+		exit(1);
+	}
+
+	// calculate own field using coord
+	pi.start_m = coord[0] * cell_per_pro[0];
+	pi.end_m = (coord[0]+1)*cell_per_pro[0] -1;
+	pi.start_n = coord[1] * cell_per_pro[1];
+	pi.end_n = (coord[1]+1) * cell_per_pro[1] -1;
+	if(pi.end_m > dim0_size - 1) {
+		pi.end_m = dim0_size - 1;
+	}
+	if(pi.end_n > dim1_size - 1) {
+		pi.end_n = dim1_size - 1;
+	}
+	pi.coord0 = coord[0];
+	pi.coord1 = coord[1];
+	pi.rank = rank;
+
+	return pi;
+}
+
+void Print_Process_Info(t_process_info pi) {
+	printf("rank: %d->(%d,%d) from (%d, %d) to (%d,%d)\n", 
+		pi.rank,
+		pi.coord0,
+		pi.coord1,
+		pi.start_m,
+		pi.start_n,
+		pi.end_m,
+		pi.end_n
+	);
+}
\ No newline at end of file

pid0.c

+#include <mpi.h> 
+#include <stdio.h>
+#include <stdlib.h>
+#include "project.h"
+
+t_process_info* Gather_Process_Info(t_process_info *pi, int rank, int cluster_size, MPI_Comm cart_comm) {
+	t_process_info *infos;
+	MPI_Datatype MPI_process_info;
+
+	Create_MPI_Type_t_process_info(&MPI_process_info);
+
+	if(rank == 0) {
+		infos = malloc(sizeof(t_process_info) * cluster_size);
+	}
+
+	if(MPI_Gather(pi, 1, MPI_process_info, infos, 1, MPI_process_info, 0, cart_comm)) {
+		fprintf(stderr, "Gather failed\n");
+		exit(1);
+	}
+
+	return infos;
+}
+
+void Send_To_Root(MPI_Comm cart_comm, int rank, int dim0_size, int dim1_size, int cart_cluster_size, int *matrix_size, t_process_info *infos, double **partial_field, double **root_field) {
+	MPI_Request send_request;
+	MPI_Isend(partial_field[0], matrix_size[0]*matrix_size[1], MPI_DOUBLE, 0, 0, cart_comm, &send_request);
+	if(rank == 0) {
+		MPI_Request *requests = malloc(sizeof(MPI_Request) * cart_cluster_size);
+		double **allocation = malloc(sizeof(double*) * cart_cluster_size);
+		for(int i = 0; i < cart_cluster_size; i++) {
+			allocation[i] = malloc(sizeof(double) * (infos[i].end_m - infos[i].start_m + 3) * (infos[i].end_n - infos[i].start_n + 3));
+			MPI_Irecv(allocation[i], 
+				(infos[i].end_m - infos[i].start_m + 3) * (infos[i].end_n - infos[i].start_n + 3), 
+				MPI_DOUBLE, 
+				infos[i].rank, 
+				MPI_ANY_TAG, 
+				cart_comm,
+				&requests[i]
+			);
+		}
+		for(int i = 0; i < cart_cluster_size; i++) {
+			int current;
+			MPI_Waitany(cart_cluster_size, requests, &current, MPI_STATUS_IGNORE);
+			Insert_Array_In_Matrix(
+				root_field, 
+				dim0_size, 
+				dim1_size, 
+				infos[current].start_m, 
+				infos[current].start_n, 
+				allocation[current], 
+				infos[current].end_m - infos[current].start_m + 3, 
+				infos[current].end_n - infos[current].start_n + 3, 
+				1, 1, 1, 1);
+			free(allocation[current]);
+		}
+		free(requests);
+		free(allocation);
+		Write_Matrix(root_field, dim0_size, dim1_size);
+	}
+	MPI_Wait(&send_request, MPI_STATUS_IGNORE);
+}
\ No newline at end of file