Commits

Anonymous committed 9e85b76

Added mpi/quadrature example and f2py examples.

Comments (0)

Files changed (10)

+
+.PHONY: clean
+
+%.so : %.f90
+	f2py -m $* -c $<
+
+clean:
+	rm -f *.so *.pyc
+
+! $CLASSHG/codes/f2py/fcn1.f90
+
+! A test function for trying out f2py:
+!   $ f2py -m fcn1 -c fcn1.f90
+!   $ python
+!   >>> import fcn1
+!   >>> fcn1.f1(1.)
+!   2.7182818284590451
+
+
+function f1(x)
+    real(kind=8), intent(in) :: x
+    real(kind=8) :: f1
+    f1 = exp(x)
+end function f1

codes/f2py/jacobi1.f90

+
+! $CLASSHG/codes/f2py/jacobi1
+
+! For use with f2py:
+!   $ f2py -m jacobi1 -c jacobi1.f90
+! Then in Python you can import jacobi1 and use it to take a fixed number of
+! iterations using Jacobi in one space dimension to solve
+!      u(i-1) - 2*u(i) + u(i+1) = -dx**2 * f(i).
+! Python call:
+!   >>> u = jacobi1.iterate(u0, iters, f)
+! where 
+!   iters is the number of iterations to take,
+!   u0(0:n+1) is the initial guess (with desired boundary values 
+!             in u0(0) and u0(n+1), these are not changed.
+!   f(0:n+1) contains values of f(x) at each grid point.
+! 
+
+subroutine iterate(u0,iters,f,u,n)
+    implicit none
+    integer, intent(in) :: n, iters
+    real(kind=8), dimension(0:n+1), intent(in) :: u0, f
+    real(kind=8), dimension(0:n+1), intent(out) :: u
+
+    ! local variables:
+    real(kind=8), dimension(:), allocatable :: uold
+    real(kind=8) :: dx
+    integer :: i, iter
+
+    ! allocate storage for boundary points too:
+    allocate(uold(0:n+1))
+
+    ! grid spacing:
+    dx = 1.d0 / (n+1.d0)
+
+    uold = u0
+
+    ! initialize boundary values of u:
+    u(0) = u0(0)
+    u(n+1) = u0(n+1)
+
+    ! take the iterations requested:
+    do iter=1,iters
+        do i=1,n
+            u(i) = 0.5d0*(uold(i-1) + uold(i+1) + dx**2*f(i))
+            enddo
+        uold = u
+        enddo
+
+end subroutine iterate

codes/f2py/plot_jacobi_iterates.py

+
+from matplotlib import pyplot as plt
+import numpy as np
+import time
+
+import jacobi1   # created from jacobi1.f90 using f2py
+
+n = 51
+
+# grid points:
+x = np.linspace(0., 1., n)
+
+# boundary conditions:
+alpha = 20.
+beta = 60.
+
+# initial guess:
+u = np.linspace(alpha, beta, n)
+
+# right hand side:
+f = 100. * np.exp(x)
+
+niters = 0
+plt.clf()
+plt.plot(x, u, 'o-')
+plt.title("After %s iterations" % niters)
+plt.ylim([10., 80.])
+
+ans = raw_input("Type return to start... ")
+
+iters_per_plot = 200   # number of iterations between plots
+nplots = 20            # number of plots to produce
+
+for nn in range(nplots):
+    u = jacobi1.iterate(u, iters_per_plot, f)
+    plt.plot(x, u, 'o-')
+    niters = niters + iters_per_plot
+    plt.title("After %s iterations" % niters)
+    plt.ylim([10., 80.])
+    plt.draw()
+    time.sleep(.5)
+    
+
+
+
+! $CLASSHG/codes/f2py/sub1.f90
+
+! A test subroutine for trying out f2py:
+!   $ f2py -m sub1 -c sub1.f90
+!   $ python
+!   >>> import sub1
+!   >>> sub1.mysub(3., 5.)
+!   (8.0, -2.0)
+
+subroutine mysub(a,b,c,d)
+    real (kind=8), intent(in) :: a,b
+    real (kind=8), intent(out) :: c,d
+    c = a+b
+    d = a-b
+end subroutine mysub

codes/mpi/quadrature/Makefile

+# $UWHPSC/codes/fortran/mpi/quadrature/Makefile
+
+OBJECTS = functions_mod.o quadrature_mod.o testquad.o
+
+NUMPROCS ?= 4
+
+.PHONY: test plot clean
+
+FLAGS ?= 
+
+test: testquad.exe
+	mpiexec -n $(NUMPROCS) ./testquad.exe
+
+testquad.exe: $(OBJECTS)
+	mpif90 $(FLAGS) $(OBJECTS) -o testquad.exe
+
+%.o : %.f90
+	mpif90 -c $(FLAGS) $< 
+
+clean:
+	rm -f *.o *.exe *.mod
+

codes/mpi/quadrature/functions_mod.f90

+! $UWHPSC/codes/mpi/quadrature/functions_mod.f90
+
+module functions_mod
+
+    ! Define the function we want to integrate.
+    ! Also define the indefinite integral (anti-derivative) of g for
+    ! use in computing the true solution to compute the error.
+
+    ! The integer gevals can be initialized to 0 in the main program
+    ! and then will count how many times the function g is called by each proc.
+    ! Note that this is automatically local to each processor!
+
+    ! k is the wave number used in g and gint.
+
+    use mpi
+    implicit none
+    integer :: gevals
+    real(kind=8) :: k
+    save
+
+contains
+
+real(kind=8) function g(x)
+    ! Test function to integrate
+    implicit none
+    real(kind=8), intent(in) :: x
+    integer :: proc_num, ierr
+
+    g = exp(0.1d0*x) + cos(k*x)
+    gevals = gevals + 1
+
+end function g
+
+real(kind=8) function gint(x)
+    ! Anti-derivative of g
+    implicit none
+    real(kind=8), intent(in) :: x
+    gint = 10.d0*exp(0.1d0*x) + sin(k*x) / k
+end function gint
+
+end module functions_mod

codes/mpi/quadrature/quadrature_mod.f90

+! $UWHPSC/codes/mpi/quadrature/quadrature_mod.f90
+
+! Quadrature rules to approximate definite integrals.
+
+module quadrature_mod
+
+    use mpi
+
+contains
+
+
+subroutine simpson(f, a, b, n, integral)
+    implicit none
+    real(kind=8), intent(in) :: a,b
+    integer, intent(in) :: n
+    real(kind=8), intent(out) :: integral
+    real(kind=8), external :: f
+    ! local variables
+    real(kind=8) :: xi, dx
+    real(kind=8) :: int_proc, a_proc, b_proc, xi_proc
+    integer :: numprocs, ierr,points_per_proc,proc_num
+    integer :: i,istart,iend
+
+    call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
+    call MPI_COMM_RANK(MPI_COMM_WORLD, proc_num, ierr)
+
+    ! Determine how many points to handle with each proc
+    points_per_proc = (n + numprocs - 1)/numprocs
+    if (proc_num == 0) then   ! Only one proc should print to avoid clutter
+        print *, "points_per_proc = ", points_per_proc
+    end if
+
+    ! Determine start and end index for this proc's points
+    istart = proc_num * points_per_proc + 1
+    iend = min((proc_num + 1)*points_per_proc, n)
+
+    dx = (b-a) / n
+    integral = 0.d0
+
+    a_proc = a + dx*(istart-1)
+    b_proc = a + dx*iend
+    int_proc = f(a_proc) + f(b_proc)
+
+
+    ! Diagnostic: tell the user which points will be handled by which proc
+    ! Note: print everything in a single print statement so all the 
+    ! info for a given process gets printed together.
+
+    print 201, proc_num, istart, iend, a_proc, b_proc
+201 format("Process ",i2," will take i = ",i9," through i = ",i9, /, &
+           "   Corresponding to a_proc = ",e13.6," and b_proc = ",e13.6)
+
+    do i=istart+1,iend
+        xi_proc = a + dx*(i-1)
+        int_proc = int_proc + 2.d0*f(xi_proc)
+        enddo
+
+    do i=istart,iend
+        xi_proc = a + dx*(i-0.5d0)
+        int_proc = int_proc + 4.d0*f(xi_proc)
+        enddo
+
+    call MPI_REDUCE(int_proc,integral,1,MPI_DOUBLE_PRECISION,MPI_SUM,0, &
+                        MPI_COMM_WORLD,ierr)
+
+    if (proc_num == 0) then
+        integral = dx*integral / 6.d0
+    endif
+
+    ! broadcast to all other processes so they all return the right thing:
+    call MPI_BCAST(integral, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+
+end subroutine simpson
+
+
+end module quadrature_mod

codes/mpi/quadrature/testquad.f90

+! $UWHPSC/codes/mpi/quadrature/testquad.f90
+! Test the quadrature module.
+
+program main
+
+    use mpi
+    use quadrature_mod, only: simpson
+    use functions_mod, only: g, gint, gevals, k
+
+    implicit none
+    real(kind=8) :: a, b, int_approx, int_true, error
+    integer :: proc_num, ierr, n
+
+    call MPI_INIT(ierr)
+
+    k = 10.d0    ! wave number
+    n = 50000    ! number of intervals
+
+    a = 0.d0
+    b = 3.d0
+    int_true = gint(b) - gint(a)
+
+    gevals = 0  ! keeps track of number of g evals
+    call simpson(g, a, b, n, int_approx)
+    error = int_approx - int_true
+
+
+    call MPI_COMM_RANK(MPI_COMM_WORLD, proc_num, ierr)
+    if (proc_num==0) then
+        print 62, int_approx, int_true, error
+        endif
+
+    print 65, gevals, proc_num
+
+62  format("   approx = ",e20.13, /, &
+           "   true =   ",e20.13, /, &
+           "   error =  ",e10.3)
+65  format("   g was evaluated ",i9," times by proc ",i1)
+
+    call MPI_FINALIZE(ierr)
+
+end program main
     Do this as follows:
 
     * Modify the main program to call `MPI_INIT` and `MPI_FINALIZE`.
+      Note that with MPI, we must call `MPI_INIT` as the first statement in
+      the main program, so every process is running the same code, and 
+      every process will call the subroutine `many_walks`.  
+      See `$UWHPSC/codes/mpi/quadrature` for an example of how Simpson's
+      method might be implemented in MPI.
 
     * In the main program, use::