Exercise on user defined datatypes

datatypes/LICENSE.txt

+
+Copyright (C) 2018  CSC - IT Center for Science Ltd.
+
+Licensed under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+Code is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+Copy of the GNU General Public License can be obtained from
+<http://www.gnu.org/licenses/>.
+

datatypes/README.md

+## Using custom datatypes 
+
+Write a program that sends the highlighted elements of a 2D array
+using user defined datatypes from one MPI task to another. Note the
+different assignments for C and Fortran, and remember that C stores
+arrays in a row-major order and Fortran in a column-major order. You can
+start from skeleton codes in [C](./c) or [Fortran](./fortran)
+
+a) 
+
+![](img/vector.png)
+
+b)
+
+![](img/indexed.png)
+
+c) Write a program that sends an array of structures (derived types in
+Fortran) between two tasks. Implement a user-defined datatype that can
+be used for sending the structured data and verify that the
+communication is performed successfully. Check the size and true
+extent of your type. A skeleton code is provided in
+[c/struct_type.c](c/struct_type.c) or
+[fortran/struct_type.F90](fortran/struct_type.F90).
+
+d) Implement the sending of structured data also by sending just a
+stream of bytes (type `MPI_BYTE`). Verify correctness and compare the
+performance of these two approaches.
+

datatypes/c/custom_type_a-b.c

+#include <stdio.h>
+#include <mpi.h>
+
+int main(int argc, char **argv)
+{
+    int rank;
+    int array[8][8];
+    //TODO: Declare a variable storing the MPI datatype
+
+    int i, j;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    // Initialize arrays
+    if (rank == 0) {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = (i + 1) * 10 + j + 1;
+            }
+        }
+    } else {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = 0;
+            }
+        }
+    }
+
+    if (rank == 0) {
+        printf("Data in rank 0\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    //TODO: Create datatype that describes one column. Use MPI_Type_vector.
+
+    //TODO: Send first column of matrix form rank 0 to rank 1
+
+    //TODO: free datatype
+
+    // Print out the result on rank 1
+    // The application is correct if the first column has the values of rank 0
+    if (rank == 1) {
+        printf("Received data\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    MPI_Finalize();
+
+    return 0;
+}

datatypes/c/solution/custom_type_a.c

+#include <stdio.h>
+#include <mpi.h>
+
+int main(int argc, char **argv)
+{
+    int rank;
+    int array[8][8];
+    MPI_Datatype columntype;
+
+    int i, j;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    // Initialize arrays
+    if (rank == 0) {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = (i + 1) * 10 + j + 1;
+            }
+        }
+    } else {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = 0;
+            }
+        }
+    }
+
+    if (rank == 0) {
+        printf("Data in rank 0\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    //TODO: Create datatype that describes one column. Use MPI_Type_vector.
+
+    // Create datatype
+    MPI_Type_vector(8, 1, 8, MPI_INT, &columntype);
+    MPI_Type_commit(&columntype);
+
+    // Send first column of matrix
+    if (rank == 0) {
+        MPI_Send(&array[0][1], 1, columntype, 1, 1, MPI_COMM_WORLD);
+    } else if (rank == 1) {
+        MPI_Recv(&array[0][1], 1, columntype, 0, 1, MPI_COMM_WORLD,
+                 MPI_STATUS_IGNORE);
+    }
+
+    // Print out the result
+    if (rank == 1) {
+        printf("Received data\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    MPI_Type_free(&columntype);
+    MPI_Finalize();
+
+    return 0;
+}

datatypes/c/solution/custom_type_b.c

+#include <stdio.h>
+#include <mpi.h>
+
+int main(int argc, char **argv)
+{
+    int rank;
+    int array[8][8];
+    MPI_Datatype indexedtype;
+    int displs[4];
+    int counts[4];
+
+    int i, j;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    // Initialize arrays
+    if (rank == 0) {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = (i + 1) * 10 + j + 1;
+            }
+        }
+    } else {
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                array[i][j] = 0;
+            }
+        }
+    }
+
+    if (rank == 0) {
+        printf("Data in rank 0\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    // Create datatype
+    for (i = 0; i < 4; i++) {
+        counts[i] = i + 1;
+        displs[i] = i + 2 * i * 8;
+    }
+
+    MPI_Type_indexed(4, counts, displs, MPI_INT, &indexedtype);
+    MPI_Type_commit(&indexedtype);
+
+    // Send first indexed of matrix
+    if (rank == 0) {
+        MPI_Send(array, 1, indexedtype, 1, 1, MPI_COMM_WORLD);
+    } else if (rank == 1) {
+        MPI_Recv(array, 1, indexedtype, 0, 1, MPI_COMM_WORLD,
+                 MPI_STATUS_IGNORE);
+    }
+
+    // Print out the result on rank 1
+    // The application is correct if the first column has the values of rank 0
+    if (rank == 1) {
+        printf("Received data\n");
+        for (i = 0; i < 8; i++) {
+            for (j = 0; j < 8; j++) {
+                printf("%3d", array[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    MPI_Type_free(&indexedtype);
+    MPI_Finalize();
+
+    return 0;
+}

datatypes/c/solution/struct_type_c.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+
+
+int main(int argc, char *argv[])
+{
+    int n = 1000, cnt = 3, reps = 10000;
+    particle particles[n];
+    int i, j, myid, ntasks, blocklen[cnt];
+    MPI_Datatype particletype, temptype, types[cnt];
+    MPI_Aint disp[cnt], dist[2], lb, extent;
+    double t1, t2;
+
+    typedef struct {
+        float coords[3];
+        int charge;
+        char label[2];
+    } particle;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+    /* fill in some values for the particles */
+    if (myid == 0) {
+        for (i = 0; i < 1000; i++) {
+            for (j = 0; j < 3; j++) {
+                particles[i].coords[j] = (float)rand() / (float)RAND_MAX * 10.0;
+            }
+            particles[i].charge = 54;
+            strcpy(particles[i].label, "Xe");
+        }
+    }
+    /* TODO: define the datatype for the struct particle */
+    types[0] = MPI_FLOAT;
+    types[1] = MPI_INT;
+    types[2] = MPI_CHAR;
+    blocklen[0] = 3;
+    blocklen[1] = 1;
+    blocklen[2] = 2;
+    MPI_Get_address(&particles[0].coords, &disp[0]);
+    MPI_Get_address(&particles[0].charge, &disp[1]);
+    MPI_Get_address(&particles[0].label, &disp[2]);
+    disp[2] -= disp[0];
+    disp[1] -= disp[0];
+    disp[0] = 0;
+    MPI_Type_create_struct(cnt, blocklen, disp, types, &particletype);
+    MPI_Type_commit(&particletype);
+
+    /* TODO: check extent (not really necessary on most platforms) */
+    MPI_Type_get_extent(particletype, &lb, &extent);
+    MPI_Get_address(&particles[0], &dist[0]);
+    MPI_Get_address(&particles[1], &dist[1]);
+    if (extent != (dist[1] - dist[0])) {
+        temptype = particletype;
+        lb = 0;
+        extent = disp[1] - disp[0];
+        MPI_Type_create_resized(temptype, lb, extent, &particletype);
+        MPI_Type_commit(&particletype);
+        MPI_Type_free(&temptype);
+    }
+
+    /* communicate using the created particletype */
+    t1 = MPI_Wtime();
+    if (myid == 0) {
+        for (i = 0; i < reps; i++) {
+            MPI_Send(particles, n, particletype, 1, i, MPI_COMM_WORLD);
+        }
+    } else if (myid == 1) {
+        for (i = 0; i < reps; i++) {
+            MPI_Recv(particles, n, particletype, 0, i, MPI_COMM_WORLD,
+                     MPI_STATUS_IGNORE);
+	}
+    }
+    t2 = MPI_Wtime();
+
+    printf("Time: %i, %e \n", myid, (t2 - t1) / (double)reps);
+    printf("Check: %i: %s %f %f %f \n", myid, particles[n - 1].label,
+           particles[n - 1].coords[0], particles[n - 1].coords[1],
+           particles[n - 1].coords[2]);
+
+    MPI_Type_free(&particletype);
+    MPI_Finalize();
+    return 0;
+}

datatypes/c/solution/struct_type_d.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+
+
+int main(int argc, char *argv[])
+{
+    int n = 1000, cnt = 3, reps = 10000;
+    particle particles[n];
+    int i, j, myid, ntasks;
+    MPI_Aint lb0, lb1, extent;
+    double t1, t2;
+
+    typedef struct {
+        float coords[3];
+        int charge;
+        char label[2];
+    } particle;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+    /* fill in some values for the particles */
+    if (myid == 0) {
+        for (i = 0; i < 1000; i++) {
+            for (j = 0; j < 3; j++) {
+                particles[i].coords[j] = (float)rand() / (float)RAND_MAX * 10.0;
+            }
+            particles[i].charge = 54;
+            strcpy(particles[i].label, "Xe");
+        }
+    }
+    /* TODO: determine the true extent of one particle struct */
+    MPI_Get_address(&particles[0], &lb0);
+    MPI_Get_address(&particles[1], &lb1);
+    extent = lb1 - lb0;
+
+    /* TODO: send and receive using the MPI_BYTE datatype */
+    t1 = MPI_Wtime();
+    if (myid == 0) {
+        for (i = 0; i < reps; i++) {
+            MPI_Send(particles, n * extent, MPI_BYTE, 1, i, MPI_COMM_WORLD);
+        }
+    } else if (myid == 1) {
+        for (i = 0; i < reps; i++) {
+            MPI_Recv(particles, n * extent, MPI_BYTE, 0, i, MPI_COMM_WORLD,
+                     MPI_STATUS_IGNORE);
+	}
+    }
+    t2 = MPI_Wtime();
+
+    printf("Time: %i, %e \n", myid, (t2 - t1) / (double)reps);
+    printf("Check: %i: %s %f %f %f \n", myid, particles[n - 1].label,
+           particles[n - 1].coords[0], particles[n - 1].coords[1],
+           particles[n - 1].coords[2]);
+
+    MPI_Finalize();
+    return 0;
+}

datatypes/c/struct_type.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+
+
+int main(int argc, char *argv[])
+{
+    int n = 1000, cnt = 3, reps = 10000;
+    particle particles[n];
+    int i, j, myid, ntasks, blocklen[cnt];
+    MPI_Datatype particletype, temptype;
+    MPI_Aint disp[cnt], dist[2], lb, extent;
+    double t1, t2;
+
+    typedef struct {
+        float coords[3];
+        int charge;
+        char label[2];
+    } particle;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+    /* fill in some values for the particles */
+    if (myid == 0) {
+        for (i = 0; i < 1000; i++) {
+            for (j = 0; j < 3; j++) {
+                particles[i].coords[j] = (float)rand() / (float)RAND_MAX * 10.0;
+            }
+            particles[i].charge = 54;
+            strcpy(particles[i].label, "Xe");
+        }
+    }
+    /* TODO (c): define the datatype for the struct particle  using MPI_Type_create_struct
+       You can use MPI_Get_address to compute offsets.
+    */
+
+    /* TODO (c): check extent (not really necessary on most platforms) That is,
+     * check that extent is identical to the distance between two consequtive
+     * structs in an array
+     * Tip, use MPI_Type_get_extent and  MPI_Get_address
+     */
+
+    if (extent != (dist[1] - dist[0])) {
+        /*TODO (c), resize particle type to correct extent */
+    }
+
+    /* communicate using the created particletype */
+    t1 = MPI_Wtime();
+    if (myid == 0) {
+        for (i = 0; i < reps; i++) {
+            MPI_Send(particles, n, particletype, 1, i, MPI_COMM_WORLD);
+        }
+    } else if (myid == 1) {
+        for (i = 0; i < reps; i++) {
+            MPI_Recv(particles, n, particletype, 0, i, MPI_COMM_WORLD,
+                     MPI_STATUS_IGNORE);
+	}
+    }
+    t2 = MPI_Wtime();
+
+    printf("Time: %i, %e \n", myid, (t2 - t1) / (double)reps);
+    printf("Check: %i: %s %f %f %f \n", myid, particles[n - 1].label,
+           particles[n - 1].coords[0], particles[n - 1].coords[1],
+           particles[n - 1].coords[2]);
+
+    //TODO: Free datatype
+
+    MPI_Finalize();
+    return 0;
+}

datatypes/fortran/custom_type_a-b.F90

+program datatype1
+  use mpi
+  implicit none
+
+  integer, dimension(8,8) :: array
+  integer :: rank, ierr
+  !TODO: declare variable for datatype
+  integer :: i, j
+
+  call mpi_init(ierr)
+  call mpi_comm_rank(MPI_COMM_WORLD, rank ,ierr)
+
+  ! initialize arrays
+  if (rank == 0) then
+     do i=1,8
+        do j=1,8
+           array(i,j) = i*10 + j
+        end do
+     end do
+  else
+     array(:,:) = 0
+  end if
+
+  if (rank == 0) then
+     write(*,*) 'Data in rank 0'
+     do i=1,8
+        write(*,'(8I3)') array(i, :)
+     end do
+  end if
+
+
+  !TODO: create datatype describing one row, use mpi_type_vector
+
+  !TODO: send first row of matrix from rank 0 to 1
+
+  ! Print out the result
+  if (rank == 1) then
+     write(*,*) 'Received data'
+     do i=1,8
+        write(*,'(8I3)') array(i, :)
+     end do
+  end if
+
+  !TODO free datatype
+
+  call mpi_finalize(ierr)
+
+end program datatype1

datatypes/fortran/solution/custom_type_a.F90

+program datatype1
+  use mpi
+  implicit none
+
+  integer, dimension(8,8) :: array
+  integer :: rank, ierr
+  integer :: rowtype
+  integer :: i, j
+
+  call mpi_init(ierr)
+  call mpi_comm_rank(MPI_COMM_WORLD, rank ,ierr)
+
+  ! initialize arrays
+  if (rank == 0) then
+     do i=1,8
+        do j=1,8
+           array(i,j) = i*10 + j
+        end do
+     end do
+  else
+     array(:,:) = 0
+  end if
+
+  if (rank == 0) then
+     write(*,*) 'Data in rank 0'
+     do i=1,8
+        write(*,'(8I3)') array(i, :)
+     end do
+  end if
+
+  ! create datatype
+  call mpi_type_vector(8, 1, 8, MPI_INTEGER, rowtype, ierr)
+  call mpi_type_commit(rowtype, ierr)
+
+  ! send first row of matrix
+  if (rank == 0) then
+     call mpi_send(array(2, 1), 1, rowtype, 1, 1, MPI_COMM_WORLD, ierr)
+  else if (rank == 1) then
+     call mpi_recv(array(2, 1), 1, rowtype, 0, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, &
+          ierr)
+  end if
+
+  ! Print out the result
+  if (rank == 1) then
+     write(*,*) 'Received data'
+     do i=1,8
+        write(*,'(8I3)') array(i, :)
+     end do
+  end if
+
+  call mpi_type_free(rowtype, ierr)
+  call mpi_finalize(ierr)
+
+end program datatype1

datatypes/fortran/solution/custom_type_b.F90

+program datatype1
+  use mpi
+  implicit none
+
+  integer, dimension(8,8) :: array
+  integer :: rank, ierr
+  integer :: indexedtype
+  integer, dimension(4) :: counts, displs
+  integer :: i, j
+
+  call mpi_init(ierr)
+  call mpi_comm_rank(MPI_COMM_WORLD, rank ,ierr)
+
+  ! initialize arrays
+  if (rank == 0) then
+     do i=1,8
+        do j=1,8
+           array(i,j) = i*10 + j
+        end do
+     end do
+  else
+     array(:,:) = 0
+