Mario Limonciello
d1aaae2539
In order for pre-commit to be useful, everything needs to meet a common baseline. Signed-off-by: Mario Limonciello (AMD) <superm1@kernel.org>
385 lines
9.0 KiB
Fortran
385 lines
9.0 KiB
Fortran
program main
|
|
|
|
!*****************************************************************************80
|
|
!
|
|
!! MAIN is the main program for INTERVALS.
|
|
!
|
|
! Discussion:
|
|
!
|
|
! INTERVALS uses MPI routines to multiprocess a computational task.
|
|
!
|
|
! We have a function F(X), an interval [XMIN,XMAX],
|
|
! and a value N.
|
|
!
|
|
! We define N equally spaced points in the interval,
|
|
!
|
|
! X(I) = ( ( N - I ) * XMIN
|
|
! + ( I - 1 ) * XMAX )
|
|
! / ( N - 1 )
|
|
!
|
|
! We thus have N-1 subintervals.
|
|
!
|
|
! We assume we have N processors available.
|
|
!
|
|
! Processor 0 is designated the master processor, assigned
|
|
! to estimating the integral of F(X) over the entire
|
|
! interval [ X(1), X(N) ].
|
|
!
|
|
! For I = 1 to N-1, processor I is assigned the subinterval
|
|
!
|
|
! [ X(I), X(I+1) ]
|
|
!
|
|
! and then estimates the integral Q(I) of F(X) over that
|
|
! subinterval.
|
|
!
|
|
! COMMUNICATION:
|
|
!
|
|
! Processor 0 communicates to processor I the endpoints of
|
|
! the interval it is assigned, and the number of sample points
|
|
! to use in that interval.
|
|
!
|
|
! Processor I communicates to processor 0 the computed value of
|
|
! Q(I).
|
|
!
|
|
! Licensing:
|
|
!
|
|
! This code is distributed under the GNU LGPL license.
|
|
!
|
|
! Modified:
|
|
!
|
|
! 17 March 2006
|
|
!
|
|
! Author:
|
|
!
|
|
! John Burkardt
|
|
!
|
|
! Reference:
|
|
!
|
|
! William Gropp, Ewing Lusk, Anthony Skjellum,
|
|
! Using MPI: Portable Parallel Programming with the
|
|
! Message-Passing Interface,
|
|
! Second Edition,
|
|
! MIT Press, 1999,
|
|
! ISBN: 0262571323.
|
|
!
|
|
! Snir, Otto, Huss-Lederman, Walker, Dongarra,
|
|
! MPI - The Complete Reference,
|
|
! Volume 1, The MPI Core,
|
|
! second edition,
|
|
! MIT Press, 1998.
|
|
!
|
|
use mpi
|
|
|
|
real ( kind = 8 ) f
|
|
real ( kind = 8 ) h
|
|
integer ( kind = 4 ) i
|
|
integer ( kind = 4 ) id
|
|
integer ( kind = 4 ) ierr
|
|
integer ( kind = 4 ) m
|
|
integer ( kind = 4 ) n
|
|
real ( kind = 8 ), parameter :: pi = 3.141592653589793238462643D+00
|
|
integer ( kind = 4 ) process
|
|
integer ( kind = 4 ) process_num
|
|
real ( kind = 8 ) q_global
|
|
real ( kind = 8 ) q_local
|
|
integer ( kind = 4 ) received
|
|
integer ( kind = 4 ) source
|
|
integer ( kind = 4 ) status(MPI_Status_size)
|
|
integer ( kind = 4 ) tag
|
|
integer ( kind = 4 ) target
|
|
real ( kind = 8 ) wtime
|
|
real ( kind = 8 ) x
|
|
real ( kind = 8 ) xb(2)
|
|
real ( kind = 8 ) :: x_max = 1.0D+00
|
|
real ( kind = 8 ) :: x_min = 0.0D+00
|
|
!
|
|
! Establish the MPI environment.
|
|
!
|
|
call MPI_Init ( ierr )
|
|
!
|
|
! Get this process's ID.
|
|
!
|
|
call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )
|
|
!
|
|
! Find out how many processes are available.
|
|
!
|
|
call MPI_Comm_size ( MPI_COMM_WORLD, process_num, ierr )
|
|
!
|
|
! Say hello (once), and shut down right away unless we
|
|
! have at least 2 processes available.
|
|
!
|
|
if ( id == 0 ) then
|
|
call timestamp ( )
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS - Master process:'
|
|
write ( *, '(a)' ) ' FORTRAN90 version'
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) ' An MPI example program.'
|
|
write ( *, '(a)' ) ' A quadrature over an interval is done by'
|
|
write ( *, '(a)' ) ' assigning subintervals to processes.'
|
|
write ( *, '(a,i8)' ) ' The number of processes is ', process_num
|
|
|
|
wtime = MPI_Wtime ( )
|
|
|
|
if ( process_num <= 1 ) then
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS - Master process:'
|
|
write ( *, '(a)' ) ' Need at least 2 processes!'
|
|
call MPI_Finalize ( ierr )
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS - Master process:'
|
|
write ( *, '(a)' ) ' Abnormal end of execution.'
|
|
stop
|
|
end if
|
|
|
|
end if
|
|
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a,i8,a)' ) 'Process ', id, ': Active!'
|
|
!
|
|
! Every process could figure out the endpoints of its interval
|
|
! on its own. But we want to demonstrate communication. So we
|
|
! assume that the assignment of processes to intervals is done
|
|
! only by the master process, which then tells each process
|
|
! what job it is to do.
|
|
!
|
|
if ( id == 0 ) then
|
|
|
|
do process = 1, process_num-1
|
|
|
|
xb(1) = ( real ( process_num - process, kind = 8 ) * x_min &
|
|
+ real ( process - 1, kind = 8 ) * x_max ) &
|
|
/ real ( process_num - 1, kind = 8 )
|
|
|
|
xb(2) = ( real ( process_num - process - 1, kind = 8 ) * x_min &
|
|
+ real ( process, kind = 8 ) * x_max ) &
|
|
/ real ( process_num - 1, kind = 8 )
|
|
|
|
target = process
|
|
tag = 1
|
|
|
|
call MPI_Send ( xb, 2, MPI_DOUBLE_PRECISION, target, tag, &
|
|
MPI_COMM_WORLD, ierr )
|
|
|
|
end do
|
|
|
|
else
|
|
|
|
tag = 1
|
|
|
|
call MPI_Recv ( xb, 2, MPI_DOUBLE_PRECISION, 0, tag, &
|
|
MPI_COMM_WORLD, status, ierr )
|
|
|
|
end if
|
|
!
|
|
! Wait here until everyone has gotten their assignment.
|
|
!
|
|
call MPI_Barrier ( MPI_COMM_WORLD, ierr )
|
|
|
|
if ( id == 0 ) then
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS - Master process:'
|
|
write ( *, '(a)' ) ' Subintervals have been assigned.'
|
|
end if
|
|
!
|
|
! Every process needs to be told the number of points to use.
|
|
! Since this is the same value for everybody, we use a broadcast.
|
|
! Again, we are doing it in this roundabout way to emphasize that
|
|
! the choice for M could really be made at runtime, by processor 0,
|
|
! and then sent out to the others.
|
|
!
|
|
m = 100
|
|
|
|
call MPI_Bcast ( m, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
|
|
!
|
|
! Now, every process EXCEPT 0 computes its estimate of the
|
|
! integral over its subinterval, and sends the result back
|
|
! to process 0.
|
|
!
|
|
if ( id /= 0 ) then
|
|
|
|
q_local = 0.0D+00
|
|
|
|
do i = 1, m
|
|
|
|
x = ( real ( 2 * m - 2 * i + 1, kind = 8 ) * xb(1) &
|
|
+ real ( 2 * i - 1, kind = 8 ) * xb(2) ) &
|
|
/ real ( 2 * m, kind = 8 )
|
|
|
|
q_local = q_local + f ( x )
|
|
|
|
end do
|
|
|
|
q_local = q_local * ( xb(2) - xb(1) ) / real ( m, kind = 8 )
|
|
|
|
tag = 2
|
|
|
|
call MPI_Send ( q_local, 1, MPI_DOUBLE_PRECISION, 0, tag, &
|
|
MPI_COMM_WORLD, ierr )
|
|
!
|
|
! Process 0 expects to receive N-1 partial results.
|
|
!
|
|
else
|
|
|
|
received = 0
|
|
q_global = 0.0D+00
|
|
|
|
do while ( received < process_num - 1 )
|
|
|
|
source = MPI_ANY_SOURCE
|
|
tag = 2
|
|
|
|
call MPI_Recv ( q_local, 1, MPI_DOUBLE_PRECISION, source, tag, &
|
|
MPI_COMM_WORLD, status, ierr )
|
|
|
|
q_global = q_global + q_local
|
|
received = received + 1
|
|
|
|
end do
|
|
|
|
end if
|
|
!
|
|
! The master process prints the answer.
|
|
!
|
|
if ( id == 0 ) then
|
|
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS - Master process:'
|
|
write ( *, '(a,g14.6)' ) ' Estimate for PI is ', q_global
|
|
write ( *, '(a,g14.6)' ) ' Error is ', q_global - pi
|
|
|
|
wtime = MPI_Wtime ( ) - wtime
|
|
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a,f14.6)' ) ' Elapsed wall clock seconds = ', wtime
|
|
|
|
end if
|
|
!
|
|
! Terminate MPI.
|
|
!
|
|
call MPI_Finalize ( ierr )
|
|
!
|
|
! Terminate.
|
|
!
|
|
if ( id == 0 ) then
|
|
write ( *, '(a)' ) ' '
|
|
write ( *, '(a)' ) 'INTERVALS:'
|
|
write ( *, '(a)' ) ' Normal end of execution.'
|
|
write ( *, '(a)' ) ' '
|
|
call timestamp ( )
|
|
end if
|
|
|
|
stop
|
|
end
|
|
function f ( x )
|
|
|
|
!*****************************************************************************80
|
|
!
|
|
!! F is the function we are integrating.
|
|
!
|
|
! Licensing:
|
|
!
|
|
! This code is distributed under the GNU LGPL license.
|
|
!
|
|
! Modified:
|
|
!
|
|
! 10 February 2000
|
|
!
|
|
! Author:
|
|
!
|
|
! John Burkardt
|
|
!
|
|
! Parameters:
|
|
!
|
|
! Input, real ( kind = 8 ) X, the argument of the function.
|
|
!
|
|
! Output, real ( kind = 8 ) F, the value of the function.
|
|
!
|
|
implicit none
|
|
|
|
real ( kind = 8 ) f
|
|
real ( kind = 8 ) x
|
|
|
|
f = 4.0D+00 / ( 1.0D+00 + x * x )
|
|
|
|
return
|
|
end
|
|
subroutine timestamp ( )
|
|
|
|
!*****************************************************************************80
|
|
!
|
|
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
|
|
!
|
|
! Example:
|
|
!
|
|
! 31 May 2001 9:45:54.872 AM
|
|
!
|
|
! Licensing:
|
|
!
|
|
! This code is distributed under the GNU LGPL license.
|
|
!
|
|
! Modified:
|
|
!
|
|
! 18 May 2013
|
|
!
|
|
! Author:
|
|
!
|
|
! John Burkardt
|
|
!
|
|
! Parameters:
|
|
!
|
|
! None
|
|
!
|
|
implicit none
|
|
|
|
character ( len = 8 ) ampm
|
|
integer ( kind = 4 ) d
|
|
integer ( kind = 4 ) h
|
|
integer ( kind = 4 ) m
|
|
integer ( kind = 4 ) mm
|
|
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
|
|
'January ', 'February ', 'March ', 'April ', &
|
|
'May ', 'June ', 'July ', 'August ', &
|
|
'September', 'October ', 'November ', 'December ' /)
|
|
integer ( kind = 4 ) n
|
|
integer ( kind = 4 ) s
|
|
integer ( kind = 4 ) values(8)
|
|
integer ( kind = 4 ) y
|
|
|
|
call date_and_time ( values = values )
|
|
|
|
y = values(1)
|
|
m = values(2)
|
|
d = values(3)
|
|
h = values(5)
|
|
n = values(6)
|
|
s = values(7)
|
|
mm = values(8)
|
|
|
|
if ( h < 12 ) then
|
|
ampm = 'AM'
|
|
else if ( h == 12 ) then
|
|
if ( n == 0 .and. s == 0 ) then
|
|
ampm = 'Noon'
|
|
else
|
|
ampm = 'PM'
|
|
end if
|
|
else
|
|
h = h - 12
|
|
if ( h < 12 ) then
|
|
ampm = 'PM'
|
|
else if ( h == 12 ) then
|
|
if ( n == 0 .and. s == 0 ) then
|
|
ampm = 'Midnight'
|
|
else
|
|
ampm = 'AM'
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
|
|
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
|
|
|
|
return
|
|
end
|