Commits

Randall LeVeque  committed dc842ed Merge

Merged in rjleveque/uwhpsc_rjl/coursera (pull request #10)

final project added

  • Participants
  • Parent commits 2690602, 471136e
  • Branches coursera

Comments (0)

Files changed (25)

File codes/mpi/matrix1norm2.f90

-! $UWHPSC/codes/mpi/matrix1norm1.f90
+! $UWHPSC/codes/mpi/matrix1norm2.f90
 !
 ! Compute 1-norm of a matrix using mpi.
 ! Process 0 is the master that sets things up and then sends a column

File notes/homeworks.rst

  * :ref:`homework4`: Week 6 
  * :ref:`homework5`: Week 8 
  * :ref:`homework6`: Week 9 
-
-The following will be available by Week 8:
-
- * Final Project: Week 10
+ * :ref:`project`: Week 10
 
 These are essentially the same assignments used in the UW course.
 

File notes/project.rst

+
+.. _project:
+
+==========================================
+Final Project
+==========================================
+
+
+
+The goals of this project are to:
+
+* Get some more experience with MPI and the master-worker paradigm.
+* Learn about Monte Carlo methods and random number generators.
+* Get experience computing in Fortran and plotting results in Python.
+
+See also:
+
+* :ref:`random`
+* :ref:`poisson`
+
+Create a new subdirectory `project` in your repository for the codes you
+create.
+
+**PART 1**
+
+Create a subdirectory `project/part1` for this part.
+
+    In the last part of :ref:`homework6`, you wrote an MPI code to compute
+    a Trapezoid approximation to an integral by splitting the original
+    interval up into `nsub = num_procs - 1` subintervals. The master Process
+    0 farmed these out to the remaining processes and collected the results.
+
+    Modify your code to make a new version `test3.f90`
+    that allows nsub to be set to a value larger (or smaller) than
+    `num_procs - 1`.  In this case the master process should send out
+    subintervals to the worker processes until all the work is done.  This
+    will require doing something similar to what is done in the sample code
+    `$UWHPSC/codes/mpi/matrix1norm2.f90` to keep track of which subintervals
+    have been handled.  
+
+    Have Process 0 read in the number of subintervals with the lines::
+
+            print *, "How many subintervals? "
+            read *, nsub
+
+    and then broadcast this value to the other processes.  
+
+    Test that your code works both when the number of subintervals is
+    smaller or larger than the number of processes specified with `mpiexec`.
+
+    If the code is run with a single process then it should halt with an
+    error message (since there are no workers to compute the integral over
+    subintervals in this case).
+
+    Provide a `Makefile` that compiles using `mpif90` and runs the code by 
+    default with 4 processes, as in Homework 6.
+
+    **Sample output**
+
+    Your code should produce output similar to this, for example... ::
+        
+        $ make test3
+        mpif90  -c  functions.f90 
+        mpif90  -c  quadrature.f90 
+        mpif90  -c  test3.f90 
+        mpif90  functions.o quadrature.o test3.o -o test3.exe
+        mpiexec -n 4    test3.exe
+        Using   4 processes
+         How many subintervals? 
+        2
+        true integral:   6.00136745954910E+00
+          
+        fevals by Process  0:             0
+        fevals by Process  1:          1000
+        fevals by Process  2:          1000
+        fevals by Process  3:             0
+        Trapezoid approximation with     2000 total points:   6.00125232481036E+00
+        Total number of fevals:       2000
+        
+        
+        
+        $ make test3
+        mpiexec -n 4    test3.exe
+        Using   4 processes
+         How many subintervals? 
+        7
+        true integral:   6.00136745954910E+00
+          
+        fevals by Process  0:             0
+        fevals by Process  1:          3000
+        fevals by Process  2:          2000
+        fevals by Process  3:          2000
+        Trapezoid approximation with     7000 total points:   6.00135820753458E+00
+        Total number of fevals:       7000
+        
+        
+        $ make test3 NUM_PROCS=1
+        mpiexec -n 1 test3.exe
+         *** Error: need to use at least two processes
+        
+
+**PART 2**
+
+Create a subdirectory `project/part2` for this part.
+
+    Monte Carlo methods are often used to estimate the values of definite
+    integrals in high dimensional spaces since traditional quadrature
+    methods based on regular grids may require too many points.  
+
+    The files in `$UWHPSC/codes/project/part2` contain much of what is
+    needed to experiment with a Monte Carlo approximation to the integral
+
+    :math:`\int_{a_1}^{b_1} \int_{a_2}^{b_2} \cdots \int_{a_d}^{b_d} g(x_1,x_2,\ldots,x_d) dx_1~dx_2~\cdots~dx_d`
+
+    over a rectangular region in :math:`d` space dimensions.  The Monte
+    Carlo approximation to the integral is given by 
+
+    :math:`\frac V N \sum_1^N g(x_1^{[k]},x_2^{[k]},\ldots,x_d^{[k]})`
+
+    where :math:`(x_1^{[k]},x_2^{[k]},\ldots,x_d^{[k]})` is the k'th
+    random point and :math:`V = (b_1-a_1)(b_2-a_2)\cdots(b_d-a_d)` is the
+    volume of the rectangular region of integration.
+
+    The main program in `test_quad_mc.f90` is set up to experiment with a
+    simple integral with varying number of Monte-Carlo points.  
+
+    What is missing is the module `quadrature_mc.f90`.  Create this module,
+    containing a function `quad_mc` with the calling sequence::
+
+        quad_mc(g, a, b, ndim, npoints)
+
+    that returns a Monte Carlo approximation to the integral, where:
+
+    * `g` is the function defining the integrand.  `g` takes two
+      arguments `x` and `ndim`, where `x` is an array of length `ndim`,
+      the number of dimensions we are integrating over.
+      (See the example in the `functions.f90` module.)
+
+    * `a` and `b` are arrays of length `ndim` that have the lower and upper
+      limits of integration in each dimension.
+
+    * `ndim` is the number of dimensions to integrate over.
+
+    * `npoints` is the number of Monte Carlo samples to use.
+
+
+    The random number generator should be called only once to generate all
+    the points needed and then the function `g` evaluated at appropriate
+    points.  Note that you will need `npoints*ndim` random numbers since
+    each point `x` has `ndim` components.
+
+    Allocate appropriate size arrays to manage this.
+
+    Note that the function :math:`g(x)` specified for this test is very
+    simple so that the true solution can be easily computed in any number of 
+    dimensions.
+
+    :math:`g(x) = x_1^2 + x_2^2 + \cdots + x_d^2`
+
+    The test program in `test_quad_mc.f90` computes the exact integral of
+    this over any rectangular region.  Convince yourself this is right.
+
+    Once you have provided a suitable module as described above,
+    running this code should give results like the following::
+
+        $ make plot
+        gfortran  -c  random_util.f90 
+        gfortran  -c  functions.f90 
+        gfortran  -c  quadrature_mc.f90 
+        gfortran  -c  test_quad_mc.f90 
+        gfortran  random_util.o functions.o quadrature_mc.o test_quad_mc.o -o
+        test.exe
+        ./test.exe
+        Testing Monte Carlo quadrature in 20 dimensions
+        True integral:   1.95734186666667E+08
+         seed1 for random number generator:       12345
+        Final approximation to integral:   1.95728471073896E+08
+         Total g evaluations:      1310720
+         python plot_mc_quad_error.py
+
+    A file `mc_quad_error.txt` should be created with the estimate of the
+    integral computed with varying number of random points and the error
+    in each.  
+
+    A plot of these results should also be created as `mc_quad_error.png`,
+    that looks like this:
+
+    .. image:: images/mc_quad_error.png
+       :width: 10cm
+
+
+    The test problem is set up to estimate a 20-dimensional integral.
+    Note that the relative error is plotted, which gives an indication
+    of the number of correct digits.  (Note that the absolute error is about 
+    2e8 times larger for this problem!)
+
+    *Note:* This problem should be quite easy; the code needed for
+    `quad_mc` should be short.  The main purpose of this problem is to
+    illustrate the basic structure of such a code, which you can follow
+    in the next problem.
+
+
+**PART 3**
+
+
+    The sample program `$UWHPSC/codes/project/part3/laplace_mc.py` 
+    can be run from
+    IPython to illustrate how a random walk on a lattice can be used to 
+    generate an approximate solution to the steady-state heat equation
+    at a single point.  This is described in more detail in the section
+    :ref:`poisson_mc`.
+    
+    Note that there is a parameter `plot_walk` that is set to `True` for
+    this demo.  If you set it to `False` and execute the code, then it will
+    take many more walks and print out the approximations as it repeatedly
+    doubles the number of walks taken.
+    
+    Using this as a model, write a Fortran code to approximate 
+    the solution to Laplace's equation at a single point :math:`(x_0,y_0)`
+    using the random walk approach.  
+
+    The module `$UWHPSC/codes/project/part3/problem_description.f90`
+    is a starting point. 
+
+    Supplement this with the following:
+
+    * A module `mc_walk.f90` containing two subroutines 
+
+      * `subroutine random_walk(i0, j0, max_steps, ub, iabort)`
+        based on the Python function `random_walk`.  
+        In the Fortran case `ub` should be an output variable with the
+        value of u at the boundary point reached, in the case when the walk
+        successfully reached the boundary.  In this case the subroutine
+        should return `iabort = 0`. If the walk did not reach the
+        boundary after `max_steps`, then `ub` can be anything, but 
+        return `iabort = 1` in this case.
+
+      * `subroutine many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)`
+        based on the Python equivalent.  In this case `u_mc` should be an
+        output variable with the average value of `u` on the boundary
+        computed based on the successful walks, and `n_success` is an output
+        variable telling how many were successful.  
+
+      * Add a module variable `nwalks` to this module that is initialized to
+        0 in the main program and incremented by one each time `random_walk`
+        is called.
+
+
+    * A main program named `laplace_mc.f90` that does something similar to
+      the main program in the Python code.  In particular it should:
+
+      * Set `x0, y0, i0, j0, max_steps` as in the Python.
+        You should `use` what's needed from the module
+        `problem_description.f90`.
+
+      * Initialize the random number generator.  You can use the 
+        `random_util.f90` module from Part 2 for this.
+        Set `seed1 = 12345`.
+
+      * Intialize `nwalks = 0` and print out at the end the value, which
+        should be the total number of times `random_walk` was called.
+
+      * Call `many_walks` first with `u_mc = 10` and then have a loop to
+        repeatedly double the number of samples and print out the
+        estimate of `u` and the relative error after each doubling.
+
+      * In addition, it should write the total number of walks, the estimate of
+        `u` and the relative error each doubling to a file named
+        `mc_laplace_error.txt` with the same format as the file
+        `mc_quad_error.txt` in Part 2.
+
+    * A python script `plot_mc_laplace_error.py` based on the plotting
+      script from Part 2 to produce a log-log plot of the results.
+
+    * A Makefile so that `make plot` will produce the `png` file.
+
+    The Fortran code does not need to include an option for plotting the
+    walks, that was just for demo purposes.
+
+    Note that the main program and each subroutine will have to `use`
+    various variables or subroutines from other modules.
+
+    **Sample output** ::
+
+        $ make plot
+        gfortran  -c  random_util.f90 
+        gfortran  -c  problem_description.f90 
+        gfortran  -c  mc_walk.f90 
+        gfortran  -c  laplace_mc.f90 
+        gfortran  random_util.o problem_description.o mc_walk.o laplace_mc.o -o
+        test.exe
+        ./test.exe
+         seed1 for random number generator:       12345
+                10  0.377000000000000E+00   0.162222E+00
+                20  0.408125000000000E+00   0.930556E-01
+                40  0.452875000000000E+00   0.638889E-02
+                80  0.436125000000000E+00   0.308333E-01
+               160  0.440656250000000E+00   0.207639E-01
+               320  0.468687500000000E+00   0.415278E-01
+               640  0.460773437500000E+00   0.239410E-01
+              1280  0.455091796874999E+00   0.113151E-01
+              2560  0.455277343749997E+00   0.117274E-01
+              5120  0.455505371093748E+00   0.122342E-01
+             10240  0.456198974609378E+00   0.137755E-01
+             20480  0.454078369140635E+00   0.906304E-02
+             40960  0.450970458984394E+00   0.215658E-02
+        Final approximation to u(x0,y0):   4.50970458984394E-01
+        Total number of random walks:      40960
+        python plot_mc_laplace_error.py
+
+    Note that with `max_steps = 100*max(nx,ny)` all of the walks
+    successfully reached the boundary.  You might try with a smaller
+    value such as `max_steps = 10` in which case many walks will fail.
+    In this case you might see results like the following
+    (*Note that the original results shown here were incorrect!*) ::
+
+         seed1 for random number generator:       12345
+                 4  0.697500000000000E+00   0.550000E+00
+                 8  0.632500000000000E+00   0.405556E+00
+                17  0.608529411764706E+00   0.352288E+00
+                31  0.623548387096774E+00   0.385663E+00
+                71  0.622042253521127E+00   0.382316E+00
+               134  0.616623134328358E+00   0.370274E+00
+               268  0.623619402985075E+00   0.385821E+00
+               553  0.620099457504520E+00   0.377999E+00
+              1092  0.623298992673990E+00   0.385109E+00
+              2184  0.622995650183145E+00   0.384435E+00
+              4416  0.624125339673914E+00   0.386945E+00
+              8765  0.625060182544217E+00   0.389023E+00
+             17623  0.624690319468906E+00   0.388201E+00
+        Final approximation to u(x0,y0):   6.24690319468906E-01
+        Total number of random walks:      40960
+
+
+    The total number of walks `nwalks` is the same, but fewer were used
+    in each estimate of the solution.  
+
+    Note that the Fortran equivalent of `floor` is `nint`, e.g.
+    `i0 = nint((x0-ax)/dx)` determines the index of the gridpoint in the
+    `x` direction closest to the point `x0`.  
+
+**PART 4**
+
+
+    Parallelize the code from Part 3 using MPI.  
+    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::
+
+        seed1 = 12345   
+        seed1 = seed1 + 97*proc_num  ! unique for each process
+        call init_random_seed(seed1)
+
+      so that each process will generate a unique set of random numbers.
+
+    * Modify subroutine `many_walks` so that Process 0 is the master
+      whose job is to farm out all of the `n_mc` walks requested
+      to each of the other processes.  Follow the master-worker paradigm for
+      this.  This is a sensible way to try to do load balancing since some
+      walks will take many more steps than others.  (It would be better to
+      ask each worker to do some number of walks greater than 1 each time so
+      that there is less communication, but let's keep it simple.)
+
+      Note that the master does not have to send any data to a worker,
+      just an empty message requesting another walk, so it could send 
+      `MPI_BOTTOM` and use `tag = 1` to indicate this is a request for
+      another walk.  Use `tag = 0` to indicate to a worker that it is done.
+
+      The worker will have to receive from the master with `MPI_ANY_TAG` and
+      then check `status(MPI_TAG)` to see what it needs to do.
+
+      If another walk is requested, the worker should call `random_walk` and
+      then send back to the Master the result as a single data value of type
+      `MPI_DOUBLE_PRECISION`.   For this message set the `tag` to the value
+      of `iabort` that was returned from the call to `random_walk` so that
+      the Master knows whether to include this walk in the accumulated 
+      Monte Carlo result.
+
+    * Recall that with MPI every process is executing the same code but that
+      all data is local to a process.   So the basic structure of the main
+      program can remain the same.  Every process will execute the loop that
+      repeatedly increases the size of `n_mc` and every process will call
+      `many_walks`.  But only the master process will return values of 
+      `u_mc` and `n_success` that are sensible, and so this process should 
+      update `u_mc_total` and print out the values to the screen and the file
+      `mc_laplace_error.txt`.
+
+    * The module variable
+      `nwalks` that is incremented in `random_walk` will be local to each
+      process. In the main program, at the end have each process print out how
+      many walks it took and use `MPI_REDUCE` to compute the total number of
+      walks taken by all processes and have Process 0 print this value.
+
+    * Create a `Makefile` that works for this by combining aspects of those
+      from Part 1 (for MPI) and Part 3 (for the targets needed).
+      
+    **Sample output** ::
+
+        $ make plot
+        mpif90  -c  random_util.f90 
+        mpif90  -c  problem_description.f90 
+        mpif90  -c  mc_walk.f90 
+        mpif90  -c  laplace_mc.f90 
+        mpif90   random_util.o problem_description.o mc_walk.o laplace_mc.o -o
+        test.exe
+        mpiexec -n 4    test.exe
+         seed1 for random number generator:       12442
+         seed1 for random number generator:       12539
+         seed1 for random number generator:       12636
+         seed1 for random number generator:       12345
+                10  0.516750000000000E+00   0.148333E+00
+                20  0.478500000000000E+00   0.633333E-01
+                40  0.425437500000000E+00   0.545833E-01
+                80  0.431562500000000E+00   0.409722E-01
+               160  0.431593750000000E+00   0.409028E-01
+               320  0.425703125000000E+00   0.539931E-01
+               640  0.426492187500000E+00   0.522396E-01
+              1280  0.427759765624999E+00   0.494227E-01
+              2560  0.430487304687498E+00   0.433615E-01
+              5120  0.443433105468749E+00   0.145931E-01
+             10240  0.449190429687505E+00   0.179905E-02
+             20480  0.449556518554698E+00   0.985514E-03
+             40960  0.451413696289083E+00   0.314155E-02
+        Final approximation to u(x0,y0):   4.51413696289083E-01
+        Total walks performed by all processes:      40960
+        Walks performed by Process  0:          0
+        Walks performed by Process  1:      12928
+        Walks performed by Process  2:      13414
+        Walks performed by Process  3:      14618
+        python plot_mc_laplace_error.py
+
+
+
+When finished
+-------------
+
+Your project directory should contain:
+
+* part1/functions.f90
+* part1/quadrature.f90
+* part1/test3.f90
+* part1/Makefile
+
+* part2/functions.f90
+* part2/quadrature_mc.f90
+* part2/random_util.f90
+* part2/test_quad_mc.f90
+* part2/plot_mc_quad_error.py
+* part2/Makefile
+
+* part3/problem_description.f90
+* part3/laplace_mc.f90
+* part3/mc_walk.f90
+* part3/random_util.f90
+* part3/plot_mc_laplace_error.py
+* part3/Makefile
+
+* part4/problem_description.f90
+* part4/laplace_mc.f90
+* part4/mc_walk.f90
+* part4/random_util.f90
+* part4/plot_mc_laplace_error.py
+* part4/Makefile
+
+Congratulations, you have completed the assignments for this course!

File solutions/project/part1/Makefile

+
+OBJECTS = functions.o quadrature.o test3.o
+FFLAGS = 
+NUM_PROCS ?= 4   # default if not specified on command line or env variable
+
+.PHONY: test3 clean 
+
+test3: test3.exe
+	mpiexec -n $(NUM_PROCS) test3.exe
+
+test3.exe: $(OBJECTS)
+	mpif90 $(FFLAGS) $(OBJECTS) -o test3.exe
+
+%.o : %.f90
+	mpif90 $(FFLAGS) -c  $< 
+
+clean:
+	rm -f *.o *.exe *.mod
+

File solutions/project/part1/functions.f90

+
+module functions
+
+    implicit none
+    integer :: fevals_proc
+    real(kind=8) :: k
+    save
+
+contains
+
+    real(kind=8) function f(x)
+        implicit none
+        real(kind=8), intent(in) :: x 
+        integer proc_num, ierr
+
+        ! keep track of number of function evaluations by
+        ! each process:
+        fevals_proc = fevals_proc + 1
+        
+        f = 1.d0 + x**3 + sin(k*x)
+        
+    end function f
+
+end module functions

File solutions/project/part1/quadrature.f90

+
+module quadrature
+
+contains
+
+real(kind=8) function trapezoid(f, a, b, n)
+
+    ! Estimate the integral of f(x) from a to b using the
+    ! Trapezoid Rule with n points.
+
+    ! Input:
+    !   f:  the function to integrate
+    !   a:  left endpoint
+    !   b:  right endpoint
+    !   n:  number of points to use
+    ! Returns:
+    !   the estimate of the integral
+     
+    implicit none
+    real(kind=8), intent(in) :: a,b
+    real(kind=8), external :: f
+    integer, intent(in) :: n
+
+    ! Local variables:
+    integer :: j
+    real(kind=8) :: h, trap_sum, xj
+
+    h = (b-a)/(n-1)
+    trap_sum = 0.5d0*(f(a) + f(b))  ! endpoint contributions
+    
+    do j=2,n-1
+        xj = a + (j-1)*h
+        trap_sum = trap_sum + f(xj)
+        enddo
+
+    trapezoid = h * trap_sum
+
+end function trapezoid
+
+end module quadrature
+

File solutions/project/part1/test3.f90

+
+program test
+
+    use mpi
+
+    use quadrature, only: trapezoid
+    use functions, only: f, fevals_proc, k
+
+    implicit none
+    real(kind=8) :: a,b,int_true, int_approx, int_sub, dx_sub
+    real(kind=8) :: ab_sub(2)   ! to hold a and b for each subinterval
+
+    integer :: i,j,jj,n_proc,proc_num, num_procs,ierr,n,fevals_total,nsub
+    integer :: nextsub, numsent, sender
+    integer, dimension(MPI_STATUS_SIZE) :: status
+
+    call MPI_INIT(ierr)
+    call MPI_COMM_SIZE(MPI_COMM_WORLD, num_procs, ierr)
+    call MPI_COMM_RANK(MPI_COMM_WORLD, proc_num, ierr)
+
+    ! All processes set these values so we don't have to broadcast:
+    k = 1.d3   ! functions module variable 
+    a = 0.d0
+    b = 2.d0
+    int_true = (b-a) + (b**4 - a**4) / 4.d0 - (1.d0/k) * (cos(k*b) - cos(k*a))
+    n = 1000
+    !nsub = 4
+
+    if (num_procs == 1) then
+        print *, "*** Error: need to use at least two processes"
+        go to 999
+        endif
+
+
+    ! Each process keeps track of number of fevals:
+    fevals_proc = 0
+
+    if (proc_num==0) then
+        int_approx = 0.d0
+        print 100, num_procs
+100     format("Using ",i3," processes")
+
+        print *, "How many subintervals? "
+        read *, nsub
+
+        print 10, int_true
+ 10     format("true integral: ", es22.14)
+        print *, " "  ! blank line
+        endif
+
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! wait for process 0 to print
+    call MPI_BCAST(nsub, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+
+    ! -----------------------------------------
+    ! code for Master (Processor 0):
+    ! -----------------------------------------
+
+    if (proc_num == 0) then
+
+      numsent = 0 ! keep track of how many subintervals sent
+
+      dx_sub = (b-a) / nsub
+
+      ! send the first batch to get all workers working:
+      do j=1,min(num_procs-1, nsub)
+        ab_sub(1) = a + (j-1)*dx_sub
+        ab_sub(2) = a + j*dx_sub
+        call MPI_SEND(ab_sub, 2, MPI_DOUBLE_PRECISION, j, j, &
+                      MPI_COMM_WORLD, ierr)
+        numsent = numsent + 1
+        enddo
+
+      ! as results come back, send out more work...
+      ! the variable sender tells who sent back a result and ready for more
+      do j=1,nsub
+        call MPI_RECV(int_sub, 1, MPI_DOUBLE_PRECISION, &
+                        MPI_ANY_SOURCE, MPI_ANY_TAG, &
+                        MPI_COMM_WORLD, status, ierr)
+        sender = status(MPI_SOURCE)
+        jj = status(MPI_TAG)
+        ! print *,"+++ int_sub, int_approx: ",int_sub, int_approx
+        int_approx = int_approx + int_sub
+
+        if (numsent < nsub) then
+            ! still more work to do, the next subinterval will be sent and
+            ! this index also used as the tag:
+            nextsub = numsent + 1 
+            ab_sub(1) = a + (nextsub-1)*dx_sub
+            ab_sub(2) = a + nextsub*dx_sub
+            
+            call MPI_SEND(ab_sub, 2, MPI_DOUBLE_PRECISION, sender, nextsub, &
+                          MPI_COMM_WORLD, ierr)
+            numsent = numsent + 1
+          else
+            ! send an empty message with tag=0 to indicate this worker
+            ! is done:
+            call MPI_SEND(MPI_BOTTOM, 0, MPI_INTEGER, &
+                            sender, 0, MPI_COMM_WORLD, ierr)
+          endif
+         enddo
+        endif
+
+    ! -----------------------------------------
+    ! code for Workers (Processors 1, 2, ...):
+    ! -----------------------------------------
+    if (proc_num /= 0) then
+
+        if (proc_num > nsub) go to 99   ! no work expected
+
+        do while (.true.)
+            ! repeat until message with tag==0 received...
+
+            call MPI_RECV(ab_sub, 2, MPI_DOUBLE_PRECISION, &
+                          0, MPI_ANY_TAG, &
+                          MPI_COMM_WORLD, status, ierr)
+
+            j = status(MPI_TAG)   ! this is the subinterval number
+
+            !print 12, proc_num, j       ! for debugging
+12          format("+++ Process ",i4,"  received message with tag ",i6)
+
+            if (j==0) go to 99    ! received "done" message
+
+            int_sub = trapezoid(f,ab_sub(1),ab_sub(2),n)
+
+            call MPI_SEND(int_sub, 1, MPI_DOUBLE_PRECISION, &
+                        0, j, MPI_COMM_WORLD, ierr)
+
+            enddo
+        endif
+
+99  continue   ! might jump to here if finished early
+
+
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! wait for all process to print
+
+    ! print the number of function evaluations by each thread:
+    print 101,  proc_num, fevals_proc
+101 format("fevals by Process ",i2,": ",i13)
+
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! wait for all process to print
+
+    call MPI_REDUCE(fevals_proc, fevals_total, 1, MPI_INTEGER, MPI_SUM, 0, &
+                    MPI_COMM_WORLD, ierr)
+
+    if (proc_num==0) then
+        print 13, nsub*n, int_approx
+13      format("Trapezoid approximation with ",i8," total points: ",es22.14)
+        print 102, fevals_total
+102     format("Total number of fevals: ",i10)
+        endif
+
+999 continue
+    call MPI_FINALIZE(ierr)
+
+end program test

File solutions/project/part2/Makefile

+
+OBJECTS = random_util.o functions.o quadrature_mc.o test_quad_mc.o
+FFLAGS = 
+
+.PHONY: test plot clean 
+
+test: test.exe
+	./test.exe
+
+test.exe: $(OBJECTS)
+	gfortran $(FFLAGS) $(OBJECTS) -o test.exe
+
+%.o : %.f90
+	gfortran $(FFLAGS) -c  $< 
+
+mc_quad_error.txt: test.exe
+	./test.exe
+
+plot: mc_quad_error.txt
+	python plot_mc_quad_error.py
+
+clean:
+	rm -f *.o *.exe *.mod
+
+clobber: clean
+	rm -f mc_quad_error.txt mc_quad_error.png

File solutions/project/part2/functions.f90

+
+module functions
+
+    implicit none
+    integer :: gevals
+    save
+
+contains
+
+    function g(x,ndim)
+
+    implicit none
+    integer, intent(in) :: ndim
+    real(kind=8), intent(in) :: x(ndim)
+    real(kind=8) :: g
+    integer :: i
+
+    g = 0.d0
+    do i=1,ndim
+        g = g + x(i)**2
+        enddo
+
+    gevals = gevals + 1
+
+    end function g
+
+end module functions

File solutions/project/part2/plot_mc_quad_error.py

+
+from pylab import *
+
+# read in three columns from file and unpack into 3 arrays:
+n,int_approx,error = loadtxt('mc_quad_error.txt',unpack=True)
+
+figure(1)
+clf()
+loglog(n,error,'-o',label='Monte-Carlo')
+loglog([1,1e7],[1,sqrt(1e-7)],'k',label='1 / sqrt(N)')
+legend()
+xlabel('number of MC points used')
+ylabel('abs(error)')
+title('Log-log plot of relative error in MC quadrature')
+savefig('mc_quad_error.png')

File solutions/project/part2/quadrature_mc.f90

+
+module quadrature_mc
+
+contains
+
+function quad_mc(g, a, b, ndim, npoints)
+
+    implicit none
+    real(kind=8),external :: g
+    integer, intent(in) :: ndim, npoints
+    real(kind=8), intent(in) :: a(ndim), b(ndim)
+    real(kind=8) :: quad_mc
+
+    integer :: i
+
+    real(kind=4), allocatable :: r(:), rr(:)
+    real(kind=8) :: gsum, volume
+    real(kind=8), allocatable :: x(:)
+
+    allocate(r(ndim*npoints))
+    allocate(x(ndim), rr(ndim))
+
+    gsum = 0.d0
+
+    volume = product(b-a)
+
+    call random_number(r)
+    do i=1,npoints
+        rr = r((i-1)*ndim + 1 : i*ndim)  ! next ndim elements of r
+        x = a + rr*(b-a)   ! sets ndim elements of x elementwise
+        gsum = gsum + g(x,ndim)
+        enddo
+    quad_mc = volume * gsum /  npoints
+
+    deallocate(x,r,rr)
+
+end function quad_mc
+
+end module quadrature_mc

File solutions/project/part2/random_util.f90

+
+module random_util
+
+contains
+
+subroutine init_random_seed(seed1)
+
+    ! Seed the random number generator.
+    ! If seed1 = 0 then set seed1 randomly using the system clock.
+    !              This will give different sequences of random numbers
+    !              from different runs, so results are not reproducible.
+    ! If seed1 > 0 then results will be reproducible since calling this
+    !              twice with the same seed will initialize the random
+    !              number generator so the same sequence is generated.
+    ! Once seed1 is set, set the other elements of the seed array by adding
+    ! multiples of 37 as suggested in the documentation.
+    ! The length of the seed array is determined by calling  
+    !    
+    integer, intent(inout) :: seed1
+
+    integer :: nseed, clock 
+    integer, dimension(:), allocatable :: seed
+
+    call random_seed(size = nseed)  ! determine how many numbers needed to seed
+    allocate(seed(nseed)) 
+
+    if (seed1 == 0) then
+        ! randomize the seed: not repeatable
+        call system_clock(count = clock)
+        seed1 = clock
+      endif
+
+    do i=1,nseed
+        seed(i) = seed1 + 37*(i-1)  
+        enddo
+
+    print *, "seed1 for random number generator:", seed1
+
+    call random_seed(put = seed)   ! seed the generator
+    deallocate(seed)
+
+end subroutine init_random_seed
+
+end module random_util

File solutions/project/part2/test_quad_mc.f90

+
+program test_quad_mc
+
+	use functions, only: g, gevals
+    use quadrature_mc, only: quad_mc
+    use random_util, only: init_random_seed
+
+    implicit none
+    integer, parameter :: ndim = 20
+    real(kind=8), dimension(ndim) :: a, b
+    real(kind=8) :: volume, int_true, int_mc, int_mc_total, error
+    integer :: i, seed1, n_total
+    integer :: num_chunks, points_per_chunk, npoints
+
+    gevals = 0
+    open(unit=25, file='mc_quad_error.txt', status='unknown')
+
+    do i=1,ndim
+        a(i) = 2.d0
+        b(i) = 4.d0
+        enddo
+
+    ! compute the true integral for special case where
+    ! g(x) = sum(x**2) over all dimensions  -- think about why this works!
+    volume = product(b-a)  ! =  product of b(i)-a(i) of ndim dimensions
+    int_true = volume * sum((b**3 - a**3) / (3.d0*(b-a)))
+
+    print '("Testing Monte Carlo quadrature in ",i2," dimensions")', ndim
+    print '("True integral: ", es22.14)', int_true
+
+
+    seed1 = 12345   ! or set to 0 for random seed
+    call init_random_seed(seed1)
+
+    ! Start with Monte Carlo using only a few points:
+    npoints = 10
+    int_mc = quad_mc(g,a,b,ndim,npoints)
+
+    ! start accumulating totals:
+    int_mc_total = int_mc
+    n_total = npoints
+
+    ! relative error gives estimate of number of correct digits:
+    error = abs(int_mc_total - int_true) / abs(int_true)
+    write(25,'(i10,e23.15,e15.6)') n_total, int_mc_total, error
+
+    ! loop to successively double the number of points used:
+
+    do i=1,17
+        ! compute approximation with new chunk of points
+        int_mc = quad_mc(g,a,b,ndim,npoints)
+
+        ! update estimate of integral based on doubling number:
+        int_mc_total = 0.5*(int_mc_total + int_mc)
+        n_total = n_total + npoints
+
+        ! relative error:
+        error = abs(int_mc_total - int_true) / abs(int_true)
+
+        write(25,'(i10,e23.15,e15.6)') n_total, int_mc_total, error
+        print '("With ",i9," points, the error is ",es15.8)', n_total,error
+
+        npoints = 2*npoints
+        enddo
+
+    print '("Final approximation to integral: ",es22.14)',int_mc_total
+    print *, "Total g evaluations: ",gevals
+
+end program test_quad_mc

File solutions/project/part3/Makefile

+
+OBJECTS = random_util.o problem_description.o mc_walk.o laplace_mc.o
+FFLAGS = 
+
+.PHONY: test plot clean clobber
+
+test: test.exe
+	./test.exe
+
+test.exe: $(OBJECTS)
+	gfortran $(FFLAGS) $(OBJECTS) -o test.exe
+
+%.o : %.f90
+	gfortran $(FFLAGS) -c  $< 
+
+
+mc_laplace_error.txt: test.exe
+	./test.exe
+
+plot: mc_laplace_error.txt
+	python plot_mc_laplace_error.py
+
+clean:
+	rm -f *.o *.exe *.mod
+
+clobber: clean
+	rm -f mc_laplace_error.txt mc_laplace_error.png

File solutions/project/part3/laplace_mc.f90

+
+program laplace_mc
+
+    use problem_description, only: utrue,nx,ny,ax,ay,dx,dy
+    use random_util, only: init_random_seed
+    use mc_walk, only: random_walk, many_walks, nwalks
+
+    implicit none
+    integer :: i,i0,j0,max_steps,seed1,n_mc,n_success,n_total
+    real(kind=8) :: x0,y0,u_true,u_mc,u_sum_old,u_sum_new, u_mc_total,error
+
+    open(unit=25, file='mc_laplace_error.txt', status='unknown')
+
+    x0 = 0.9d0
+    y0 = 0.6d0
+    i0 = nint((x0-ax)/dx)
+    j0 = nint((y0-ay)/dy)
+
+    ! shift (x0,y0) to a grid point if it wasn't already:
+    x0 = ax + i0*dx
+    y0 = ay + j0*dy
+
+    max_steps = 100*max(nx,ny)
+    !max_steps = 10
+
+    nwalks = 0   ! to keep track of calls to random_walk
+
+    u_true = utrue(x0, y0)
+
+    seed1 = 12345   ! or set to 0 for random seed
+    call init_random_seed(seed1)
+
+    n_mc = 10
+    call many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+
+    ! start accumulating totals:
+    u_mc_total = u_mc
+    n_total = n_success
+
+    ! relative error gives estimate of number of correct digits:
+    error = abs(u_mc_total - u_true) / abs(u_true)
+    write(25,'(i10,e23.15,e15.6)') n_total, u_mc_total, error
+    print '(i10,e23.15,e15.6)', n_total, u_mc_total, error
+
+	
+    ! loop to add successively double the number of points used:
+
+    do i=1,12
+        ! compute approximation with new chunk of points
+        call many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+
+        ! update estimate of u based on doubling number:
+        ! but not that not all were successful:
+        u_sum_old = u_mc_total * n_total
+        u_sum_new = u_mc * n_success
+        n_total = n_total + n_success
+        u_mc_total = (u_sum_old + u_sum_new) / n_total
+
+        ! relative error:
+        error = abs(u_mc_total - u_true) / abs(u_true)
+
+        write(25,'(i10,e23.15,e15.6)') n_total, u_mc_total, error
+        print '(i10,e23.15,e15.6)', n_total, u_mc_total, error
+
+        n_mc = 2*n_mc
+        enddo
+
+    print '("Final approximation to u(x0,y0): ",es22.14)',u_mc_total
+    print '("Total number of random walks: ",i10)',nwalks
+    
+
+end program laplace_mc

File solutions/project/part3/mc_walk.f90

+
+module mc_walk
+
+implicit none
+integer :: nwalks
+
+contains
+
+
+subroutine random_walk(i0, j0, max_steps, ub, iabort)
+
+    use problem_description, only: nx,ny,ax,ay,dx,dy,uboundary
+
+    implicit none
+    integer, intent(in) :: i0,j0,max_steps
+    real(kind=8), intent(out) :: ub
+    integer, intent(out) :: iabort
+
+    real(kind=4), allocatable :: r(:)
+    integer i,j,istep
+    real(kind=4) :: rr
+    real(kind=8) :: xb,yb
+
+    allocate(r(max_steps))
+
+    i = i0
+    j = j0
+    call random_number(r)
+
+    do istep=1,max_steps
+        rr = r(istep)
+        if (rr < 0.25d0) then
+            i = i-1
+          else if (rr < 0.5d0) then
+            i = i+1
+          else if (rr < 0.75d0) then
+            j = j-1
+          else 
+            j = j+1
+          endif
+
+        if (i*j*(nx+1-i)*(ny+1-j) == 0) then
+            xb = ax + i*dx
+            yb = ay + j*dy
+            ub = uboundary(xb, yb)
+            exit  ! end the walk
+            endif
+        enddo
+
+    
+    if (istep == max_steps+1) then
+        iabort = 1
+        ub = 0.d0  ! this value should not be used
+    else
+        iabort = 0
+    endif
+
+    nwalks = nwalks + 1
+        
+end subroutine random_walk
+
+
+subroutine many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+    implicit none
+    integer, intent(in) :: i0,j0,max_steps,n_mc
+    integer, intent(out) :: n_success
+    real(kind=8), intent(out) :: u_mc
+
+    integer :: k, iabort, nb_sum
+    real(kind=8) :: ub,ub_sum
+
+    ub_sum = 0.d0
+    nb_sum = 0
+
+    do k=1,n_mc
+        call random_walk(i0, j0, max_steps, ub, iabort)
+        if (iabort == 0) then
+            ub_sum = ub_sum + ub
+            nb_sum = nb_sum + 1
+            endif
+        enddo
+
+    u_mc = ub_sum / nb_sum
+    n_success = nb_sum
+    
+end subroutine many_walks
+
+
+end module mc_walk

File solutions/project/part3/plot_mc_laplace_error.py

+
+from pylab import *
+
+# read in three columns from file and unpack into 3 arrays:
+n,int_approx,error = loadtxt('mc_laplace_error.txt',unpack=True)
+
+figure(1)
+clf()
+loglog(n,error,'-o',label='Monte-Carlo')
+loglog([1,1e7],[1,sqrt(1e-7)],'k',label='1 / sqrt(N)')
+legend()
+xlabel('number of random walks taken')
+ylabel('abs(error)')
+title('Log-log plot of relative error in MC Laplace')
+savefig('mc_laplace_error.png')

File solutions/project/part3/problem_description.f90

+
+module problem_description
+
+    implicit none
+    real(kind=8), parameter :: ax = 0.0d0
+    real(kind=8), parameter :: bx = 1.0d0
+    real(kind=8), parameter :: ay = 0.4d0
+    real(kind=8), parameter :: by = 1.0d0
+    integer, parameter :: nx = 19
+    integer, parameter :: ny = 11
+    real(kind=8), parameter :: dx = (bx - ax) / (nx+1)
+    real(kind=8), parameter :: dy = (by - ay) / (ny+1)
+
+contains
+
+function utrue(x, y)
+
+    ! True solution for comparison, if known.
+
+    implicit none
+    real(kind=8), intent(in) :: x,y
+    real(kind=8) :: utrue
+
+    utrue = x**2 - y**2
+
+end function utrue
+
+function uboundary(x, y)
+
+    ! Return u(x,y) assuming (x,y) is a boundary point.
+
+    implicit none
+    real(kind=8), intent(in) :: x,y
+    real(kind=8) :: uboundary
+
+    if ((x-ax)*(x-bx)*(y-ay)*(y-by) .ne. 0.d0) then
+        print *, "*** Error -- called uboundary at non-boundary point"
+        stop
+        endif
+
+    uboundary = utrue(x,y)   ! assuming we know this
+
+end function uboundary
+    
+
+end module problem_description

File solutions/project/part3/random_util.f90

+
+module random_util
+
+contains
+
+subroutine init_random_seed(seed1)
+
+    ! Seed the random number generator.
+    ! If seed1 = 0 then set seed1 randomly using the system clock.
+    !              This will give different sequences of random numbers
+    !              from different runs, so results are not reproducible.
+    ! If seed1 > 0 then results will be reproducible since calling this
+    !              twice with the same seed will initialize the random
+    !              number generator so the same sequence is generated.
+    ! Once seed1 is set, set the other elements of the seed array by adding
+    ! multiples of 37 as suggested in the documentation.
+    ! The length of the seed array is determined by calling  
+    !    
+    integer, intent(inout) :: seed1
+
+    integer :: nseed, clock 
+    integer, dimension(:), allocatable :: seed
+
+    call random_seed(size = nseed)  ! determine how many numbers needed to seed
+    allocate(seed(nseed)) 
+
+    if (seed1 == 0) then
+        ! randomize the seed: not repeatable
+        call system_clock(count = clock)
+        seed1 = clock
+      endif
+
+    do i=1,nseed
+        seed(i) = seed1 + 37*(i-1)  
+        enddo
+
+    print *, "seed1 for random number generator:", seed1
+
+    call random_seed(put = seed)   ! seed the generator
+    deallocate(seed)
+
+end subroutine init_random_seed
+
+end module random_util

File solutions/project/part4/Makefile

+
+OBJECTS = random_util.o problem_description.o mc_walk.o laplace_mc.o
+FFLAGS = 
+NUM_PROCS ?= 4   # default if not specified on command line or env variable
+
+
+.PHONY: test plot clean clobber
+
+test: test.exe
+	mpiexec -n $(NUM_PROCS) test.exe
+
+test.exe: $(OBJECTS)
+	mpif90  $(FFLAGS) $(OBJECTS) -o test.exe
+
+%.o : %.f90
+	mpif90 $(FFLAGS) -c  $< 
+
+
+mc_laplace_error.txt: test.exe
+	mpiexec -n $(NUM_PROCS) test.exe
+
+plot: mc_laplace_error.txt
+	python plot_mc_laplace_error.py
+
+clean:
+	rm -f *.o *.exe *.mod
+
+clobber: clean
+	rm -f mc_laplace_error.txt mc_laplace_error.png

File solutions/project/part4/laplace_mc.f90

+
+program laplace_mc
+
+    use mpi
+
+    use problem_description, only: utrue,nx,ny,ax,ay,dx,dy
+    use random_util, only: init_random_seed
+    use mc_walk, only: random_walk, many_walks, nwalks
+
+    implicit none
+    integer :: i,i0,j0,max_steps,seed1,n_mc,n_success,n_total,nwalks_total
+    real(kind=8) :: x0,y0,u_true,u_mc,u_sum_old,u_sum_new, u_mc_total,error
+
+    integer :: num_procs, proc_num, ierr
+
+    call MPI_INIT(ierr)
+    call MPI_COMM_SIZE(MPI_COMM_WORLD, num_procs, ierr)
+    call MPI_COMM_RANK(MPI_COMM_WORLD, proc_num, ierr)
+
+
+    open(unit=25, file='mc_laplace_error.txt', status='unknown')
+
+    x0 = 0.9d0
+    y0 = 0.6d0
+    i0 = nint((x0-ax)/dx)
+    j0 = nint((y0-ay)/dy)
+
+    ! shift (x0,y0) to a grid point if it wasn't already:
+    x0 = ax + i0*dx
+    y0 = ay + j0*dy
+
+    max_steps = 100*max(nx,ny)
+
+    u_true = utrue(x0, y0)
+
+    seed1 = 12345   ! or set to 0 for random seed
+    seed1 = seed1 + 97*proc_num  ! unique for each process
+    call init_random_seed(seed1)
+
+    nwalks = 0  ! to keep track of how many walks each process does
+
+    n_mc = 10
+    call many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+
+    if (proc_num == 0) then
+
+        ! start accumulating totals:
+        u_mc_total = u_mc
+        n_total = n_success
+
+        ! relative error gives estimate of number of correct digits:
+        error = abs(u_mc_total - u_true) / abs(u_true)
+        print '(i10,e23.15,e15.6)', n_total, u_mc_total, error
+        write(25,'(i10,e23.15,e15.6)') n_total, u_mc_total, error
+        endif
+
+    ! loop to add successively double the number of points used:
+
+    do i=1,12
+        ! compute approximation with new chunk of points
+        call many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+
+        if (proc_num == 0) then
+            ! update estimate of u based on doubling number:
+            ! but not that not all were successful:
+            u_sum_old = u_mc_total * n_total
+            u_sum_new = u_mc * n_success
+            n_total = n_total + n_success
+            u_mc_total = (u_sum_old + u_sum_new) / n_total
+
+            ! relative error:
+            error = abs(u_mc_total - u_true) / abs(u_true)
+
+            write(25,'(i10,e23.15,e15.6)') n_total, u_mc_total, error
+            print '(i10,e23.15,e15.6)', n_total, u_mc_total, error
+            endif
+        n_mc = 2*n_mc
+        enddo
+
+    call MPI_REDUCE(nwalks,nwalks_total,1,MPI_INTEGER,MPI_SUM,0, &
+                 MPI_COMM_WORLD,ierr)
+
+    if (proc_num == 0) then
+        print '("Final approximation to u(x0,y0): ",es22.14)',u_mc_total
+        print '("Total walks performed by all processes: ",i10)', nwalks_total
+        endif
+
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
+
+    print '("Walks performed by Process ",i2,": ",i10)', &
+        proc_num,nwalks
+    
+    call MPI_FINALIZE(ierr)
+
+
+end program laplace_mc

File solutions/project/part4/mc_walk.f90

+
+module mc_walk
+
+implicit none
+integer :: nwalks
+
+contains
+
+
+subroutine random_walk(i0, j0, max_steps, ub, iabort)
+
+    use problem_description, only: nx,ny,ax,ay,dx,dy,uboundary
+
+    implicit none
+    integer, intent(in) :: i0,j0,max_steps
+    real(kind=8), intent(out) :: ub
+    integer, intent(out) :: iabort
+
+    real(kind=4), allocatable :: r(:)
+    integer i,j,istep
+    real(kind=4) :: rr
+    real(kind=8) :: xb,yb
+
+    allocate(r(max_steps))
+
+    i = i0
+    j = j0
+    call random_number(r)
+
+    do istep=1,max_steps
+        rr = r(istep)
+        if (rr < 0.25d0) then
+            i = i-1
+          else if (rr < 0.5d0) then
+            i = i+1
+          else if (rr < 0.75d0) then
+            j = j-1
+          else 
+            j = j+1
+          endif
+
+        if (i*j*(nx+1-i)*(ny+1-j) == 0) then
+            xb = ax + i*dx
+            yb = ay + j*dy
+            ub = uboundary(xb, yb)
+            exit  ! end the walk
+            endif
+        enddo
+
+    
+    if (istep == max_steps+1) then
+        iabort = 1
+        ub = 0.d0  ! this value won't be used
+    else
+        iabort = 0
+    endif
+
+    nwalks = nwalks + 1
+        
+end subroutine random_walk
+
+
+subroutine many_walks(i0, j0, max_steps, n_mc, u_mc, n_success)
+
+    use mpi
+
+    implicit none
+    integer, intent(in) :: i0,j0,max_steps,n_mc
+    integer, intent(out) :: n_success
+    real(kind=8), intent(out) :: u_mc
+
+    integer :: k, iabort, nb_sum
+    real(kind=8) :: ub,ub_sum
+
+    integer, dimension(MPI_STATUS_SIZE) :: status
+    integer :: num_procs, proc_num, ierr, j, numsent, sender
+
+    call MPI_COMM_SIZE(MPI_COMM_WORLD, num_procs, ierr)
+    call MPI_COMM_RANK(MPI_COMM_WORLD, proc_num, ierr)
+
+    ! -----------------------------------------
+    ! code for Master (Processor 0):
+    ! -----------------------------------------
+
+
+    if (proc_num == 0) then
+      ub_sum = 0.d0
+      nb_sum = 0
+
+      numsent = 0
+      ! send the first batch to get all workers working:
+      do j=1,min(num_procs-1, n_mc)
+          call MPI_SEND(MPI_BOTTOM, 0, MPI_INTEGER, j, 1, &
+                      MPI_COMM_WORLD, ierr)
+          numsent = numsent + 1
+          enddo
+
+      ! as results come back, send out more work...
+      ! the variable sender tells who sent back a result and ready for more
+      do k=1,n_mc
+        call MPI_RECV(ub, 1, MPI_DOUBLE_PRECISION, &
+                        MPI_ANY_SOURCE, MPI_ANY_TAG, &
+                        MPI_COMM_WORLD, status, ierr)
+        sender = status(MPI_SOURCE)
+        iabort = status(MPI_TAG)
+        if (iabort == 0) then
+            ub_sum = ub_sum + ub
+            nb_sum = nb_sum + 1
+            endif
+        if (numsent < n_mc) then
+            ! still more work to do
+            call MPI_SEND(MPI_BOTTOM, 0, MPI_INTEGER, sender, 1, &
+                      MPI_COMM_WORLD, ierr)
+            numsent = numsent + 1
+          else
+            ! send an empty message with tag=0 to indicate this worker
+            ! is done:
+            call MPI_SEND(MPI_BOTTOM, 0, MPI_INTEGER, &
+                            sender, 0, MPI_COMM_WORLD, ierr)
+          endif
+         enddo
+
+        u_mc = ub_sum / nb_sum
+        n_success = nb_sum
+        endif
+
+    ! -----------------------------------------
+    ! code for Workers (Processors 1, 2, ...):
+    ! -----------------------------------------
+    if (proc_num /= 0) then
+
+        if (proc_num > n_mc) go to 99   ! no work expected
+
+        do while (.true.)
+            ! repeat until message with tag==0 received...
+
+            call MPI_RECV(MPI_BOTTOM, 0, MPI_INTEGER, &
+                          0, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
+
+            j = status(MPI_TAG)   ! this tells worker whether to stop
+
+            if (j==0) go to 99    ! received "done" message
+
+            call random_walk(i0, j0, max_steps, ub, iabort)
+            call MPI_SEND(ub, 1, MPI_DOUBLE_PRECISION, &
+                        0, iabort, MPI_COMM_WORLD, ierr)
+
+            enddo
+        endif
+
+99  continue   ! jump here when done
+
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! wait for all processes to finish
+
+end subroutine many_walks
+
+
+end module mc_walk

File solutions/project/part4/plot_mc_laplace_error.py

+
+from pylab import *
+
+# read in three columns from file and unpack into 3 arrays:
+n,int_approx,error = loadtxt('mc_laplace_error.txt',unpack=True)
+
+figure(1)
+clf()
+loglog(n,error,'-o',label='Monte-Carlo')
+loglog([1,1e7],[1,sqrt(1e-7)],'k',label='1 / sqrt(N)')
+legend()
+xlabel('number of random walks taken')
+ylabel('abs(error)')
+title('Log-log plot of relative error in MC Laplace')
+savefig('mc_laplace_error.png')

File solutions/project/part4/problem_description.f90

+
+module problem_description
+
+    implicit none
+    real(kind=8), parameter :: ax = 0.0d0
+    real(kind=8), parameter :: bx = 1.0d0
+    real(kind=8), parameter :: ay = 0.4d0
+    real(kind=8), parameter :: by = 1.0d0
+    integer, parameter :: nx = 19
+    integer, parameter :: ny = 11
+    real(kind=8), parameter :: dx = (bx - ax) / (nx+1)
+    real(kind=8), parameter :: dy = (by - ay) / (ny+1)
+
+contains
+
+function utrue(x, y)
+
+    ! True solution for comparison, if known.
+
+    implicit none
+    real(kind=8), intent(in) :: x,y
+    real(kind=8) :: utrue
+
+    utrue = x**2 - y**2
+
+end function utrue
+
+function uboundary(x, y)
+
+    ! Return u(x,y) assuming (x,y) is a boundary point.
+
+    implicit none
+    real(kind=8), intent(in) :: x,y
+    real(kind=8) :: uboundary
+
+    if ((x-ax)*(x-bx)*(y-ay)*(y-by) .ne. 0.d0) then
+        print *, "*** Error -- called uboundary at non-boundary point"
+        stop
+        endif
+
+    uboundary = utrue(x,y)   ! assuming we know this
+
+end function uboundary
+    
+
+end module problem_description

File solutions/project/part4/random_util.f90

+
+module random_util
+
+contains
+
+subroutine init_random_seed(seed1)
+
+    ! Seed the random number generator.
+    ! If seed1 = 0 then set seed1 randomly using the system clock.
+    !              This will give different sequences of random numbers
+    !              from different runs, so results are not reproducible.
+    ! If seed1 > 0 then results will be reproducible since calling this
+    !              twice with the same seed will initialize the random
+    !              number generator so the same sequence is generated.
+    ! Once seed1 is set, set the other elements of the seed array by adding
+    ! multiples of 37 as suggested in the documentation.
+    ! The length of the seed array is determined by calling  
+    !    
+    integer, intent(inout) :: seed1
+
+    integer :: nseed, clock 
+    integer, dimension(:), allocatable :: seed
+
+    call random_seed(size = nseed)  ! determine how many numbers needed to seed
+    allocate(seed(nseed)) 
+
+    if (seed1 == 0) then
+        ! randomize the seed: not repeatable
+        call system_clock(count = clock)
+        seed1 = clock
+      endif
+
+    do i=1,nseed
+        seed(i) = seed1 + 37*(i-1)  
+        enddo
+
+    print *, "seed1 for random number generator:", seed1
+
+    call random_seed(put = seed)   ! seed the generator
+    deallocate(seed)
+
+end subroutine init_random_seed
+
+end module random_util