Add derived type and I/O lab

lab3/README.md

@@ -13,6 +13,7 @@ Three hours
 # Source Codes
 
 - MPI I/O. Serial hello world in C and Fortran ([hello_mpi.c](hello_mpi.c) and [hello_mpi.f90](hello_mpi.f90))
+- MPI Derived types and I/O. Serial STL file reader in C and Fortran ([mpi_derived_types.c](mpi_derived_types.c) and [mpi_derived_types.f90](mpi_derived_types.f90)
 - MPI Latency: C and Fortran ([mpi_latency.c](mpi_latency.c) and [mpi_latency.f90](mpi_latency.f90))
 - MPI Bandwidth : C and Fortran ([mpi_bandwidth.c](mpi_bandwidth.c) and [mpi_bandwidth.f90](mpi_bandwidth.f90))
 - MPI Bandwidth Non-Blocking: C and Fortran ([mpi_bandwidth-nonblock.c](mpi_bandwidth-nonblock.c) 
@@ -29,7 +30,19 @@ MPI I/O is used so that results can be written to the same file in parallel. Tak
 
 The simplest solution is likely to be for you to create a character buffer, and then use the MPI_File_write_at function.
 
-# Exercises 2 - Bandwidth and latency between nodes
+# Exercise 2 - MPI I/O and dervied types
+
+Take the serial stl reader and modify it such that the stl file is read and written in parallel using collective MPI I/O.
+
+The simplest solution is likely to create a derived type for each triangle, and then use the MPI_File_write_at_all function. A correct solution will have the same MD5 hash for both stl models (input and output)
+
+```
+md5sum out.stl data/sphere.stl
+822aba6dc20cc0421f92ad50df95464c  out.stl
+822aba6dc20cc0421f92ad50df95464c  data/sphere.stl
+```
+
+# Exercises 3 - Bandwidth and latency between nodes
 
 Use `mpi_wtime` to compute latency and bandwidth with the bandwidth and latency codes above
 

lab3/mpi_derived_types.c

+/*
+  STL file format 
+
+  UINT8[80] – Header
+  UINT32 – Number of triangles
+
+  foreach triangle
+  REAL32[3] – Normal vector
+  REAL32[3] – Vertex 1
+  REAL32[3] – Vertex 2
+  REAL32[3] – Vertex 3
+  UINT16 – Attribute byte count
+  end
+
+  (see https://en.wikipedia.org/wiki/STL_(file_format)
+  
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <string.h>
+#include <mpi.h>
+
+#define STL_HDR_SIZE 80
+
+typedef struct {
+  float n[3];			/* Normal vector */
+  float v1[3];			/* Vertex 1 */
+  float v2[3];			/* Vertex 2 */
+  float v3[3];			/* Vertex 3 */
+  uint16_t attrib;		/* Attribute byte count */
+} __attribute__((packed)) stl_triangle_t;
+
+typedef struct {
+  char hdr[STL_HDR_SIZE];	/* Header */
+  uint32_t n_tri;		/* Number of triangles */
+  stl_triangle_t *tri;		/* Triangles */
+} stl_model_t;
+
+
+void stl_read(const char *fname, stl_model_t *model) {
+  FILE *fp;
+  int pe_size, pe_rank;
+
+  MPI_Comm_size(MPI_COMM_WORLD, &pe_size);
+  MPI_Comm_rank(MPI_COMM_WORLD, &pe_rank);  
+
+  fp = fopen(fname, "r");
+
+  if (pe_rank == 0) printf("Reading STL file: %s\n", fname);
+
+  /* Read STL header */
+  fread(model->hdr, sizeof(char), STL_HDR_SIZE, fp);
+
+  /* Make sure it's a binary STL file */
+  if (strncmp(model->hdr, "solid", 5) == 0) {
+    fprintf(stderr, "ASCII STL files not supported!\n");
+    exit(-1);
+  }
+
+  /* Read how many triangles the file contains */
+  fread(&model->n_tri, sizeof(uint32_t), 1, fp);
+  if (pe_rank == 0) printf("Found: %d triangles\n", model->n_tri);
+
+  /* Allocate memory for triangles, and read them */
+  model->tri = malloc(model->n_tri * sizeof(stl_triangle_t));
+  fread(model->tri, sizeof(stl_triangle_t), model->n_tri, fp);
+
+  fclose(fp);
+  if (pe_rank == 0) printf("Done\n");
+
+}
+
+void stl_write(const char *fname, stl_model_t *model) {
+  FILE *fp;
+  int pe_size, pe_rank;
+
+  MPI_Comm_size(MPI_COMM_WORLD, &pe_size);
+  MPI_Comm_rank(MPI_COMM_WORLD, &pe_rank);
+  
+  fp = fopen(fname, "w");
+  if (pe_rank == 0) printf("Writing STL file: %s\n", fname);
+
+  /* Write STL header */
+  fwrite(model->hdr, sizeof(char), STL_HDR_SIZE, fp);
+
+  /* Write number of triangles */
+  fwrite(&model->n_tri, sizeof(uint32_t), 1, fp);
+
+  /* Write all triangles */
+  fwrite(model->tri, sizeof(stl_triangle_t), model->n_tri, fp);
+
+  fclose(fp);
+  if (pe_rank == 0) printf("Done\n");
+}
+
+int main(int argc, char **argv) {
+  stl_model_t model;
+  
+  MPI_Init(&argc, &argv);
+  
+  stl_read("./data/sphere.stl", &model);
+  stl_write("out.stl", &model);
+  free(model.tri);
+
+  MPI_Finalize();
+  
+  return 0;
+}
+
+

lab3/mpi_derived_types.f90

+!
+!  STL file format 
+!
+!  UINT8[80] – Header
+!  UINT32 – Number of triangles
+!
+!  foreach triangle
+!  REAL32[3] – Normal vector
+!  REAL32[3] – Vertex 1
+!  REAL32[3] – Vertex 2
+!  REAL32[3] – Vertex 3
+!  UINT16 – Attribute byte count
+!  end
+!
+!  (see https://en.wikipedia.org/wiki/STL_(file_format)
+!
+module types
+  integer, parameter :: dp = kind(0.0d0)
+  integer, parameter :: sp = kind(0.0)
+end module types
+
+module stl
+  use mpi
+  use types
+  implicit none
+
+  type :: stl_triangle_t
+     real(kind=sp) :: n(3)
+     real(kind=sp) :: v1(3)
+     real(kind=sp) :: v2(3)
+     real(kind=sp) :: v3(3)
+     integer(kind=2) :: attrib
+  end type stl_triangle_t
+
+  type :: stl_model_t
+     character(len=80) :: hdr
+     integer :: n_tri
+     type(stl_triangle_t), allocatable :: tri(:)
+  end type stl_model_t
+  
+contains
+
+  subroutine stl_read(fname, model)
+    character(len=*), intent(in) :: fname
+    type(stl_model_t), intent(inout) :: model
+    integer :: pe_size, pe_rank, ierr
+
+    call MPI_Comm_size(MPI_COMM_WORLD, pe_size, ierr)
+    call MPI_Comm_rank(MPI_COMM_WORLD, pe_rank, ierr)
+
+    open(42, file=fname, access='stream', form='unformatted')
+    
+    if (pe_rank .eq. 0) write(*,*) 'Reading STL file: ', fname
+    ! Read STL header
+    read(42, pos=1) model%hdr
+    
+    ! Make sure it's a binary STL file 
+    if (model%hdr(1:6) .eq. 'solid') then
+       write(*,*) 'ASCII STL files not supported!'
+       stop 
+    end if
+    
+    ! Read how many triangles the file contains 
+    read(42, pos=81) model%n_tri
+    if (pe_rank .eq. 0) write(*,*) 'Found: ', model%n_tri,' triangles'
+
+    ! Allocate memory for triangles, and read them
+    allocate(model%tri(model%n_tri))    
+    read(42, pos=85) model%tri
+
+    close(42)    
+    if (pe_rank .eq. 0) write(*,*) 'Done'
+
+  end subroutine stl_read
+
+  subroutine stl_write(fname, model)
+    character(len=*), intent(in) :: fname
+    type(stl_model_t), intent(in) :: model
+    integer :: pe_size, pe_rank, ierr
+
+    call MPI_Comm_size(MPI_COMM_WORLD, pe_size, ierr)
+    call MPI_Comm_rank(MPI_COMM_WORLD, pe_rank, ierr)
+    
+    open(17, file=fname, access='stream', form='unformatted')
+    
+    if (pe_rank .eq. 0) write(*,*) 'Writing STL file: ', fname
+
+    ! Write STL header
+    write(17, pos=1) model%hdr
+
+    ! Write number of triangles
+    write(17, pos=81) model%n_tri
+    
+    ! Write all triangles
+    write(17, pos=85) model%tri
+
+    close(17)
+    if (pe_rank .eq. 0) write(*,*) 'Done'
+
+  end subroutine stl_write
+end module stl
+
+program mpi_derived_types
+  use stl
+  use mpi
+  implicit none
+  type(stl_model_t) :: model
+  integer :: ierr
+  call MPI_Init(ierr)
+
+  call stl_read("./data/sphere.stl", model)
+  call stl_write("out.stl", model)
+  deallocate(model%tri)
+
+  call MPI_Finalize(ierr)
+
+end program mpi_derived_types