SColson82 · June 23, 2025 20:23 · Jun 23, 2025
diff --git a/monte_carlo_pi_hybrid.f90 b/monte_carlo_pi_hybrid.f90
@@ -0,0 +1,88 @@
+! https://stackoverflow.com/a/24077338
+module str2int_mod
+contains
+
+  elemental subroutine str2int(str,int,stat)
+    implicit none
+    ! Arguments
+    character(len=*),intent(in) :: str
+    integer*8,intent(out)       :: int
+    integer,intent(out)         :: stat
+
+    read(str,*,iostat=stat)  int
+  end subroutine str2int
+
+end module
+
+program monte_carlo_pi_hybrid
+  use omp_lib
+  use mpi
+  use str2int_mod
+
+  implicit none
+  real :: x, y, r
+  integer*8 :: i, close_enough = 0, total
+  integer*8 :: total_close_enough, normal_bin_size, bigger_bin_size, bigger_bin_count, my_size
+  integer :: ierr
+! integer can go up to                                          2,147,483,647
+! integer*8 can go up to                            9,223,372,036,854,775,807
+! integer*16 can go up to 170,141,183,460,469,231,731,687,303,715,884,105,727
+! if total = 1000000000 takes 18 seconds on 1 core, the maximum integer*8
+! would take 9.2*10^9 times longer (around 5200 years), and the maximum
+! integer*16 would take 1.7*10^29 times longer (around 9.7*10^22 years)
+  integer :: comm, rank, world_size
+  character(len=80) str_total
+  double precision :: time_start, time_end
+
+  call get_command_argument(1, str_total)
+  call str2int(str_total, total, ierr)
+
+  call random_seed()
+
+  call mpi_init(ierr)
+  call mpi_comm_size(MPI_COMM_WORLD, world_size, ierr)
+  call mpi_comm_rank(MPI_COMM_WORLD, rank, ierr)
+
+  if (rank==0) then
+    print *, 'total is', total, ', world_size is', world_size, 'and total % world_size is', modulo(total, world_size)
+    time_start = omp_get_wtime()
+  end if
+
+  normal_bin_size = (total / world_size)
+  bigger_bin_size = normal_bin_size + 1
+  bigger_bin_count = modulo(total, world_size)
+  if (rank < bigger_bin_count) then
+    my_size = bigger_bin_size
+  else
+    my_size = normal_bin_size
+  end if
+
+!$omp parallel
+  print *, "I am thread number", omp_get_thread_num(), "in rank", rank
+!$omp end parallel
+
+!$omp parallel do private(i, x, y, r) shared(total) reduction(+:close_enough)
+  do i=1, my_size
+    call random_number(x)
+    call random_number(y)
+    r = sqrt(x*x + y*y)
+    if (r<=1.0) then
+      close_enough = close_enough + 1
+    end if
+  end do
+
+
+  call MPI_Reduce(close_enough, total_close_enough, 1, mpi_integer8, mpi_sum, 0, MPI_COMM_WORLD, ierr)  
+
+  call mpi_finalize(ierr)
+  print *, "I am rank", rank, ", my number of close enoughs is", close_enough, "my_size is", my_size
+
+  if (rank == 0) then
+    print *, "Total close enoughs: ", total_close_enough, "total:", total
+    print *, "Pi is approximately", 4.0*total_close_enough/total
+    time_end = omp_get_wtime()
+    print *, "time,", time_end-time_start, ", nodes, ", world_size
+  end if
+
+
+end program monte_carlo_pi_hybrid