mpi first draft

Matrix.txt

-4
-4
-
-8.0000000000
-8.0000000000
-8.0000000000
-8.0000000000
-
-8.0000000000
-8.0000000000
-8.0000000000
-8.0000000000
-
-8.0000000000
-8.0000000000
-8.0000000000
-8.0000000000
-
-8.0000000000
-8.0000000000
-8.0000000000
-8.0000000000
+8
+8
+
+0.0000000000
+0.1428571429
+0.2857142857
+0.4285714286
+0.5714285714
+0.7142857143
+0.8571428571
+1.0000000000
+
+0.1428571429
+0.1571330166
+0.1926800492
+0.2544795945
+0.3424186881
+0.4656438772
+0.6390629773
+0.8571428571
+
+0.2857142857
+0.1926800492
+0.1411981312
+0.1397585075
+0.1844451573
+0.2859846957
+0.4656438772
+0.7142857143
+
+0.4285714286
+0.2544795945
+0.1397585075
+0.0946683354
+0.1080432823
+0.1844451573
+0.3424186881
+0.5714285714
+
+0.5714285714
+0.3424186881
+0.1844451573
+0.1080432823
+0.0946683354
+0.1397585075
+0.2544795945
+0.4285714286
+
+0.7142857143
+0.4656438772
+0.2859846957
+0.1844451573
+0.1397585075
+0.1411981312
+0.1926800492
+0.2857142857
+
+0.8571428571
+0.6390629773
+0.4656438772
+0.3424186881
+0.2544795945
+0.1926800492
+0.1571330166
+0.1428571429
+
+1.0000000000
+0.8571428571
+0.7142857143
+0.5714285714
+0.4285714286
+0.2857142857
+0.1428571429
+0.0000000000
 

proj.c

@@ -60,8 +60,6 @@ int main(int argc, char **argv)
 			fprintf(stderr, "PID %s: Can't allocate root_field !\n", pid);
 			exit(1);
 		}
-
-		Init_Matrix(root_field, m, n, 0);
 	}
 
 	if(pid == 0) {
@@ -167,12 +165,25 @@ int main(int argc, char **argv)
 		exit(1);
 	}
 
+	Init_Matrix(partial_field, matrix_size[0], matrix_size[1], 0);
+
 	double **partial_field_tmp = New_Matrix(matrix_size[0], matrix_size[1]);
+	Init_Matrix(partial_field_tmp, matrix_size[0], matrix_size[1], 0);
 	double **swap;
 	double hx = 1.0/(double)m;
 	double hy = 1.0/(double)n;
 	double hx_square = hx * hx;
 	double hy_square = hy * hy;
+	double maxdiff;
+
+	double *dim1_own_edge_values = malloc(sizeof(double) * (2*(matrix_size[0]-2)));
+	double *dim1_neighbor_egde_values = malloc(sizeof(double) * (2*(matrix_size[0]-2)));
+	MPI_Request sync_requests[9];	// 2 for each edge + 1 completion
+
+	// set to MPI null for waitany -- only needed once as others are overwritten with every iteration
+	for(i = 0; i < 9; i++) {
+		sync_requests[i] = MPI_REQUEST_NULL;
+	}
 
 	double max_delta_t = 0.25*((min(hx,hy))*(min(hx,hy)))/alpha;  /* minimaler Wert für Konvergenz */
 	if (delta_t > max_delta_t) { 
@@ -181,25 +192,170 @@ int main(int argc, char **argv)
 			printf ("Info: delta_t set to %.10lf.\n", delta_t);
 	}
 
-	for(i = 1; i < pi.end_m - pi.start_m + 2; i++) {		// catch edges
-		for(j = 1; j < pi.end_n - pi.start_n + 2; j++) {	// catch edges
-			/*
-			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_tmp[i][j] = partial_field[i][j] + delta_a;
-			
-			if(delta_a > maxdiff)
-				maxdiff = delta_a;
-			*/
-			partial_field_tmp[i][j] = 8;
-		}
-	}
-	swap = partial_field_tmp;
-	partial_field_tmp = partial_field;
-	partial_field = swap;
-	
+	int neighbor_dim0_left, neighbor_dim0_right, neighbor_dim1_left, neighbor_dim1_right;
+	if(MPI_Cart_shift(cart_comm, 0, 1, &neighbor_dim0_left, &neighbor_dim0_right)) {
+		fprintf(stderr, "Shift failed\n");		
+		exit(1);
+	}
+	if(MPI_Cart_shift(cart_comm, 1, 1, &neighbor_dim1_left, &neighbor_dim1_right)) {
+		fprintf(stderr, "Shift failed\n");		
+		exit(1);
+	}
+
+	//init edges
+	if(neighbor_dim1_left == MPI_PROC_NULL) {
+		for(i = pi.start_m; i <= pi.end_m; i++) {
+			partial_field[i - pi.start_m + 1][1] = (double)i / (m-1);
+			partial_field_tmp[i - pi.start_m + 1][1] = (double)i / (m-1);
+		}
+	}
+	if(neighbor_dim1_right == MPI_PROC_NULL) {
+		for(i = pi.start_m; i <= pi.end_m; i++) {
+			partial_field[i - pi.start_m + 1][matrix_size[1]-2] = 1 - (double)i / (m-1);
+			partial_field_tmp[i - pi.start_m + 1][matrix_size[1]-2] = 1 - (double)i / (m-1);
+		}
+	}
+	if(neighbor_dim0_left == MPI_PROC_NULL) {
+		for(i = pi.start_n; i <= pi.end_n; i++) {
+			partial_field[1][i - pi.start_n + 1] = (double)i / (n-1);
+			partial_field_tmp[1][i - pi.start_n + 1] = (double)i / (n-1);
+		}
+	}
+	if(neighbor_dim0_right == MPI_PROC_NULL) {
+		for(i = pi.start_n; i <= pi.end_n; i++) {
+			partial_field[matrix_size[0] - 2][i - pi.start_n + 1] = 1 - (double)i / (n-1);
+			partial_field_tmp[matrix_size[0] - 2][i - pi.start_n + 1] = 1 - (double)i / (n-1);
+		}
+	}
+	int *completions = malloc(sizeof(int) * num_p);
+	int k = 0;
+
+	/*
+	*
+	* START ITERATION
+	*
+	*/
+
+	while (1) {			// iterate until break;
+		k++;
+		maxdiff = 0;
+		for(
+			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++
+		) {	
+			for(
+				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++
+			) {
+				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_tmp[i][j] = partial_field[i][j] + delta_a;
+				
+				if(delta_a > maxdiff)
+					maxdiff = delta_a;
+			}
+		}
+		swap = partial_field_tmp;
+		partial_field_tmp = partial_field;
+		partial_field = swap;
+
+		printf("comp pid %d: max: %.10lf eps: %.10lf\n", pid, maxdiff, eps);
+		int completion = maxdiff <= eps;
+
+		printf("pid %d completion %d\n", pid, completion);
+		if(MPI_Allgather(&completion, 1, MPI_INT, completions, 1, MPI_INT, cart_comm)) {
+			fprintf(stderr, "Alltoall failed\n");		
+			exit(1);
+		}
+		for(i = 0; i < num_p; i++) {
+			printf("pid %d: %d -> %d \n", pid, i, completions[i]);
+		}
+		int all_completed = 1;
+		for(i = 0; i < num_p; i++) {
+			if(!completions[i]) {
+				all_completed = 0;
+				break;
+			}
+		}
+		if(all_completed) {
+			printf("pid: %d: break after %d iterations\n", pid, k);
+			break;
+		}
+
+		// Sync edges
+
+		// copy own edges in send buffer
+		if(neighbor_dim0_left != MPI_PROC_NULL) {
+			MPI_Isend(&(partial_field[1][1]), (matrix_size[1] - 2), MPI_DOUBLE, neighbor_dim0_left, 0, cart_comm, &(sync_requests[0]));
+			MPI_Irecv(&(partial_field[0][1]), (matrix_size[1] - 2), MPI_DOUBLE, neighbor_dim0_left, 0, cart_comm, &(sync_requests[4]));
+			//memcpy(&(dim1_own_edge_values[2*(matrix_size[0]-2)]), &(partial_field[1][1]), sizeof(double) * (matrix_size[1] - 2));
+		}
+		if(neighbor_dim0_right != MPI_PROC_NULL) {
+			MPI_Isend(&(partial_field[matrix_size[0]-2][1]), (matrix_size[1] - 2), MPI_DOUBLE, neighbor_dim0_right, 0, cart_comm, &(sync_requests[1]));
+			MPI_Irecv(&(partial_field[matrix_size[0]-1][1]), (matrix_size[1] - 2), MPI_DOUBLE, neighbor_dim0_right, 0, cart_comm, &(sync_requests[5]));
+			//memcpy(&(dim1_own_edge_values[2*(matrix_size[0]-2) + matrix_size[1] - 2]), &(partial_field[matrix_size[0]-2][1]), sizeof(double) * (matrix_size[1] - 2));
+		}
+		if(neighbor_dim1_left != MPI_PROC_NULL) {
+			for(i = 0; i < matrix_size[0] - 2; i++) {
+				dim1_own_edge_values[i] = partial_field[i+1][1];
+			}
+			MPI_Isend(&(dim1_own_edge_values[0]), matrix_size[0] - 2, MPI_DOUBLE, neighbor_dim1_left, 0, cart_comm, &(sync_requests[2]));
+			MPI_Irecv(&(dim1_neighbor_egde_values[0]), matrix_size[0] -2, MPI_DOUBLE, neighbor_dim1_left, 0, cart_comm, &(sync_requests[6]));
+		}
+		if(neighbor_dim1_right != MPI_PROC_NULL) {
+			int right_edge_index = matrix_size[1]-2;
+			for(i = 0; i < matrix_size[0] - 2; i++) {
+				dim1_own_edge_values[matrix_size[0]-2+i] = partial_field[i+1][right_edge_index];
+			}
+			MPI_Isend(&(dim1_own_edge_values[matrix_size[0]-2]), matrix_size[0] - 2, MPI_DOUBLE, neighbor_dim1_right, 0, cart_comm, &(sync_requests[3]));
+			MPI_Irecv(&(dim1_neighbor_egde_values[matrix_size[0]-2]), matrix_size[0] -2, MPI_DOUBLE, neighbor_dim1_right, 0, cart_comm, &(sync_requests[7]));
+		}
+
+		while(1) {
+			int current;
+			MPI_Waitany(9, sync_requests, &current, MPI_STATUS_IGNORE);
+			if(current == MPI_UNDEFINED) {
+				break;
+			}
+			if(current == 6) {
+				for(i = 0; i < matrix_size[0] - 2; i++) {
+					partial_field[i+1][0] = dim1_neighbor_egde_values[i];
+				}
+			}
+			if(current == 7) {
+				int right_edge_index = matrix_size[1] - 1;
+				for(i = 0; i < matrix_size[0] - 2; i++) {
+					partial_field[i+1][right_edge_index] = dim1_neighbor_egde_values[matrix_size[0]-2+i];
+				}
+			}
+			if(current == 8) {
+				int all_completed = 1;
+				for(i = 0; i < num_p; i++) {
+					if(!completions[i]) {
+						all_completed = 0;
+						break;
+					}
+				}
+				if(all_completed) {
+					printf("pid: %d: break after %d iterations\n", pid, k);
+					break;
+				}
+			}
+			printf("waitany...\n");
+		}
+		printf("waitany completed\n");
+	}
+
+	/*
+	*
+	* END ITERATION
+	*
+	*/
+
 	//Send_To_Root(partial_field, pi.end_m - pi.start_m + 2, pi.end_n - pi.start_n + 2);
 	MPI_Send(partial_field[0], matrix_size[0]*matrix_size[1], MPI_DOUBLE, 0, 0 ,cart_comm);
 	if(pid == 0) {
@@ -219,6 +375,11 @@ int main(int argc, char **argv)
 		for(i = 0; i < num_p; i++) {
 			int current;
 			MPI_Waitany(num_p, requests, &current, MPI_STATUS_IGNORE);
+			for (j = 0; j < 16; ++j)
+			{
+				//printf("%.10lf\n", allocation[current][j]);
+			}
+			printf("\n");
 			Insert_Array_In_Matrix(
 				root_field, 
 				m, 
@@ -236,9 +397,6 @@ int main(int argc, char **argv)
 		Write_Matrix(root_field, m, n);
 	}
 
-	// write neighbar_com function
-
-	//MPI_Gather()
 	MPI_Finalize();
 
 	return 0;