MPI

Apps/scone.f90

@@ -1,8 +1,9 @@
 program scone
 
   use numPrecision
-  use genericProcedures,          only : printStart
+  use display_func,               only : printStart, statusMsg
   use openmp_func,                only : ompSetNumThreads
+  use mpi_func,                   only : mpiInit, mpiFinalise
   use commandLineUI,              only : getInputFile, clOptionIsPresent, addClOption, getFromCL
   use dictionary_class,           only : dictionary
   use dictParser_func,            only : fileToDict
@@ -12,6 +13,7 @@ program scone
   use timer_mod                 , only : registerTimer, timerStart, timerStop, timerTime, secToChar
 
   implicit none
+
   type(dictionary)                  :: input
   class(physicsPackage),allocatable :: core
   character(:),allocatable          :: inputPath
@@ -29,6 +31,9 @@ program scone
   ! Get path to input file
   call getInputFile(inputPath)
 
+  ! Initialize MPI
+  call mpiInit()
+
   ! Set Number of threads
   if (clOptionIsPresent('--omp')) then
     call getFromCL(cores, '--omp', 1)
@@ -56,6 +61,10 @@ program scone
   call core % run()
 
   call timerStop(timerIdx)
-  print *, 'Total calculation time: ', trim(secToChar(timerTime(timerIdx)))
-  print *, 'Have a good day and enjoy your result analysis!'
+
+  call mpiFinalise()
+
+  call statusMsg('Total calculation time: ' // trim(secToChar(timerTime(timerIdx))))
+  call statusMsg('Have a good day and enjoy your result analysis!')
+
 end program scone

CMakeLists.txt

@@ -19,6 +19,8 @@ option(LTO         "Enables link-time optimisation"     ON)
 option(COVERAGE    "Collect Coverage Info"              OFF)
 option(DEBUG       "Enable extra run-time checks"       OFF)
 option(OPENMP      "Enable parallelism with OpenMP"     ON)
+option(MPI         "Enable parallelism with MPI"        ON)
+
 
 # Include local cmake modules. TODO: Test on WINDOWS!!!
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake")
@@ -77,6 +79,12 @@ if(OPENMP)
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}")
 endif()
 
+# Set up MPI
+if(MPI)
+  find_package(MPI COMPONENTS Fortran REQUIRED)
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -DMPI")
+endif()
+
 
 ####################################################################################################
 # CHECK FOR DEPENDENCIES
@@ -137,7 +145,10 @@ get_property(SRCS GLOBAL PROPERTY SRCS_LIST)
 # Compile library
 add_library(scone STATIC ${SRCS})
 target_compile_options(scone PRIVATE ${scone_extra_flags} )
-target_link_libraries(scone  ${LAPACK_LIBRARIES} )
+target_link_libraries(scone PUBLIC ${LAPACK_LIBRARIES} )
+if(MPI)
+  target_link_libraries(scone PUBLIC MPI::MPI_Fortran)
+endif()
 
 if(LTO)
   set_property(TARGET scone PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE)
@@ -163,7 +174,8 @@ if(BUILD_TESTS)
 
   add_pfunit_ctest(unitTests
     TEST_SOURCES ${UNIT_TESTS_RELATIVE}
-    LINK_LIBRARIES scone ${LAPACK_LIBRARIES}
+    LINK_LIBRARIES scone
+    MAX_PES 1
     )
 
   # pFUnit may have a bug which causes a unused variable `class(Test), allocatable :: t` be
@@ -184,7 +196,7 @@ if(BUILD_TESTS)
 
   add_pfunit_ctest(integrationTests
     TEST_SOURCES ${INTEGRATION_TESTS_RELATIVE}
-    LINK_LIBRARIES scone ${LAPACK_LIBRARIES}
+    LINK_LIBRARIES scone
     WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
     )
 
@@ -206,6 +218,7 @@ message(STATUS "  Link-time optimisation:     " ${LTO})
 message(STATUS "  Code coverage logging:      " ${COVERAGE})
 message(STATUS "  Extra runtime debug checks: " ${DEBUG})
 message(STATUS "  OpenMP parallelism:         " ${OPENMP})
+message(STATUS "  MPI parallelism:            " ${MPI})
 message(STATUS "  Fortran compiler:     " ${CMAKE_Fortran_COMPILER})
 message(STATUS "  Compiler version:     " ${CMAKE_Fortran_COMPILER_VERSION})
 message(STATUS "  OpenMP version:       " ${OpenMP_CXX_VERSION})

DataStructures/CMakeLists.txt

@@ -4,13 +4,15 @@ add_sources( ./dictionary_class.f90
              ./charMap_class.f90
              ./stack_class.f90
              ./dynArray_class.f90
-             ./dictParser_func.f90)
+             ./dictParser_func.f90
+             ./heapQueue_class.f90)
 
 # Add Unit Tests to a global List
 add_unit_tests( ./Tests/dictionary_test.f90
                 ./Tests/intMap_test.f90
                 ./Tests/charMap_test.f90
                 ./Tests/dynArray_test.f90
-                ./Tests/dictParser_test.f90)
+                ./Tests/dictParser_test.f90
+                ./Tests/heapQueue_test.f90)
 
 add_integration_tests( ./Tests/dictParser_iTest.f90)

DataStructures/Tests/heapQueue_test.f90

@@ -0,0 +1,64 @@
+module heapQueue_test
+  use numPrecision
+  use heapQueue_class, only: heapQueue
+  use funit
+
+  implicit none
+
+contains
+
+!!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+!! PROPER TESTS BEGIN HERE
+!!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+
+  !!
+  !! Test with simple sequence without reaching maximum size
+  !!
+@Test
+  subroutine testBelowMaximum()
+    type(heapQueue)   :: hq
+    integer(shortInt) :: i
+    real(defReal), dimension(*), parameter :: seq = [2.0_defReal, 3.0_defReal, 1.0_defReal, 4.0_defReal, 5.0_defReal]
+
+    call hq % init(8)
+
+    @assertEqual(hq % getSize(), 0)
+
+    do i = 1, size(seq)
+      call hq % pushReplace(seq(i))
+    end do
+
+    ! Check that the maximum value is the maximum value in the sequence
+    @assertEqual(hq % maxValue(), maxval(seq))
+    @assertEqual(hq % getSize(), size(seq))
+
+  end subroutine testBelowMaximum
+
+  !!
+  !! Test the intended use case
+  !!
+  !! It is not a very good test but I had no better idea at the moment [MAK]
+  !!
+  @Test
+  subroutine testAboveMaximum()
+    type(heapQueue)   :: hq
+    integer(shortInt) :: i
+    real(defReal)     :: val
+    real(defReal), dimension(*), parameter :: seq = [2.0_defReal, 3.0_defReal, 1.0_defReal, 1.4_defReal, 5.0_defReal]
+
+    call hq % init(3)
+
+    ! Push upper bound
+    call hq % pushReplace(1000.0_defReal)
+
+    do i = 1, size(seq)
+      val = seq(i)
+      if (val < hq % maxValue()) call hq % pushReplace(seq(i))
+    end do
+
+    ! Check that the threshold is correct
+    @assertEqual(hq % maxValue(), 2.0_defReal)
+
+  end subroutine testAboveMaximum
+
+end module heapQueue_test

DataStructures/heapQueue_class.f90

@@ -0,0 +1,185 @@
+module heapQueue_class
+
+  use numPrecision
+  use genericProcedures, only : swap
+  use errors_mod,        only : fatalError
+
+  implicit none
+  private
+
+  !!
+  !! A fixed size heap queue designed to store N smallest values of a sample
+  !!
+  !! The queue is implemented as a binary heap, that is a binary tree with the
+  !! heap property. The heap property is that the value of a parent node is
+  !! larger than the value of both its children. This means that the largest
+  !! value is always at the root of the tree.
+  !!
+  !! The queue is implemented as an array with the root at index 1. The children
+  !! of a node at index i are at 2*i and 2*i + 1. Note that for the even number
+  !! of elements the last node will only have one child.
+  !!
+  !! This data structure is intended to be used for sampling without replacement
+  !! as it allows to find a value of a threshold that selects K values from a
+  !! stream of the N random numbers.
+  !!
+  !! Interface:
+  !!   init        -> Initialise the queue to a given size
+  !!   pushReplace -> Add a value to the queue either growing the size or replacing the largest element
+  !!                  if maximum size was reached
+  !!   maxValue    -> Returns the largest value in the queue
+  !!   getSize     -> Returns the current size of the queue
+  !!
+  type, public :: heapQueue
+    private
+    real(defReal), dimension(:) , allocatable :: heap
+    integer(shortInt)                         :: size
+  contains
+    procedure :: init
+    procedure :: pushReplace
+    procedure :: maxValue
+    procedure :: getSize
+
+    procedure, private :: push
+    procedure, private :: replace
+  end type heapQueue
+
+contains
+
+  !!
+  !! Initialise the queue to a given size
+  !!
+  !! Args:
+  !!  maxSize [in] -> Maximum size of the queue
+  !!
+  subroutine init(self, maxSize)
+    class(heapQueue), intent(out) :: self
+    integer(shortInt), intent(in) :: maxSize
+
+    self % size = 0
+    allocate(self % heap(maxSize))
+
+  end subroutine init
+
+  !!
+  !! Add a value to to queue either growing the size or replacing the largest
+  !!
+  !! Args:
+  !!  value [in] -> Value to add to the queue
+  !!
+  subroutine pushReplace(self, val)
+    class(heapQueue), intent(inout) :: self
+    real(defReal), intent(in)       :: val
+
+    if (self % size < size(self % heap)) then
+      call self % push(val)
+    else
+      call self % replace(val)
+    end if
+
+  end subroutine pushReplace
+
+  !!
+  !! Add a value to the queue
+  !!
+  !! Assumes enough space is available
+  !!
+  !! Args:
+  !!  value [in] -> Value to add to the queue
+  !!
+  subroutine push(self, val)
+    class(heapQueue), intent(inout) :: self
+    real(defReal), intent(in)       :: val
+    integer(shortInt)               :: parent, child
+
+    ! Increase the size of the queue and add the new value
+    self % size = self % size + 1
+    self % heap(self % size) = val
+
+    ! Shift the new value up the heap to restore the heap property
+    child = self % size
+    parent = child / 2
+
+    do while (child > 1 .and. self % heap(parent) < self % heap(child))
+      call swap(self % heap(parent), self % heap(child))
+      child = parent
+      parent = child / 2
+    end do
+
+  end subroutine push
+
+  !!
+  !! Replaces the largest value in the queue with a new value
+  !!
+  !! Args:
+  !!  value [in] -> Value to add to the queue
+  !!
+  subroutine replace(self, val)
+    class(heapQueue), intent(inout) :: self
+    real(defReal), intent(in)       :: val
+    integer(shortInt)               :: parent, child
+
+    self % heap(1) = val
+
+    parent = 1
+    child = 2
+
+    ! Shift down the new value until heap property is restored be comparing
+    ! with the largest child and swapping if necessary
+    do while (child <= self % size)
+
+      ! We need to consider the case where there is only one child in a node
+      ! (can happen when the size is even). If child is the last element it is the
+      ! node with no sibling
+      if (child /= self % size .and. self % heap(child) < self % heap(child + 1)) then
+        child = child + 1
+      end if
+
+      ! If the parent is larger than the larger child we are done
+      if (self % heap(parent) >= self % heap(child)) then
+        return
+      end if
+
+      ! Otherwise swap the parent with the larger child and continue
+      ! the recursion
+      call swap(self % heap(parent), self % heap(child))
+      parent = child
+
+      ! Child points to the start of next level
+      child = parent * 2
+    end do
+
+  end subroutine replace
+
+  !!
+  !! Returns the largest value in the queue
+  !!
+  !! Errors:
+  !!   fatal error if the queue is empty
+  !!
+  function maxValue(self) result(val)
+    class(heapQueue), intent(in) :: self
+    real(defReal)                :: val
+    character(100), parameter :: Here = "maxValue (heapQueue_class.f90)"
+
+    if (self % size == 0) then
+      call fatalError(Here, "The queue is empty!")
+    end if
+
+    val = self % heap(1)
+
+  end function maxValue
+
+  !!
+  !! Get the current size of the queue
+  !!
+  pure function getSize(self) result(size)
+    class(heapQueue), intent(in) :: self
+    integer(shortInt)            :: size
+
+    size = self % size
+
+  end function getSize
+
+
+end module heapQueue_class