Created
June 23, 2025 20:23
-
-
Save SColson82/72d1e02c6c727b3cafee4931643a5a9a to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ! 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment