\documentclass{article}
\usepackage{html}
\begin{document}
\title{The LAM companion to ``Using MPI...''}
\author{Zdzislaw Meglicki\\
        Centre for Information Science Research\\
        The Australian National University\\
        Canberra, ACT, 0200, Australia\\
        phone: +61-6-249-0158\\
        fax: +61-6-249-0747\\
        e-mail: Zdzislaw.Meglicki@cisr.anu.edu.au}

\maketitle

\section{About this tutorial}

   This tutorial is based on an excellent book by William Gropp, Ewing Lusk, and
   Anthony Skjellum, ``Using MPI -- Portable Parallel Programming with the
   Message-Passing Interface''. The book was published by The MIT Press in
   1994, ISBN 0-262-57104-8, and it costs \$US24.95, if bought directly from
   the publisher, whom you can contact by electronic mail.

   However, you should be able to go through the tutorial even if you don't
   have the book. The tutorial is self-contained and references are provided
   to the basic MPI definition document which you can obtain by ftp anonymous. 

   \subsection{MPICH versus LAM}

   ``Using MPI...'' is based on an MPI implementation developed 
   jointly by Argonne
   National Laboratory and Mississippi State University called MPICH. The ``CH''
   ending derives from ``Chameleon'', a system which MPICH is based on. MPICH
   is really quite a towering construction. It can be built either on top of
   PVM or on top of P4. Once you have any of the systems compiled,
   you build ``Chameleon'' on top of it, and then you build MPI
   library on top of ``Chameleon''. 

   This complex installation procedure is fully automated and, in particular if
   you decide to base MPICH on P4, it can be achieved just with 
   one {\tt make} command. The P4 implementation of MPICH is actually very
   fast, but it has one flaw. The startup time of MPICH/P4 applications on a
   workstation farm can be very long, because P4 spawns jobs on nodes 
   other than the console by {\tt rsh}. If your system does a lot of DES based
   checking and automounting, you may have to wait up to several minutes for
   your application to come up.

   This can be a nerve wrecking experience if you are developing a code or
   learning about MPI.

   For this reason this course is based on the LAM implementation of MPI from
   Ohio Supercomputer Center. LAM, like PVM, is based on a network of daemons
   running on a farm. It takes a while to start up the LAM engine, but once it
   is running, the applications come up very quickly: the daemons are there
   to pick them up and run. For this reason LAM is good for application 
   development and learning. LAM has many other advantages and the LAM 
   based MPI implementation is stable and thorough. Section,
   \ref{s:tut-and-book}, describes how to
   obtain and install LAM and how to operate the ``LAM engine''. 

   \subsection{The tutorial and the book}
   \label{s:tut-and-book}

      The tutorial itself covers examples from chapters 3 through 5 of
      ``Using MPI...''. I also provide detailed references to the MPI
      definition document about which more below. I have tweaked the examples
      so as to get them going on the ``LAM engine''. Some examples had to be
      modified more than others. In particular Argonne MPICH {\em assumes\/}
      that the console process always has rank 0. But this is not specified in
      the MPI standard, and LAM does it differently. This has an impact on how
      you do I/O. These issues have been left out of the MPI definition because
     they are implementation and vendor dependent. For example on a
      workstation farm you can always rely on standard {\sc unix} I/O. Even if
      you cannot write directly to your VDU, you can always open a file in
      {\tt /tmp} or {\tt /usr/tmp} and write there. But this may not be
      possible if you run your MPI programs on a machine which does not support
     full {\sc unix} on every node.
 
      Lack of I/O specification is a shortcoming of MPI library, acknowledged
      in ``Using MPI...'' (section 10.4, page 245) and in the MPI definition
      document (section 1.5, page 3). But this shortcoming is at present
      insurmountable due to the fact that the idea of installing a full
      {\sc unix} I/O on every node has not been universally
      accepted by all supercomputer vendors.
 
      Another specific feature of this tutorial is that it is restricted to
      pure MPI. You will find quite a few references to an extension package
      available with MPICH in ``Using MPI...''. Those extensions, which are
      enclosed in the MPE library, are not necessarily supported by other
      MPI implementations. LAM specific features, such as using the LAM frame
      buffer and LAM LED facility will be covered in the companion to this
      tutorial later.

   \subsection{The MPI definition}

      The MPI standard has been defined in
      ``MPI: A Message-Passing Interface
      Standard''. That document was produced by the
      Message Passing Interface Forum,
      and is available from the World Wide Web server at
      the University of Tennessee at Knoxville,
      {\tt http://netlib.cs.utk.edu/}. When you will have made
      the connection, go to the
      Netlib Repository, and then select {\tt mpi}.
      The file to transfer is {\tt mpi/mpi-report.ps}, which contains the MPI
      Standard specification as of May 5, 1994. The report is 228 pages long.
      When printing the report you may try to put 4 pages on one sheet of paper
     using the {\tt multi.ps} package for PostScript from Ross Cartlidge,
      {\tt rossc@extro.ucc.su.oz.au}.
 
      This document is a must, if you are really serious about using MPI.
      It is quite readable and contains numerous examples and explanations.

      MPI Standard document is being converted to html and it is already 
      available from {\tt www.mcs.anl.gov}. If you press
      \htmladdnormallink{here}
      {http://www.mcs.anl.gov/mpi/mpi-report/mpi-report.html}, the magic
      of Mosaic will take you all the way from Canberra in Australia 
      (or from wherever you read this document now) to
      Argonne near Chicago in the US, and you will be able to view the
      MPI Standard document on line.

 
\section{Installing, configuring, and operating the LAM engine}

If you already have LAM installed and configured on your system, and if you
know how to operate it, you can skip right to section \ref{s:mpi-lam}, and
then go to the first lesson. If you don't have LAM or if you haven't used
it before, read on. The remainder of this section provides you with sufficient
pointers, I hope, to obtain and install LAM. It also tells you how to operate
the LAM engine. Only the minimum of LAM related knowledge is covered 
necessary to run simple MPI jobs, so that we can focus on MPI itself as quickly
as possible.

   \subsection{Obtaining LAM}

      What to do if you don't have LAM? LAM is available by ftp anonymous from
      {\tt tbag.osc.edu}, directory {\tt /pub/lam}. You can find out 
      more about LAM
      by accessing the World Wide Web server at the Ohio Supercomputer 
      Center, {\tt http://www.osc.edu/}.

   \subsection{Installation of LAM}

      Once you have the software transferred, extract the documentation and
      read the installation manual. You should be able to find it in 
      {\tt ~lam52/doc} subdirectory.

      A word of warning. The document is in PostScript and it requires that you
      have an American ``letter'' tray in your printer. I could not print it on
      Apple Laser Printers with A4 trays. I also had some problems viewing it
      with ghostview, but finally managed to read this and other LAM documents
      using pageview (all on a Sun running Solaris 2.3).

      The installation itself is quite trivial. There are 12 patches to apply
      to the source, then configure one of the files and type {\tt make}.
      It's all explained in the installation manual which is only 10 pages
      long. 

      LAM is available for Sun 4 running SunOS and Solaris, SGI, IBM RS6000,
      Alphas running OSF/1, and HPs.

      If you haven't done any installations before, and if you have a system
      administrator, ask her to do it for you. I'm sure she'll be delighted to
      be of assistance. You may have to have a large bar of chocolate ready.

   \subsection{Configuring LAM}

      You should again read LAM documentation about configuring your farm. I
      went for the simplest solution. My LAM is installed under Solaris 2.3
      in {\tt /opt/lam}. There is a subdirectory there 
      called {\tt /opt/lam/boot}. In that directory there is a file 
      {\tt bhost.def}, and that file looks as follows
\begin{verbatim}
gustav@cisr:/opt/lam/boot 236 $ cat bhost.def
#
#       Ohio Trollius
#
#       Function:       - bhost.def
#                       - CISR LAM system
#
cisr
mercury
venus
earth
mars
jupiter
saturn
uranus
bacchus
ceres
diana
minerva
vesta
gustav@cisr:/opt/lam/boot 237 $ 
\end{verbatim}
      It simply lists the names of machines in our Solaris farm. That's all you
      need to do. 

      You can have a more elaborate set up with weights and topologies. This is
      also described in LAM documents.

   \subsection{Starting up the LAM engine}

   Assuming that your LAM machine has been installed and properly configured,
   the operation of the system is quite easy. Before you start up the engine,
   issue the command {\tt recon}:
\begin{verbatim}
gustav@cisr:~/mpi/tutorial/mpi 259 $ recon -v
recon: testing n0 (cisr)
recon: testing n1 (mercury)
recon: testing n2 (venus)
recon: testing n3 (earth)
recon: testing n4 (mars)
recon: testing n5 (jupiter)
recon: testing n6 (saturn)
recon: testing n7 (uranus)
recon: testing n8 (bacchus)
recon: testing n9 (ceres)
recon: testing n10 (diana)
recon: testing n11 (minerva)
recon: testing n12 (vesta)
gustav@cisr:~/mpi/tutorial/mpi 260 $ 
\end{verbatim}
   This command checks if every node of the LAM engine
   is correctly configured and available. If {\tt recon} complains that
   something is wrong, don't try to start LAM. Consult your system
   administrator, if you can't fix the problem yourself. 

   Assuming that everything is fine, start the LAM engine by issuing the
   command {\tt lamboot}:
\begin{verbatim}
gustav@jupiter:~ 201 $ lamboot -v
 
LAM 5.2 - Ohio Supercomputer Center
  
hboot n0 (cisr)...
hboot n1 (mercury)...
hboot n2 (venus)...
hboot n3 (earth)...
hboot n4 (mars)...
hboot n5 (jupiter)...
hboot n6 (saturn)...
hboot n7 (uranus)...
hboot n8 (bacchus)...
hboot n9 (ceres)...
hboot n10 (diana)...
hboot n11 (minerva)...
hboot n12 (vesta)...
   
topology
n0...n1...n2...n3...n4...n5...n6...n7...n8...n9...n10...n11...n12...
gustav@jupiter:~ 202 $ 
\end{verbatim}
   This will take a long while, but once LAM is up, normal MPI applications
   will start up without such delay.

   \subsection{Testing and inspecting LAM}

   With the LAM engine running, you can inspect it in various ways. For example
   you can {\tt tping} all nodes:
\begin{verbatim}
gustav@jupiter:~ 202 $ tping -c 3 N
  1 byte from 13 nodes: 0.434 secs
  1 byte from 13 nodes: 0.213 secs
  1 byte from 13 nodes: 0.213 secs
 
3 messages, 3 bytes (0.003K), 0.860 secs (0.007K/sec)
roundtrip min/avg/max: 0.213/0.287/0.434
gustav@jupiter:~ 203 $ 
\end{verbatim}
   You can also check what runs on the nodes by issuing the command {\tt
   state}:
\begin{verbatim}
gustav@jupiter:~ 203 $ state N
NODE      INDEX  PID      KPRI     KSTATE               PROGRAM
gustav@jupiter:~ 204 $ 
\end{verbatim}

   The command line argument {\tt N} stands for {\em all nodes\/}. You can
   inspect only selected nodes with LAM commands. For example:
\begin{verbatim}
gustav@jupiter:~ 204 $ tping -c 3 n0-5
  1 byte from 6 nodes: 0.190 secs
  1 byte from 6 nodes: 0.098 secs
  1 byte from 6 nodes: 0.132 secs
 
3 messages, 3 bytes (0.003K), 0.420 secs (0.014K/sec)
roundtrip min/avg/max: 0.098/0.140/0.190
gustav@jupiter:~ 205 $ 
\end{verbatim}

   \subsection{Recovering from LAM crashes}

   If at any stage your LAM program hangs or crashes, you should clean the
   LAM machine. The command to do that is {\tt lamclean}:
\begin{verbatim}
gustav@jupiter:~ 205 $ lamclean -v
killing processes: 0...1...2...3...4...5...6...7...8...9...10...11...12...
sweeping messages: 0...1...2...3...4...5...6...7...8...9...10...11...12...
closing files: 0...1...2...3...4...5...6...7...8...9...10...11...12...
removing groups...
gustav@jupiter:~ 206 $ 
\end{verbatim}

   \subsection{Running MPI jobs under LAM}
   \label{s:mpi-lam}

   To run an MPI job you must use the command {\tt mpirun}. I will explain that
   command later. Here you can see the example of its usage on a simple 
   MPI program which will be our first example:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 213 $ mpirun -w N `pwd`/first -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
gustav@jupiter:~/mpi/simplempi 214 $ 
\end{verbatim}

   \subsection{Shutting down the LAM engine}

   Finally, how to terminate your LAM session? It is not enough just to log out.
   LAM daemons will be left running in the background. You should shutdown the
   LAM engine with the command {\tt wipe}:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 214 $ wipe -v
Enter boot schema filename (bhost format): /opt/lam/boot/bhost.def
tkill n0 (cisr)...
tkill n1 (mercury)...
tkill n2 (venus)...
tkill n3 (earth)...
tkill n4 (mars)...
tkill n5 (jupiter)...
tkill n6 (saturn)...
tkill n7 (uranus)...
tkill n8 (bacchus)...
tkill n9 (ceres)...
tkill n10 (diana)...
tkill n11 (minerva)...
tkill n12 (vesta)...
gustav@jupiter:~/mpi/simplempi 215 $ 
\end{verbatim}
   Observe that in this one case you will have to specify the name of the boot
   schema file explicitly. All other LAM commands 
   read {\tt /opt/lam/boot/bhost.def} by default.

\section{Lesson 1 -- exploring MPI environment, sending and receiving messages}

   Our first MPI example is already quite involved, but this is a hurdle we
   must overcome right at the beginning in order to be able to write messages
   to the display in a reasonably portable manner. This program is shown in
   the {\em verbatim\/} listing in the next section, \ref{s:first}. 
   Like all other examples in this tutorial
   you can obtain it by ftp anonymous from {\tt cisr.anu.edu.au}, directory
   {\tt /pub/papers/meglicki/mpi/simplempi}. The file is called {\tt first.c}.

   This program does not come from the Argonne ``Using MPI...'' book and you
   will not find it in the Argonne ftp-anonymous directories. Their first
   program is much simpler than this one, but it makes certain assumptions
   about I/O, which are valid for MPICH, but which are not valid for LAM, or
   some other MPI implementations. So we must take this detour before we can
   merge with ``Using MPI...''.

   \subsection{The listing of {\tt first.c}}
   \label{s:first}

\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include "mpi.h"

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

    MPI_Init ( &argc, &argv );
    MPI_Comm_size ( MPI_COMM_WORLD, &pool_size ); 
    MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
    MPI_Get_processor_name ( my_node_name, &node_name_length );
    MPI_Attr_get ( MPI_COMM_WORLD, MPI_HOST, (void**) &host_rank_ptr,
                   &found_flag );
    if ( (int) host_rank_ptr != MPI_PROC_NULL ) host_rank = *host_rank_ptr; 
    else {
       MPI_Attr_get ( MPI_COMM_WORLD, MPI_IO, (void **) &host_rank_ptr,
                      &found_flag );
       if ( (int) host_rank_ptr != MPI_ANY_SOURCE ) host_rank = *host_rank_ptr;
       else {
          strcpy ( console_name, argv[1] );
          if ( 0 == strcmp ( my_node_name, console_name ) ) {
             host_rank = my_rank;
             i_am_the_master = B_TRUE;
          }
       }
    }
    if (i_am_the_master) {
       int_buffer[0] = host_rank;
       for (destination = 0; destination < pool_size; destination++)
          if (destination != my_rank) 
             MPI_Send(int_buffer, 1, MPI_INT, destination, 1001, 
                      MPI_COMM_WORLD);
    }
    else {
       MPI_Recv(int_buffer, 1, MPI_INT, MPI_ANY_SOURCE, 1001, 
                MPI_COMM_WORLD, &status);
       host_rank = int_buffer[0];
    }
    if (i_am_the_master) strcpy (char_buffer, my_node_name);
    MPI_Bcast(char_buffer, BUFSIZ, MPI_CHAR, host_rank, 
              MPI_COMM_WORLD);
    strcpy (host_name, char_buffer);
    if (! i_am_the_master) {
       sprintf (char_buffer, "Greetings to the master (%s, %d) from (%s, %d)",
                host_name, host_rank, my_node_name, my_rank);
       MPI_Send (char_buffer, strlen(char_buffer) + 1, MPI_CHAR, host_rank, 
                 2002, MPI_COMM_WORLD);
    }
    else {
       for (source = 0; source < pool_size; source++) {
          if (source != my_rank) {
             MPI_Recv(char_buffer, BUFSIZ, MPI_CHAR, source, 2002,
                      MPI_COMM_WORLD, &status);
             printf ("%s\n", char_buffer);
          }
       }
    }
    MPI_Finalize ();
}
\end{verbatim}

   \subsection{Compiling and running {\tt first.c} under LAM}

      I now assume that you have {\tt first.c} in your directory.

      Before we go through this program in detail, let us compile it and 
      run first.  To compile {\tt first.c} simply say:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 240 $ hcc -o first first.c -lmpi
gustav@jupiter:~/mpi/simplempi 241 $ 
\end{verbatim}

      {\tt hcc} is a convenience wrapper 
      for the local native C compiler. It is a part of LAM system.
      On our systems it seems to prefer GNU C.
      It knows about the location of LAM libraries which is why you don't have
      to specify {\tt -L/opt/lam/lib}. But note that if you compile an MPI
      program you {\em do\/} have to specify {\tt -lmpi}.

      To run this program invoke it through the {\tt mpirun} command. For
      example:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 243 $ mpirun -w N `pwd`/first -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
gustav@jupiter:~/mpi/simplempi 244 $ 
\end{verbatim}

      The {\tt -w} switch tells {\tt mpirun} not to exit before all MPI
      processes complete. The capital {\tt N} tells {\tt mpirun} to run the job
      on all nodes of the LAM engine. 
      We have to specify the full pathname of the binary if it is
      not located on the normal command search path. This pathname must be
      identical on all machines running LAM, unless a {\tt -s} switch is used,
      in which case the binary is loaded from your local node and distributed
      to other nodes through LAM channels. Finally we have the double hyphen,
      \verb|--|, which 
      marks the beginning of command line arguments referring 
      to our MPI program. The only argument we pass to our first
      MPI program is the name of the machine
      from which the program is invoked. This name is generated dynamically by
      {\sc unix} command {\tt uname -n}.

      Like with all other LAM commands you can obtain more information about
      {\tt mpirun} and {\tt hcc} by issuing the commands
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 244 $ man mpirun
\end{verbatim}
      and
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 245 $ man hcc
\end{verbatim}

   \subsection{The synopsis of {\tt first.c}}

      Program {\tt first} demonstrates the use of MPI inquiry functions,
      and the use of point to point communications as well as the use of
      broadcast. The purpose of the activities within this program
      is to configure one of the
      nodes as the one through which all messages will be written to 
      your VDU. 

      Text written using {\tt printf} on nodes other than the node from
      which LAM was started are sent to {\tt /dev/null}. Only text printed
      on the {\em local\/} node will appear on your display. This behaviour is
      LAM specific. MPICH, for example, channels output from all nodes to
      the window from which the job was started. PVM based UNIFY, on the other
      hand, writes standard output to special files 
      which PVM creates automatically in
      {\tt /tmp}. On other systems some nodes may not support any I/O at all,
      and there may be one predefined {\em host\/} node.

      Our job is to find out if there is a predefined host node and
      redirect all information from other nodes
      to that host, if there is.
      If there is no host, we should check which nodes can 
      do I/O at all, and select one of those nodes as our host. 
      Finally, if there is no host and if all processes can do I/O we
      have to determine our host in some other way. For example we can
      take the name of the host from the command line. The possibility that
      there is no host and that no processes can do I/O is too horrible to
      contemplate. A system like that would be useless anyway. 

      Once we have determined, which of the processes is going 
      to channel the I/O,
      we have to inform all other processes about our choice. When other
      processes know where to send the information they send homage to the
      master process, which redirects the messages to {\tt stdout}.

   \subsection{{\tt first.c} in detail}
   \label{s:detail}

      Let us go through {\tt first.c} step by step. The program begins with an
      opening
\begin{verbatim}
MPI_Init ( &argc, &argv );
\end{verbatim}
      and ends with
\begin{verbatim}
MPI_Finalize ();
\end{verbatim}

      All MPI programs must do that. You can have some other stuff before
      \verb|MPI_Init| and after \verb|MPI_Finalize|, but that must not be MPI.
      You should not assume that your job becomes parallel only after you issue
      the \verb|MPI_Init| call though. What happens before \verb|MPI_Init| and
      after \verb|MPI_Finalize| is system dependent. LAM will start all
      processes on the parallel machine from the very beginning.

      \subsubsection{Environmental enquiries}
      \label{s:environ}

      After MPI has been initialised with \verb|MPI_Init| we find out about the
      number of processes available for computation. That number is returned in
      \verb|pool_size| after the call to 
\begin{verbatim}
MPI_Comm_size ( MPI_COMM_WORLD, &pool_size ); 
\end{verbatim}

      The first argument to \verb|MPI_Comm_size| is \verb|MPI_COMM_WORLD|. This
      is the default {\em communicator,\/} i.e., a group of processes which
      comprises all processes in the pool. All MPI operations occur within the
      communicators. We will discuss this concept in more detail later. 

      Now, note that every process carries out all these operations on its own.
      So that when the next enquiry is made:
\begin{verbatim}
MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
\end{verbatim}
      each process ends up with a different number in \verb|my_rank|. That
      number is the process identity within the communicator 
      \verb|MPI_COMM_WORLD|.
      Processes can have several rank numbers corresponding to different
      communicators. The rank number allows processes to find out which is
      which, and what is their own place in the pool.

      Those three calls
\begin{verbatim}
MPI_Init ( &argc, &argv );
MPI_Comm_size ( MPI_COMM_WORLD, &pool_size ); 
MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
\end{verbatim}
      are a standard beginning of every MPI program. They are not wrapped up
      into one call though, because the functions \verb|MPI_Comm_size| and
      \verb|MPI_Comm_rank| may have to be called again, whenever you create a
      new {\em communicator}.

      In MPI each processor must also have a name. On a farm this name is
      simply the host name of the machine the process runs on. On a tightly
      coupled multiprocessor, such as Cray T3D or the Connection Machine, it
      may be the machine number of the processor, for example 
      something like \verb|node_17|. To find out about the name of the
      processor the function \verb|MPI_Get_processor_name| is called:
\begin{verbatim}
MPI_Get_processor_name ( my_node_name, &node_name_length );
\end{verbatim}
      We must reserve appropriate amount of space for \verb|my_node_name|.
      \verb|BUFSIZ|, which is either 512 bytes or 1024 bytes should do for even
      the most extravagantly named machines.

      Now we check if our process pool has a predefined host. This is done with
      the call
\begin{verbatim}
MPI_Attr_get ( MPI_COMM_WORLD, MPI_HOST, (void**) &host_rank_ptr,
               &found_flag );
\end{verbatim}
      This function checks if an attribute defined by \verb|MPI_HOST| exists.
      If it does, \verb|found_flag| is set to 1, otherwise \verb|found_flag| is
      zero. Finding that the attribute exists, does not mean that there is a
      host. The answer may still be \verb|MPI_PROC_NULL|. Here we don't
      inspect the flag, because the attribute \verb|MPI_HOST| is one of the 
      three attributes which {\em must\/} be attached 
      to the \verb|MPI_COMM_WORLD| communicator when MPI is initialised. 
      The other two attributes are
      \verb|MPI_TAG_UB| and \verb|MPI_IO|. The first one refers to the upper
      bound for the message tag value (more about message tags soon). 
      In principle we should check that too, but the
      upper bound {\em must\/} be at least 32,767, so as long as our message
      tags are going to be below that, we're safe. The second attribute, 
      \verb|MPI_IO| will tell us which processes can support I/O.

      When the function \verb|MPI_Attr_get| returns, we first check what the
      returned pointer \verb|host_rank_ptr| points to. If it points to 
      \verb|MPI_PROC_NULL| it means that there is no host. Otherwise there is a
      host, and its rank number within the \verb|MPI_COMM_WORLD| can be found
      by inspecting the first entry in the \verb|host_rank_ptr| array.

      The MPI environmental management functions are very important, yet they
      seem to have been neglected by many early implementations of MPI to the
      extent that some of them used to return a dangling pointer which would
      crash the system the moment you asked
\begin{verbatim}
if ( (int) host_rank_ptr != MPI_PROC_NULL ) ...
\end{verbatim}
      You can find more about MPI environmental management in 
      \htmladdnormallink{Chapter 7}
      {http://www.mcs.anl.gov/mpi/mpi-report/node140.html#Node140}
      (page 189) of ``MPI: A Message-Passing Interface Standard''. The function
      \verb|MPI_Attr_get| is described on page 169 of the Standard in section
      5.7.

      In the case of the LAM engine there will not be a predefined host. So,
      we still have to continue our exploration of the MPI environment.

      We continue by checking which processes can support I/O. We repeat the
      call to \verb|MPI_Attr_get|, but the attribute this time is
      \verb|MPI_IO|.
\begin{verbatim}
MPI_Attr_get ( MPI_COMM_WORLD, MPI_IO, (void **) &host_rank_ptr,
               &found_flag );
\end{verbatim}
      When the call comes back we again check what \verb|host_rank_ptr| points
      to. As we remarked in the synopsis \verb|MPI_PROC_NULL| is too terrible
      to contemplate, but it may turn out, as indeed it will for {\sc
      unix} farms, that {\em all\/} processes can do I/O. In this
      case \verb|host_rank_ptr| will point to \verb|MPI_ANY_SOURCE|. 

      If only some
      processes can do I/O, there is little choice left. The ranks of those
      processes will be returned in the array \verb|host_rank_ptr|, and 
      in \verb|first.c| we choose the first of those processes. That will
      work, for example, under MPICH, because under MPICH process number 0 is
      always the one talking to your VDU:
\begin{verbatim}
if ( (int) host_rank_ptr != MPI_ANY_SOURCE ) host_rank = *host_rank_ptr;
\end{verbatim}

      In general you will have to be very careful at this stage. 
      This first process may not be the one which has access to your VDU.
      In LAM, for example, a process running on the machine from which LAM was
      started talks to the VDU, and that process will have a rank corresponding
      to the position that machine has in /opt/lam/boot/bhost.def.

      If the returned pointer is \verb|MPI_ANY_SOURCE|, it means that all
      processes can do I/O, although not necessarily to your VDU. In this case
      \verb|first.c| gives up and simply transfers the name of the machine
      from the command line to \verb|console_name|:
\begin{verbatim}
else {
   strcpy ( console_name, argv[1] );
   if ( 0 == strcmp ( my_node_name, console_name ) ) {
      host_rank = my_rank;
      i_am_the_master = B_TRUE;
   }
}
\end{verbatim}
      Now every process compares that name to its own name,
      \verb|my_node_name|, and the one for which this condition is satisfied
      sets \verb|i_am_the_master| to \verb|B_TRUE|, whereas for all other
      processes it remains, as it had been initialised originally,  
      \verb|B_FALSE|.

      The procedure described in this section is by no means perfect.
      I can imagine several situations in which it will fail. If you work with
      LAM you may be tempted to skip all this mess with 
      \verb|MPI_Attr_get| and take the name of the node from the command line
      at the very beginning. If you work with MPICH you know that the process
      you are looking for will always have rank 0.

      But knowing about environmental enquiries you can also design quite an
      elaborate procedure which will check for the availability of predefined
      hosts, I/O, etc, and make appropriate choices using whatever other
      information it can obtain.

      \subsubsection{Sending and receiving messages}
      \label{s:send-receive}

      At this stage the situation is as follows: one of the processes 
      already knows who is the master, namely the master process itself.
      But no other process knows that. The master must inform other
      processes about the choice that has been made. The master does it by
      sending a message to every other process:
\begin{verbatim}
if (i_am_the_master) {
   int_buffer[0] = host_rank;
   for (destination = 0; destination < pool_size; destination++)
      if (destination != my_rank) 
         MPI_Send(int_buffer, 1, MPI_INT, destination, 1001, 
                  MPI_COMM_WORLD);
}
\end{verbatim}

      In order to send a message you must have a preallocated buffer. 
      In our case it is
      an integer array \verb|int_buffer|. The contents of the message is
      transferred to the buffer, and the buffer is passed to 
      \verb|MPI_Send|. The second argument in the call to 
      \verb|MPI_Send| is the number of items in the buffer and the third
      argument is the type of those items. You must not use C types here. MPI
      provides its own constants, such as 
      \verb|MPI_INT|, 
      \verb|MPI_CHAR|, 
      \verb|MPI_DOUBLE|, which
      specify the type of transmitted data. The fourth argument is the rank of
      the distination process within the communicator which is specified by the
      sixth argument. The fifth argument is the message tag. This is a specific
      number which is assigned to each message. These numbers can be used for
      ordering or narrowing the scope of messages. 

      When messages arrive at the addressee, they will be read only if the
      addressee wants to access the message with the same tag number.

      The function \verb|MPI_Send| is described in section 
      \htmladdnormallink{3.2.1}
      {http://www.mcs.anl.gov/mpi/mpi-report/node30.html#Node30}, 
      page 16, of ``MPI: A Message-Passing Interface Standard''.

      The part of the code which refers to receiving processes is shown below:
\begin{verbatim}
else {
   MPI_Recv(int_buffer, 1, MPI_INT, MPI_ANY_SOURCE, 1001, 
            MPI_COMM_WORLD, &status);
   host_rank = int_buffer[0];
}
\end{verbatim}
      
      In order to receive a message, you also must have a buffer into which
      the message will be transferred. You must specify the length of the
      buffer in the second argument to \verb|MPI_Recv| and the type of expected
      data item in the third argument. 
      The fourth argument specifies the rank of the process from
      which the message is to come. In this case we don't care where the
      message comes from, as long as it has a tag 1001, and as long as it
      comes within the \verb|MPI_COMM_WORLD| communicator.
      The constant \verb|MPI_ANY_SOURCE| is used to convey this
      indifference to \verb|MPI_Recv|. The last argument is a pointer to 
      \verb|status|, which is a structure of type \verb|MPI_Status|.
      Special MPI functions are available to read that structure, see
      section \ref{s:fourth-master-part}.

      The function \verb|MPI_Recv| is described in section 
      \htmladdnormallink{3.2.4}
      {http://www.mcs.anl.gov/mpi/mpi-report/node33.html#Node33}, 
      page 19, of ``MPI: A Message-Passing Interface Standard''.

      After all processes with the exception of the sender
      have received the message, the whole process pool knows
      the rank of the master process.

      The functions \verb|MPI_Send| and \verb|MPI_Recv| are fundamental.
      With these two functions and with environmental enquiry facilities
      covered in section \ref{s:environ}, you can already build serious
      applications. John Salmon's famous N-body octal tree code, for example,
      has been designed to use message 
      {\em send\/} and {\em receive\/} facilities
      {\em only\/} in order to ensure easy portability of
      the program between various message passing systems.

      This section also illustrates how you can make various processes
      execute different parts of the code depending on their rank.
      This is a very common situation in message passing programs.

      \subsubsection{Collective communication, broadcast}
      \label{s:broadcast}

      Once all processes within the \verb|MPI_COMM_WORLD| communicator
      know the rank of the master process we can attempt a collective
      communication. This is accomplished in the following three lines of the
      code:
\begin{verbatim}
if (i_am_the_master) strcpy (char_buffer, my_node_name);
MPI_Bcast(char_buffer, BUFSIZ, MPI_CHAR, host_rank,
          MPI_COMM_WORLD);
strcpy (host_name, char_buffer);
\end{verbatim}
      The call to function \verb|MPI_Bcast| is issued by all processes.
      But in this operation only one process defined by the fourth argument,
      i.e., the process whose rank is \verb|host_rank|, sends the 
      information originally contained in \verb|char_buffer| of that
      process. All other processes receive that information into their
      own \verb|char_buffer| arrays. The length of that buffer is defined by
      the second argument, and the type of the data to be received is defined
      by the third argument.
      The broadcast is restricted to the \verb|MPI_COMM_WORLD| communicator.

      Note that it is dangerous to use here something like
\begin{verbatim}
strlen(char_buffer) + 1
\end{verbatim}
      in place of BUFSIZ, because only on the master process will 
      \verb|strlen(char_buffer)| return anything meaningful.

      Note also that we could not use broadcast until all processes knew
      the rank of the host process. For this reason we had to communicate that
      rank using \verb|MPI_Send|.

      Collective communications are described in 
      \htmladdnormallink{Chapter 4} 
      {http://www.mcs.anl.gov/mpi/mpi-report/node62.html#Node62}
      of ``MPI: A
      Message-Passing Interface Standard''. And \verb|MPI_Bcast| is described
      on page 93, section 
      \htmladdnormallink{4.4}
      {http://www.mcs.anl.gov/mpi/mpi-report/node66.html#Node66}.

      In ``Using MPI...'' broadcast is the first communication operation
      covered by the first example in Chapter 3, page 23. This can be done
      under MPICH because the host process has a rank number 0 by default and
      all other processes know that from the beginning.

      \subsubsection{Paying homage to the master}

      The last part of the code does not introduce any new elements.
      It is basically a reversal of the first {\em send-receive} operation.
      This time the slave processes send some information about themselves
      to the master, and the master having received that information 
      prints it on standard output.
\begin{verbatim}
if (! i_am_the_master) {
   sprintf (char_buffer, "Greetings to the master (%s, %d) from (%s, %d)",
            host_name, host_rank, my_node_name, my_rank);
   MPI_Send (char_buffer, strlen(char_buffer) + 1, MPI_CHAR, host_rank, 
             2002, MPI_COMM_WORLD);
}
else {
   for (source = 0; source < pool_size; source++) {
      if (source != my_rank) {
         MPI_Recv(char_buffer, BUFSIZ, MPI_CHAR, source, 2002,
                  MPI_COMM_WORLD, &status);
         printf ("%s\n", char_buffer);
      }
   }
}
\end{verbatim}
      Observe that in this case the host process receives the messages 
      in a specific
      order. The order corresponds to the rank of the sender. If the messages
      arrive out of order, they will be spooled until the reading process
      is ready to accept them. 

      \subsection{Summary: the end of the detour}

      In this first MPI lesson we have covered a large number of MPI
      issues. We will summarise by listing MPI functions used in
      {\tt first.c}:

      \begin{itemize}
         \item \verb|MPI_Init| initialises the MPI system; explained in
            section \ref{s:detail}
         \item \verb|MPI_Finalize| closes the MPI system; explained in section
            \ref{s:detail}
         \item \verb|MPI_Comm_size| finds the number of processes in the pool;
            explained in section \ref{s:environ}
         \item \verb|MPI_Comm_rank| finds the rank number of a calling process;
            explained in section \ref{s:environ}
         \item \verb|MPI_Get_processor_name| finds the name of a processor the
            calling process runs on; explained in section \ref{s:environ}
         \item \verb|MPI_Attr_get| obtains the value of a communicator
            attribute, for example, the host process rank or the ranks of
            processes which can do I/O; explained in section \ref{s:environ}
         \item \verb|MPI_Send| sends a message; explained in section
            \ref{s:send-receive}
         \item \verb|MPI_Recv| receives a message; explained in section
            \ref{s:send-receive}
         \item \verb|MPI_Bcast| broadcasts a message; explained in section
            \ref{s:broadcast}
      \end{itemize}

      In following lessons all example programs will begin by peforming the
      operations described in {\tt first.c}. We will be including a file
      which contains those operations using the cpp \verb|#include|
      directive, so as to avoid cluttering the
      code and the lessons with a repetitive material that has been covered
      in this first lesson. See section \ref{s:preamble} for the listing
      of that file.

\section{Lesson 2 -- collective communications}

   Now we are almost ready to merge with ``Using MPI...''. In all our future
   programs we will include as a preamble the body of the program developed 
   in the first section. This insert is called {\tt preamble.c} and, like all
   other programs of this tutorial, it can be obtained by ftp-anonymous from
   {\tt cisr.anu.edu.au}, directory {\tt pub/papers/meglicki/mpi/simplempi},
   file {\tt preamble.c}.

   We also provide a listing of {\tt preamble.c} 
   in the next subsection, in case you cannot 
   access our ftp area. There is nothing new in this listing, so you can skip 
   it, if you already have {\tt preamble.c} transferred to your directory.

   In this lesson we will calculate $\pi$ using the formula
   \begin{equation}
      \pi = 4 \int_0^1 \frac{1}{1 + x^2} \,\hbox{d} x
   \end{equation}

   We will illustrate the use of simple collective communications and timing
   routines.

   The program is called {\tt second.c} and it can be also found in
   {\tt pub/papers/meglicki/mpi/simplempi}.

   As in the first lesson, we shall list {\tt second.c} first. 
   Next, we shall present its synopsis, we'll compile it and run it on the LAM
   engine.
   Then we shall discuss it in detail. And, finally, we shall summarise 
   the new MPI issues covered in this section.

   \subsection{The listing of {\tt preamble.c}}
   \label{s:preamble}

\begin{verbatim}
    MPI_Init ( &argc, &argv );
    MPI_Comm_size ( MPI_COMM_WORLD, &pool_size ); 
    MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
    MPI_Get_processor_name ( my_node_name, &node_name_length );
    MPI_Attr_get ( MPI_COMM_WORLD, MPI_HOST, (void**) &host_rank_ptr,
                   &found_flag );
    if ( (int) host_rank_ptr != MPI_PROC_NULL ) host_rank = *host_rank_ptr; 
    else {
       MPI_Attr_get ( MPI_COMM_WORLD, MPI_IO, (void **) &host_rank_ptr,
                      &found_flag );
       if ( (int) host_rank_ptr != MPI_ANY_SOURCE ) host_rank = *host_rank_ptr;
       else {
          strcpy ( console_name, argv[1] );
          if ( 0 == strcmp ( my_node_name, console_name ) ) {
             host_rank = my_rank;
             i_am_the_master = B_TRUE;
          }
       }
    }
    if (i_am_the_master) {
       int_buffer[0] = host_rank;
       for (destination = 0; destination < pool_size; destination++)
          if (destination != my_rank) 
             MPI_Send(int_buffer, 1, MPI_INT, destination, 1001, 
                      MPI_COMM_WORLD);
    }
    else {
       MPI_Recv(int_buffer, 1, MPI_INT, MPI_ANY_SOURCE, 1001, 
                MPI_COMM_WORLD, &status);
       host_rank = int_buffer[0];
    }
    if (i_am_the_master) strcpy (char_buffer, my_node_name);
    MPI_Bcast(char_buffer, BUFSIZ, MPI_CHAR, host_rank, 
              MPI_COMM_WORLD);
    strcpy (host_name, char_buffer);
    if (! i_am_the_master) {
       sprintf (char_buffer, "Greetings to the master (%s, %d) from (%s, %d)",
                host_name, host_rank, my_node_name, my_rank);
       MPI_Send (char_buffer, strlen(char_buffer) + 1, MPI_CHAR, host_rank, 
                 2002, MPI_COMM_WORLD);
    }
    else {
       for (source = 0; source < pool_size; source++) {
          if (source != my_rank) {
             MPI_Recv(char_buffer, BUFSIZ, MPI_CHAR, source, 2002,
                      MPI_COMM_WORLD, &status);
             printf ("%s\n", char_buffer);
          }
       }
    }
\end{verbatim}

   \subsection{The listing of {\tt second.c}}
   \label{s:second}

\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include "mpi.h"
 
double f(a)
double a;
{
   return (4.0 / (1.0 + a*a));
}
 
int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 
 
#include "preamble.c"
 
    { 
       int n, i;
       double mypi, pi, h, sum, x, a;
 
       if (i_am_the_master) {
          printf("Enter the number of intervals: ");
          scanf("%d",&n);
          if (n==0) n=100;
       }
 
       MPI_Bcast(&n, 1, MPI_INT, host_rank, MPI_COMM_WORLD);
 
       h   = 1.0 / (double) n;
       sum = 0.0;
       for (i = my_rank + 1; i <= n; i += pool_size) {
          x = h * ((double)i - 0.5);
          sum += f(x);
       }
       mypi = h * sum;
 
       MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, host_rank,
                  MPI_COMM_WORLD);
 
       if (i_am_the_master) printf("pi is approximately %.16f\n", pi);
    }
 
    MPI_Finalize ();
}
\end{verbatim}

   \subsection{The synopsis of {\tt second.c}}

      The integral 
      \begin{equation}
         \int_0^1 \frac{4}{1 + x^2} \,\hbox{d}x
      \end{equation}
      is calculated by dividing the interval [0,1] into a number of
      subintervals and approximating the integrated function within each
      subinterval by a constant value which corresponds to the middle 
      of the subinterval. The master process obtains the number of subdivisions
      from standard input and broadcasts it to other processes.

      The processes perform summations each on its own portion of work. 
      At the end the partial results are added in another collective MPI
      operation called {\em reduction}.

      The answer is written on standard output by the master process,
      which, otherwise, participates equally in all other computations.

   \subsection{Compiling and running {\tt second.c} under LAM}

      To compile {\tt second.c} type:
\begin{verbatim}
gustav@cisr:~/mpi/simplempi 260 $ hcc -o second second.c -lmpi
gustav@cisr:~/mpi/simplempi 261 $ 
\end{verbatim}
      The file {\tt preamble.c} must exist in the same directory.

      In order to run the program, type the following in the same window from
      which the LAM machine was started. Enter a number of intervals
      when the program prompts for input:
\begin{verbatim}
gustav@cisr:~/mpi/simplempi 269 $ mpirun -w N `pwd`/second -- `uname -n`
Greetings to the master (cisr, 0) from (mercury, 1)
Greetings to the master (cisr, 0) from (venus, 2)
Greetings to the master (cisr, 0) from (earth, 3)
Greetings to the master (cisr, 0) from (mars, 4)
Greetings to the master (cisr, 0) from (jupiter, 5)
Greetings to the master (cisr, 0) from (saturn, 6)
Greetings to the master (cisr, 0) from (uranus, 7)
Greetings to the master (cisr, 0) from (bacchus, 8)
Greetings to the master (cisr, 0) from (ceres, 9)
Greetings to the master (cisr, 0) from (diana, 10)
Greetings to the master (cisr, 0) from (minerva, 11)
Greetings to the master (cisr, 0) from (vesta, 12)
Enter the number of intervals: 400
pi is approximately 3.1415931744231274
gustav@cisr:~/mpi/simplempi 270 $ 
\end{verbatim}

   \subsection{{\tt second.c} in detail -- data reductions}
   \label{s:reduction}


      We have already discussed \verb|MPI_Bcast| in section \ref{s:broadcast}.
      Thus the functioning of the program up to \verb|MPI_Bcast| should be
      clear. The computation itself is straightforward. Every process finds
      the length of the subinterval, initialises \verb|sum| to 0, and adds
      up the values of \verb|f(x)| which correspond to its subintervals.
      The only slighly confusing thing in this computation may be that instead
      of working on one simply connected subset of the domain, each process
      jumps over \verb|pool_size| of other subintervals. Thus the subintervals
      a given process works on are interleaved with subintervals belonging to
      other processes.

      At the end of this computation each process puts its own area under the
      curve of {\em f(x)\/} in \verb|mypi|. All processes then call 
\begin{verbatim}
MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, host_rank,
           MPI_COMM_WORLD);
\end{verbatim}
      
      This operation is discussed in ``Using MPI...'' on page 25, and
      in ``MPI: A Message-Passing Interface Standard'' on page 111, section
      \htmladdnormallink{4.9}
      {http://www.mcs.anl.gov/mpi/mpi-report/node76.html#Node76}. 
      
      All processes which issue this operation, pack their data into
      \verb|mypi| and send it out to a process whose rank is given by the 6th
      argument: \verb|host_rank| in our case. The type of data put in
      \verb|mypi| is given by the 4th argument and the number of data items
      in \verb|mypi| is given by the third argument.
      
      The process number
      \verb|host_rank| receives the data from all processes (including 
      itself), and performs a reduction operation on
      that data. The reduction operation is defined by the 5th argument. 
      Here it is \verb|MPI_SUM|, i.e., a summation. The result is placed
      in the second argument, \verb|pi|. Only processes which belong
      to the communicator \verb|MPI_COMM_WORLD|, specified by the last
      argument to \verb|MPI_Reduce|, participate in the transaction.

      MPI does not actually specify who is to perform the 
      final summation. On some supercomputers, there may be a special 
      circuitry within the network itself for doing such things.
      For example on the Connection Machine CM5, integer reductions are
      performed by the network, whereas floating point reductions are
      peformed by the destination node. CM5 network reductions are much faster
      than node reductions. But on the farms, it is almost certain that the
      reduction will be performed by the destination node.

      There is a large number of predefined reduction operations, for example
      \verb|MPI_SUM| used above, \verb|MPI_PROD|, \verb|MPI_MIN|,
      \verb|MPI_MAX|, etc. These are discussed on page 113 of ``MPI: A
      Message-Passing Interface Standard'', section 
      \htmladdnormallink{4.9.2}
      {http://www.mcs.anl.gov/mpi/mpi-report/node77.html#Node77}. 
      The programmer can
      define her own reduction operations using function 
      \verb|MPI_Op_create|, discussed in section 
      \htmladdnormallink{4.9.4}
      {http://www.mcs.anl.gov/mpi/mpi-report/node79.html#Node79}, 
      page 118.

   \subsection{{\tt third.c}, or {\tt second.c} with timing}

      We shall now modify our {\tt second.c} by adding the timing calls to the
      code. The modified code lives, like {\tt second.c}, in our ftp anonymous
      area on {\tt cisr.anu.edu.au}, directory 
      {\tt pub/papers/meglicki/mpi/simplempi}, file {\tt third.c}.
      The listing of {\tt third.c} is provided in the next subsection.

      \subsubsection{The listing of {\tt third.c}}
      \label{s:third}
\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include "mpi.h"

double f(a)
double a;
{
   return (4.0 / (1.0 + a*a));
}

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    { 
       int n, i;
       double mypi, pi, h, sum, x, a, start_time, end_time,
              communication_time, computation_time;

       while (1) {
          if (i_am_the_master) {
             printf("Enter the number of intervals: (0 quits) ");
             scanf("%d",&n);
          }

          start_time = MPI_Wtime();
          MPI_Bcast(&n, 1, MPI_INT, host_rank, MPI_COMM_WORLD);
          end_time = MPI_Wtime();
          communication_time = end_time - start_time;

          if (n == 0) break;
          else {

             start_time = MPI_Wtime();

             h   = 1.0 / (double) n;
             sum = 0.0;
             for (i = my_rank + 1; i <= n; i += pool_size) {
                x = h * ((double)i - 0.5);
                sum += f(x);
             }
             mypi = h * sum;

             end_time = MPI_Wtime();
             computation_time = end_time - start_time;
      
             start_time = MPI_Wtime();
             MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, host_rank,
                        MPI_COMM_WORLD);
             end_time = MPI_Wtime();
             communication_time = communication_time + end_time - start_time;
      
             if (i_am_the_master) {
                printf("pi is approximately %.16f\n", pi);
                printf("computation time %.16f\n", computation_time);
                printf("communication time %.16f\n", communication_time);
             }
          }
       }
    }

    MPI_Finalize ();
}
\end{verbatim}

      \subsubsection{Compilation and running of {\tt third.c}}

         In order to compile {\tt third.c} type
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 228 $ hcc -o third third.c -lmpi
gustav@jupiter:~/mpi/simplempi 229 $ 
\end{verbatim}
         To run it type:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 245 $ mpirun -w N `pwd`/third -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
Enter the number of intervals: (0 quits) 100000
pi is approximately 3.1415926535981265
computation time 0.0446860790252686
communication time 0.2531449794769287
Enter the number of intervals: (0 quits) 1000000
pi is approximately 3.1415926535898757
computation time 0.4425549507141113
communication time 1.1623179912567139
Enter the number of intervals: (0 quits) 10000000
pi is approximately 3.1415926535897714
computation time 4.1224819421768188
communication time 11.3426218032836914
Enter the number of intervals: (0 quits) 0
gustav@jupiter:~/mpi/simplempi 246 $ 
\end{verbatim}

         Observe that this time we have modified the logic of the program 
         somewhat, so that now we can repeat the
         computations without having to go through selection of the master
         process. The program also times both computations and communication
         operations.

      \subsubsection{About barrier synchronisation}
      \label{s:barrier}

         The first thing you should observe here is that the communication
         time increases as the number of computations increases. Why should it
         be so, if, clearly, in this case no more data is being transferred
         for n = 10,000,000 than for n = 10,000? 

         The answer to this question is that a reduction and broadcast
         operations are barrier operations. This means that if some processors
         are slower than others, as is the case on our farm, you will have to
         wait for the slow ones to finish their computations before the
         reduction operation can take place. Thus even though there is no more
         data transfer for n = 100,000 than for n = 10,000,000 in our program,
         there is more computation to be carried out, which results in slow
         processors lagging behind the fast ones more, 
         and this, in turn, increases
         the waiting time, whenever a barrier set by a broadcast or a reduction
         operation is encountered.

      \subsubsection{Timing MPI programs}
      \label{s:timing}

         The function \verb|MPI_Wtime| returns a floating point number of
         seconds, representing {\em elapsed\/} wall-clock time since some time
         in the past. Note that, unlike {\sc unix} function {\tt clock}, 
         \verb|MPI_Wtime| does not return the CPU time. The time returned by 
         \verb|MPI_Wtime| may be different for different nodes. Thus, the times
         returned by {\tt third.c} refer only to the master node.

         In the next lesson we will demonstrate how to {\em arm\/} a program
         and obtain information about the behaviour of other nodes as well.

   \subsection{Summary}

      In our second MPI lesson we met two new functions:

      \begin{itemize}
         \item \verb|MPI_Reduce| performs the reduction operation on 
               data sent to a designated process; explained in section
               \ref{s:reduction}
         \item \verb|MPI_Wtime| times the process with respect to wall clock;
               explained in section \ref{s:timing}
      \end{itemize}

      We also demonstrated how to distribute a computational task between
      processors, and how we may end up producing a code which does not
      scale with the size of the problem if the task is partitioned
      equally, regardless of the differences in the speed of the nodes
      participating in the cluster.

   \section{Lesson 3 -- the ``bank queue'' model}

   In this lesson we shall consolidate what we have learnt so far about MPI. We
   shall also cover a parallel computational paradigm, which is variously
   referred to as ``hungry puppies'', ``bank queue'', ``master-slaves'', ``task
   farm'', ``job queue'', and ``self-scheduling''. 
   My preferred name is ``bank queue'', because
   this is how most banks operate in Australia. But on my recent
   visit to Germany, for
   example, I have noticed that their banks don't use this practice, so
   the analogy may be lost on them.

   This paradigm suits especially
   heterogeneous farms, where we may have to utilise a number of different
   processors. As we have seen in the previous example, distributing tasks
   of equal size on such a farm may result in great delays on every barrier
   operation, and the execution of the program will be slowed down to the speed
   of the slowest processor. The ``bank queue'' paradigm overcomes this
   difficulty.

   \subsection{The ``bank queue'' paradigm}

   Now, how do banks operate in Australia? If you go, for example, to the
   Commonwealth Bank, there is just one queue to the tellers, regardless of
   the number of the tellers. Customers wait in
   the queue until a clerk in one of the tellers finishes with his or her 
   work and asks the next customer in the queue. This is
   an amazing organisational innovation, because slower clerks, or clerks who
   have more work to do will not slow down the processing of the queue.
   In other words, the queue itself is always processed as fast as the
   resources available to the bank branch allow.

   We can design a parallel computer program in a similar way. We may have
   a master process which organises the queue of tasks. Slave processes
   correspond to bank clerks working at the tellers. A slave process
   asks the master for a new task as soon as it has finished with the old one. 
   This way faster processes do not have to wait for the slower ones, 
   and the execution of the program proceeds as fast as possible, 
   given available resources.

   \subsection{How to obtain {\tt fourth.c}}

   Our example program is called {\tt fourth.c} and it is available
   by ftp anonymous from {\tt cisr.anu.edu.au}, directory
   {\tt /pub/papers/meglicki/mpi/simplempi}, file {\tt fourth.c}.
   I also provide a full listing of this program in the next section.
   The program is based on the matrix vector multiplication example
   in ``Using MPI...'', page 29, section 3.5, but our version has been
   written in C and is ``armed''. 
   
   Remember that you will also need {\tt preamble.c} (section \ref{s:preamble})
   in order to compile {\tt fourth.c}.

   \subsection{``Arming'' a program}

   ``Arming'' a program means inserting special diagnostic statements so that
   you can see what the program is doing during execution. 
   The simplest example of such a diagnostic statement is
\begin{verbatim}
#ifdef DEBUG
   fprintf(stderr, "...", ...);
#endif
\end{verbatim}
   MPI programs running under P4 based MPICH can be armed with special
   statements which will result in dumping an execution log. That log
   can be used later with a special program called {\tt upshot} in order
   to observe communications which took place 
   between processes during program execution.
   A similar facility called LAMLED is also available under LAM. 

   However, a great deal of information can be obtained using a simple
   {\tt printf} and we shall see in {\tt fourth.c} how this can be done
   in case of a parallel program which runs on a farm. 

   Like armed grenades, ``armed'' programs can sometimes explode spewing
   enormous amount of information, which is of limited interest
   or use to anyone but the author of the program. End user applications
   should not be ``armed''.

   \subsection{The listing of {\tt fourth.c}}
   \label{s:fourth}

\begin{verbatim}

#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include "mpi.h"

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    {

#define COLS 100
#define ROWS 100

        int a[ROWS][COLS], b[ROWS], c[ROWS], i, j;

        if (i_am_the_master) {

           int row, count, sender;

           for (i = 0; i < COLS; i++) {
              b[i] = 1;
              for (j = 0; j < ROWS; j++) a[i][j] = i;
           }

           MPI_Bcast(b, ROWS, MPI_INT, host_rank, MPI_COMM_WORLD);

           count = 0;
           for (destination = 0; destination < pool_size; destination++) {
              if (destination != my_rank) {
                 for (j = 0; j < COLS; j++) int_buffer[j] = a[count][j];
                 MPI_Send(int_buffer, COLS, MPI_INT, destination, count,
                          MPI_COMM_WORLD);
                 printf("sent row %d to %d\n", count, destination);
                 count = count + 1;
              }
           }

           for (i = 0; i < ROWS; i++) {
              MPI_Recv (int_buffer, BUFSIZ, MPI_INT, MPI_ANY_SOURCE, 
                        MPI_ANY_TAG, MPI_COMM_WORLD, &status);
              sender = status.MPI_SOURCE;
              row = status.MPI_TAG;
              c[row] = int_buffer[0];
              printf("\treceived row %d from %d\n", row, sender);
              if (count < ROWS) {
                 for (j = 0; j < COLS; j++) int_buffer[j] = a[count][j];
                 MPI_Send(int_buffer, COLS, MPI_INT, sender, count,
                          MPI_COMM_WORLD);
                 printf("sent row %d to %d\n", count, sender);
                 count = count + 1;
              }
              else {
                 MPI_Send(0, 0, MPI_INT, sender, ROWS, MPI_COMM_WORLD);
                 printf("terminated process %d with tag %d\n", sender, ROWS);

                 /*
                  Exercise: why will the master process hang if we use
                  0 instead of ROWS to terminate the slave, as suggested 
                  in the MPI book?
                  */
              }
           }
        }
        else { /* I am not the master */

           int sum, row;
           FILE *log_file;

           log_file = fopen ("/tmp/gustav_log", "w");

           MPI_Bcast(b, COLS, MPI_INT, host_rank, MPI_COMM_WORLD);
           fprintf(log_file, "received broadcast from %d\n", host_rank);
           fflush(log_file);
           MPI_Recv(int_buffer, COLS, MPI_INT, host_rank, MPI_ANY_TAG,
                    MPI_COMM_WORLD, &status);
           fprintf(log_file, "received a message from %d, tag %d\n",
                   status.MPI_SOURCE, status.MPI_TAG);
           fflush(log_file);
           while (status.MPI_TAG != ROWS) { /* The job is not finished */
              row = status.MPI_TAG; sum = 0;
              for (i = 0; i < COLS; i++) sum = sum + int_buffer[i] * b[i];
              int_buffer[0] = sum;
              MPI_Send (int_buffer, 1, MPI_INT, host_rank, row, MPI_COMM_WORLD);
              fprintf(log_file, "sent row %d to %d\n", row, host_rank);
              fflush(log_file);
              MPI_Recv (int_buffer, COLS, MPI_INT, host_rank, MPI_ANY_TAG,
                        MPI_COMM_WORLD, &status);
              fprintf(log_file, "received a message from %d, tag %d\n",
                      status.MPI_SOURCE, status.MPI_TAG);
              fflush(log_file);
           }
           fprintf(log_file, "exiting on  tag %d\n", status.MPI_TAG);
           fflush(log_file);
        }
    }

    MPI_Finalize ();
}
\end{verbatim}

   \subsection{The synopsis of {\tt fourth.c}}

      This program performs a multiplication of a 100 x 100 matrix by a vector
      using the ``bank queue'' paradigm. 
      
      First the master process broadcasts
      the content of the vector to all processes. Then the first
      \verb|pool_size - 1| rows of the matrix are sent to the workers.
      Each worker multiplies its row of the matrix by the vector and
      returns the result to the master. Upon having received the result
      the master sends the next row of the matrix to the worker until
      the whole matrix has already been handed out. In the latter case,
      the master sends a termination signal to the worker.

      All communication operations are logged by the master on standard output.
      The workers log the communication operations on special log files which
      they open in {\tt /tmp}. 

   \subsection{Modifying and compiling {\tt fourth.c}}

      Before you compile {\tt fourth.c} you should edit the file and change the
      name of the log file, which the worker processes will open in their 
      ``/tmp'' directories, to something more suitable for your system. 
      In this program the
      name of the log file is hardwired -- not an elegant procedure 
      in general, but simple enough to fit on one line of code.

      Search for the line:
\begin{verbatim}
log_file = fopen ("/tmp/gustav_log", "w");
\end{verbatim}
      and change \verb|"/tmp/gustav_log"|, for example to 
      \verb|"/tmp/<whatever your login name is>_log"|. The name of the file
      should be unique to your {\tt uid}, in case more than one user
      works with this tutorial on the same system at the same time. You must
      also write that file in \verb|/tmp|, so that processes running on
      different machines would not overwrite one another's log files. 
      In particular, do not
      open the log file in your home directory, if it is mounted on other
      machines participating in the farm under the same name.

      After you have modified {\tt fourth.c}, compile it as follows:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 241 $ hcc -o fourth fourth.c -lmpi
gustav@jupiter:~/mpi/simplempi 242 $ 
\end{verbatim}

      Running {\tt fourth.c} will produce a copious amount of output. For this
      reason we shall discuss it in a separate section.

   \subsection{Running {\tt fourth.c} under LAM}

      To run {\tt fourth.c} under LAM type
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 203 $ mpirun -w N `pwd`/fourth -- `uname -n`
\end{verbatim}
      in the same window, in which you have run \verb|lamboot -v|.

      If everything works, you should see an output similar to the one shown
      in the next section.

      \subsubsection{The output of {\tt fourth.c}}
\begin{verbatim}
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
sent row 0 to 0
sent row 1 to 1
sent row 2 to 2
sent row 3 to 3
sent row 4 to 4
sent row 5 to 6
sent row 6 to 7
sent row 7 to 8
sent row 8 to 9
sent row 9 to 10
sent row 10 to 11
sent row 11 to 12
        received row 0 from 0
sent row 12 to 0
        received row 1 from 1
sent row 13 to 1
        received row 2 from 2
sent row 14 to 2

   ...

        received row 85 from 9
sent row 97 to 9
        received row 86 from 2
sent row 98 to 2
        received row 87 from 3
sent row 99 to 3
        received row 88 from 4
terminated process 4 with tag 100
        received row 89 from 6
terminated process 6 with tag 100
        received row 90 from 7
terminated process 7 with tag 100
        received row 91 from 8
terminated process 8 with tag 100
        received row 92 from 10
terminated process 10 with tag 100
        received row 93 from 11
terminated process 11 with tag 100
        received row 94 from 12
terminated process 12 with tag 100
        received row 95 from 0
terminated process 0 with tag 100
        received row 96 from 1
terminated process 1 with tag 100
        received row 98 from 2
terminated process 2 with tag 100
        received row 97 from 9
terminated process 9 with tag 100
        received row 99 from 3
terminated process 3 with tag 100
\end{verbatim}

      \subsubsection{The discussion}

      This output can be divided in 4 parts. The first part covers 
      the workers paying homage to the master. That's our usual output
      from all programs. 

      Then we can see the master process sending the first 11
      rows to worker processes. The master does not send a job to itself
      (typical, isn't it?).

      The third distinct part is the longest. Here we see the master process
      receiving answers from worker processes and immediately
      sending the next row to each process that comes back with the answer.

      Finally, when the master will have sent row 99, it begins to terminate
      processes returning with answers by sending them empty messages
      with tag value of 100. 

      \subsubsection{Workers' logs}

      In order to see workers' logs, type something similar to:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 209 $ rsh diana cat /tmp/gustav_log
received broadcast from 5
received a message from 5, tag 9
sent row 9 to 5
received a message from 5, tag 20
sent row 20 to 5
received a message from 5, tag 32
sent row 32 to 5
received a message from 5, tag 44
sent row 44 to 5
received a message from 5, tag 56
sent row 56 to 5
received a message from 5, tag 68
sent row 68 to 5
received a message from 5, tag 80
sent row 80 to 5
received a message from 5, tag 92
sent row 92 to 5
received a message from 5, tag 100
exiting on  tag 100
gustav@jupiter:~/mpi/simplempi 210 $ 
\end{verbatim}

      You can see that the worker process
      is informed about the matrix row number, which it has received
      in the message, through the message tag value. When the answer is sent
      back to the master process the same tag value is used, so that the master
      can put the answer in an appropriate slot.

      When the worker process receives a message with tag value 100 (the
      highest row number of the matrix is 99), it interprets it 
      as a termination signal and exits.

   \subsection{{\tt fourth.c} in detail}
   \label{s:fourth-detail}

      This program has a typical master-worker structure:
\begin{verbatim}
   if (i_am_the_master) {
      <initialise the problem>;
      <initialise workers>;
      <send initial batch of the tasks to the workers>;
      REPEAT {
         <receive an answer from a worker>;
         <send another task to the worker if the job is not finished>;
            or
         <send a termination signal to the worker if the job is done>;
      } UNTIL <all workers are dead>;
   } 
   else {
      <receive initialisation from the master>;
      REPEAT {
         <receive a task from the master>;
         <work on the task>;
         <send the answer to the master>;
      } UNTIL <terminated by the master>;
   }
\end{verbatim}

      \subsubsection{The master part}
      \label{s:fourth-master-part}

      The initialisation is trivial.

      Once the vector \verb|b| and the matrix \verb|a[i][j]| are initialised,
      the master broadcasts the vector to all workers using
\begin{verbatim}
MPI_Bcast(b, ROWS, MPI_INT, host_rank, MPI_COMM_WORLD);
\end{verbatim}
      This ``initialises'' the workers.

      Then the master process sends a row of the matrix to each of the workers.
      The number of the row is conveyed as the message tag:
\begin{verbatim}
for (j = 0; j < COLS; j++) int_buffer[j] = a[count][j];
MPI_Send(int_buffer, COLS, MPI_INT, destination, count,
         MPI_COMM_WORLD);
\end{verbatim}
      In this manner, the master sends the initial batch of the tasks to
      all workers.

      Now we have the \verb|for| loop, which collects replies from all 
      workers for all rows, both the ones that have been sent already 
      and the ones that will be sent within this loop.

      The master process waits for any message with any tag from any source.
      The rank of the sender process and the row of the matrix the answer
      refers to are obtained from structure \verb|status|:
\begin{verbatim}
   MPI_Recv (int_buffer, BUFSIZ, MPI_INT, MPI_ANY_SOURCE,
             MPI_ANY_TAG, MPI_COMM_WORLD, &status);
   sender = status.MPI_SOURCE;
   row = status.MPI_TAG;
\end{verbatim}
      This is a new aspect of \verb|MPI_Recv| which we haven't covered yet.
      The structure \verb|MPI_Status| is discussed in section 
      \htmladdnormallink{3.2.5}
      {http://www.mcs.anl.gov/mpi/mpi-report/node34.html#Node34}, 
      page 20
      of ``MPI: A Message-Passing Interface Standard''. This structure
      contains two mandatory fields named \verb|MPI_SOURCE| and 
      \verb|MPI_TAG|, and {\em may\/} also contain additional fields. The
      length of the received message can be also extracted from
      \verb|MPI_Status| using function \verb|MPI_Get_count|. The latter
      is most useful in combination with another function
      \verb|MPI_Probe|, which can be used to obtain some information
      about the message, and allocating appropriate amount of space
      for the message, before the message is actually read into the buffer.
      \verb|MPI_Probe| is discussed in section
      \htmladdnormallink{3.8}
      {http://www.mcs.anl.gov/mpi/mpi-report/node49.html#Node49} 
      of MPI Standard.

      After having read and processed the answer, the master checks if
      the whole matrix has been finished
\begin{verbatim}
if (count < ROWS) {...
\end{verbatim}
      and if there is still some work to be done, another row of the
      matrix is sent to the same worker:
\begin{verbatim}
for (j = 0; j < COLS; j++) int_buffer[j] = a[count][j];
MPI_Send(int_buffer, COLS, MPI_INT, sender, count,
         MPI_COMM_WORLD);
\end{verbatim}
      Otherwise, i.e., if the matrix has already been finished, the master
      sends to the worker a null message with tag \verb|ROWS|. Because
      in C the matrix will never have a row number \verb|ROWS|, and
      because all workers know about the dimensions of the matrix, receiving
      a message with that tag is a signal to the worker that it's time
      to exit:
\begin{verbatim}
   MPI_Send(0, 0, MPI_INT, sender, ROWS, MPI_COMM_WORLD);
\end{verbatim}

      In ``Using MPI...'' you will find that the termination signal in
      their example is 0. But if we did so here, we would terminate process
      rank 0 the moment we send the zeroth row of the matrix to it. The 
      master process
      would then hang, because process rank 0 would not have sent the
      result of multiplication of the zeroth row of the matrix by vector b
      back to the master. 
      
      There are
      two important differences between the example from ``Using MPI...'' and
      our example. First, their example is in Fortran, and in Fortran, unlike
      in C, the first row of a matrix has index 1, not 0. The second difference
      is that in ``Using MPI...'' the master process always has a rank 0,
      whereas in our case it can have any rank, depending on the machine you
      start your LAM engine from.

      \subsubsection{The worker part}

      The worker begins by receiving the broadcast from the master:
\begin{verbatim}
   MPI_Bcast(b, COLS, MPI_INT, host_rank, MPI_COMM_WORLD);
\end{verbatim}
      Then all workers receive their first batch to work on:
\begin{verbatim}
   MPI_Recv(int_buffer, COLS, MPI_INT, host_rank, MPI_ANY_TAG,
            MPI_COMM_WORLD, &status);
\end{verbatim}
      At this stage every worker enters a loop at the top of
      which it checks if the tag of the received message is
      equal to \verb|ROWS|. If it is, the worker exits.

      If the tag of the message is different from \verb|ROWS|,
      the worker performs the computation:
\begin{verbatim}
   sum = 0;
   for (i = 0; i < COLS; i++) sum = sum + int_buffer[i] * b[i];
\end{verbatim}
      and sends the answer back to the master using the row number 
      as the tag:
\begin{verbatim}
   int_buffer[0] = sum;
   MPI_Send (int_buffer, 1, MPI_INT, host_rank, row, MPI_COMM_WORLD);
\end{verbatim}
      At this stage every worker knows that it should receive either
      the next job or the termination tag from the master:
\begin{verbatim}
   MPI_Recv (int_buffer, COLS, MPI_INT, host_rank, MPI_ANY_TAG,
             MPI_COMM_WORLD, &status);
\end{verbatim}
      When this message is received, the worker goes back to the top
      of the loop, where the tag is inspected.

   \subsection{Summary}

      We have covered one new element of MPI in this section,
      namely the structure of the receive status, and the ways
      of accessing information contained in it. This issue
      is discussed in section \ref{s:fourth-master-part}.

   \section{Lesson 4 -- working with communicators}

      In this lesson we shall talk about manipulating communicators. Our
      program will calculate $\pi$ using a Monte Carlo method. We will have one
      separate process serving random numbers to all other processes. That
      process will be excluded from a new communicator, which will be used by
      remaining processes for reduction operations. The master process, i.e.,
      the process responsible for writing to the VDU, must not be the same
      process as the process that serves random numbers. So in this program we
      have two special processes: the random number server and 
      the process which writes to the VDU.

      Our program is called {\tt fifth.c}. You can obtain it by 
      ftp anonymous from {\tt cisr.anu.edu.au}, directory
      {\tt /pub/papers/meglicki/mpi/simplempi}, file {\tt fifth.c}.

      You must also have a file {\tt preamble.c}, which is listed in
      section \ref{s:preamble}. That file can be also found in the
      same directory on our ftp server.

      This program is also discussed in section 3.7 of ``Using MPI...''.
      {\tt fifth.c} is not identical with the Argonne program though.
      I had to make a number of small changes in order to get it going
      under LAM on our farm.

      The next section presents the listing of {\tt fifth.c}.

      \subsection{The listing of {\tt fifth.c}}
      \label{s:fifth}

\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <math.h>
#include <limits.h>
#include "mpi.h"

#define CHUNKSIZE 1000
#define REQUEST 1
#define REPLY 2

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    {
       MPI_Comm world, workers;
       MPI_Group world_group, worker_group;
       int ranks[1], server, request, rands[CHUNKSIZE], done, in, out,
           max, my_worker_rank, iter, totalin, totalout, i;
       double x, y, epsilon, error, Pi;
       boolean_t i_am_the_server = B_FALSE;

       if (i_am_the_master) {
          printf ( "epsilon = " );
          scanf ( "%lf", &epsilon );
       }

       MPI_Bcast ( &epsilon, 1, MPI_DOUBLE, host_rank, MPI_COMM_WORLD );

       world = MPI_COMM_WORLD;

       /* Server must not be the same process as the host. We must have
          sufficient number of processes so as not to hit a negative
          rank number here! */

       server = (host_rank == pool_size - 1)? pool_size - 2 : pool_size - 1;
       if (my_rank == server) i_am_the_server = B_TRUE;

       MPI_Comm_group ( MPI_COMM_WORLD, &world_group );
       ranks[0] = server;
       MPI_Group_excl ( world_group, 1, ranks, &worker_group );
       MPI_Comm_create ( world, worker_group, &workers );
       MPI_Group_free ( &worker_group );
       MPI_Group_free ( &world_group );

       if (i_am_the_server) {
          do {
             MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST, world,
                      &status);
             if (request) {
                for (i = 0; i < CHUNKSIZE; i++) rands[i] = random();
                MPI_Send(rands, CHUNKSIZE, MPI_INT, status.MPI_SOURCE,
                         REPLY, world);
             }
          } while ( request > 0 );
       }
       else {
          request = 1;
          done = in = out = 0;
          max = INT_MAX;
          MPI_Send ( &request, 1, MPI_INT, server, REQUEST, world );
          MPI_Comm_rank ( workers, &my_worker_rank );
          iter = 0;
          while (! done) {
             iter++;
             request = 1;
             MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY, world,
                       &status );
             for (i = 0; i < CHUNKSIZE;) {
                x = (((double) rands[i++])/max) * 2 - 1;
                y = (((double) rands[i++])/max) * 2 - 1;
                if (x * x + y * y < 1.0) in++;
                else out++;
             }
             MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM, workers);
             MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM, workers);
             Pi = (4.0*totalin)/(totalin + totalout);
             error = fabs( Pi - 3.141592653589793238462643);
             done = ((error < epsilon) || ((totalin + totalout) > 1000000));
             request = (done) ? 0 : 1;
             if (i_am_the_master) {
                printf("\rpi = %23.20lf", Pi); fflush(stdout);
                MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
             }
             else {
                if (request)
                   MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
             }
          }
          if (i_am_the_master) {
             printf("\npoints: %d\nin: %d, out: %d\n", totalin+totalout, 
                    totalin, totalout);
          }
          MPI_Comm_free(&workers);
       }
    }

    MPI_Finalize ();
}
\end{verbatim}

   \subsection{The synopsis of {\tt fifth.c}}

      $\pi$ is computed by generating random points within a square, whose
      dimensions are 2 x 2, and checking how many of these 
      points end up within a circle of radius 1. The ratio of that
      number to the total number of points generated by the random number
      server should be close to $\pi/4$ for a sufficiently large number of
      points.

      The program begins by asking the user for required accuracy. 
      The computation of $\pi$ will go on until either the required
      accuracy is reached, or until the number of points generated
      by the random number server exceeds certain hardwired limit.

      The process responsible for writing to the VDU updates the displayed
      value of $\pi$ every time the results are collected from
      the workers. 

   \subsection{Compiling and running {\tt fifth.c} under LAM}

      In order to compile {\tt fifth.c} type
\begin{verbatim}
<mpi 235 $ hcc -o fifth fifth.c -lmpi /usr/ucblib/libucb.a
gustav@jupiter:~/mpi/simplempi 236 $
\end{verbatim}

      Observe that my bash command line scrolled to the left, as the input line
      became longer. This time we also link the program with 
      {\tt /usr/ucblib/libucb.a} in order to attach the BSD function
      {\tt random()}. This will work under Solaris 2.3, but you may have
      to modify this compilation procedure if you have a different system.

      To run the program type:
\begin{verbatim}
gustav@jupiter:~/mpi/simplempi 238 $ mpirun -w N `pwd`/fifth -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
epsilon = 0.0001
pi =  3.14161904761904775000
points: 126000
in: 98961, out: 27039
gustav@jupiter:~/mpi/simplempi 239 $ 
\end{verbatim}

   \subsection{{\tt fifth.c} in detail}

   After obtaining the value of $\epsilon$ from the user, the master broadcasts
   it to other processes. Afterwards we select a random number server. This
   process, as I have remarked before, should not be the same as the master
   process:
\begin{verbatim}
server = (host_rank == pool_size - 1)? pool_size - 2 : pool_size - 1;
if (my_rank == server) i_am_the_server = B_TRUE;
\end{verbatim}

   Once the server process is selected, we have to create a new communicator,
   which does not contain that process.

   \subsubsection{Extraction, manipulation, and creation of process groups
                  and communicators}
   \label{s:groups}

   We begin the construction of the new communicator by
   extracting a process group, associated with the
   \verb|MPI_COMM_WORLD| communicator and call it \verb|world_group|:
\begin{verbatim}
MPI_Comm_group ( MPI_COMM_WORLD, &world_group );
\end{verbatim}

   The next operation creates a new group, which will be called
   \verb|worker_group|, by {\em excluding\/} processes declared
   in the array \verb|ranks| from the \verb|world_group|:
\begin{verbatim}
ranks[0] = server;
MPI_Group_excl ( world_group, 1, ranks, &worker_group );
\end{verbatim}

   These two new functions, \verb|MPI_Comm_group| and 
   \verb|MPI_Group_excl|, are described in 
   \htmladdnormallink{Chapter 5}
   {http://www.mcs.anl.gov/mpi/mpi-report/node86.html#Node86}
   ``Groups, Contexts, and
   Communicators'', of ``MPI: A Message-Passing Interface Standard''. 

   The number of processes to be excluded from the new group
   must be passed as the second parameter
   to function \verb|MPI_Group_excl|. In our case we exclude only one
   process, the \verb|server|.

   Having created the new group, we can now create the communicator 
   associated with that group. We shall call that communicator \verb|workers|:
\begin{verbatim}
MPI_Comm_create ( world, worker_group, &workers );
\end{verbatim}

   Once the communicator has been created we can release 
   {\em the names\/} of groups \verb|worker_group| and \verb|world_group|:
\begin{verbatim}
MPI_Group_free ( &worker_group );
MPI_Group_free ( &world_group );
\end{verbatim}
   It is important to realise that this operation only frees {\em the names\/}
   of the groups, not the groups themselves.
   The groups survive inside the communicators, and if needed, they can 
   be retrieved again by calling \verb|MPI_Comm_group|.

   Process groups and communicators are a very important part of MPI. They
   make it possible to write parallel libraries, which can be accessed by user
   programs without colliding with their own nomenclature, and their own
   messages. The idea of process groups also lies at the basis of fault
   tolerant and virtually synchronous programming. But MPI, which is designed
   specifically for scientific and engineering computations, doesn't take it
   that far. 

      \subsubsection{The server process}

      The job of the server process is simple: receive a request from any
      worker within the \verb|world = MPI_COMM_WORLD| communicator, and 
      return an array of \verb|CHUNKSIZE| random numbers to that process if
      the \verb|request > 0|. Otherwise, exit:
\begin{verbatim}
do {
   MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST, world,
            &status);
   if (request) {
      for (i = 0; i < CHUNKSIZE; i++) rands[i] = random();
      MPI_Send(rands, CHUNKSIZE, MPI_INT, status.MPI_SOURCE,
               REPLY, world);
   }
} while ( request > 0 );
\end{verbatim}

      \subsubsection{The worker process}
      \label{s:allreduce}

      Worker processes have a harder life. They begin by sending a request
      for a \verb|CHUNKSIZE| of random numbers to the random number server:
\begin{verbatim}
MPI_Send ( &request, 1, MPI_INT, server, REQUEST, world );
\end{verbatim}

      Then they find their new rank number within the new \verb|workers|
      communicator:
\begin{verbatim}
MPI_Comm_rank ( workers, &my_worker_rank );
\end{verbatim}
      This number does not have to be the same as their rank number within
      the \verb|MPI_COMM_WORLD| communicator. We don't make any use of 
      that number in this program. You can print it on a log file
      if you would like to see it. Here we just want to make
      the point that a process can have several rank numbers associated
      with different communicators.

      Now the workers enter the \verb|while| loop. They receive a batch
      of random numbers from the server:
\begin{verbatim}
MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY, world,
          &status );
\end{verbatim}
      scale them and convert them into points which fit 
      into the 2 x 2 square. Then for every point they check if it 
      belongs to the
      circle, or if it is left outside. Depending on the result they
      increment either the \verb|in| or the \verb|out| counters.
\begin{verbatim}
for (i = 0; i < CHUNKSIZE;) {
   x = (((double) rands[i++])/max) * 2 - 1;
   y = (((double) rands[i++])/max) * 2 - 1;
   if (x * x + y * y < 1.0) in++;
   else out++;
}
\end{verbatim}

      Now they exchange their information with all other processes within
      the \verb|workers| communicator. They use a new MPI function, which
      we haven't seen yet:
\begin{verbatim}
MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM, workers);
MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM, workers);
\end{verbatim}
      This function works the same way as \verb|MPI_Reduce|, but this time
      the reduced data ends up {\em not\/} with one selected process, but
      with {\em all\/} processes participating in the communicator.
      Function \verb|MPI_Allreduce| is discussed in section 
      \htmladdnormallink{4.9.5}
      {http://www.mcs.anl.gov/mpi/mpi-report/node81.html#Node81}, 
      page 122, ``MPI: A Message-Passing Interface Standard''.

      The data to be contributed to the pool
      is placed in the first slot of the function, and the number
      of items is given by its third argument. The reduced data is written
      on the second argument. The MPI type (not the C type!) of the data
      is given by the fourth argument, and the type of the operation
      to be performed by the data is given by the fifth argument.

      Now {\em every\/} process can evaluate $\pi$, the accuracy of the
      result, and decide if it is time to quit:
\begin{verbatim}
Pi = (4.0*totalin)/(totalin + totalout);
error = fabs( Pi - 3.141592653589793238462643);
done = ((error < epsilon) || ((totalin + totalout) > 1000000));
\end{verbatim}
      They all should arrive at identical results, because they
      work with identical data, thanks to \verb|MPI_Allreduce|.

      If the job is done, they set the request to be sent to
      the random numbers server to 0, otherwise the request is set to 1:
\begin{verbatim}
request = (done) ? 0 : 1;
\end{verbatim}

      Observe that only the master process will send the termination
      signal to the random numbers server:
\begin{verbatim}
if (i_am_the_master) {
   printf("\rpi = %23.20lf", Pi); fflush(stdout);
   MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
}
else {
   if (request)
      MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
}
\end{verbatim}

      Before we quit the program, we release the \verb|workers| 
      communicator:
\begin{verbatim}
MPI_Comm_free(&workers);
\end{verbatim}

   \subsection{Summary}

   In this lesson we have encountered several new operations: group and
   communicator manipulation functions, and the all-reduce function:
   \begin{itemize}
      \item \verb|MPI_Comm_group| extracts a group associated with a
         communicator; discussed in section \ref{s:groups}
      \item \verb|MPI_Group_excl| creates a new group by exluding named
         processes from a given group; discussed in section \ref{s:groups}
      \item \verb|MPI_Comm_create| creates a new communicator associated
         with a given group; discussed in section \ref{s:groups}
      \item \verb|MPI_Group_free| releases the name associated with
         a given group; the group itself may survive inside a communicator
         if such was associated with it before; discussed in section
         \ref{s:groups}
      \item \verb|MPI_Comm_free| releases a communicator; discussed in section
         \ref{s:allreduce}
      \item \verb|MPI_Allreduce| performs a reduction operation on data
         contributed by participating processes and delivers the result to each
         process; discussed in section \ref{s:allreduce}
   \end{itemize}

   \section{Lesson 5 -- Cartesian grids}

      This is going to be a long lesson. It covers much of the material
      presented in Chapter 4, ``Intermediate MPI'', of ``Using MPI...''. The
      example program of Chapter 4 is in Fortran 77. Our example program is a
      little different, written in C, but it illustrates the same point. 

      We shall explore in this lesson the issue of process topologies
      and, in particular, the Cartesian topology. The Cartesian topology of
      process space is something you may wish to use for data parallel
      problems, for example solving partial differential equations
      on rectangular grids. 

      The data parallel approach usually is unsuitable for heterogeneous farms
      with widely differing CPU performance and load, because of the need
      to wait at the barrier for the slowest process
      whenever in-step operations are carried out.
      We have already encountered such a situation in 
      section \ref{s:barrier}. For such systems the ``bank queue'' model
      or something similarly flexible seems to be the best approach. But if you
      happen to have an Alpha farm, or an SP2, or just a bunch of unused 
      Jedis on a fast switch collecting
      dust somewhere in your garden shed, and if you 
      happen to have the whole farm to yourself so that
      all machines are equally loaded, you may seriously
      consider using such a system for a data parallel problem. And 
      being able to organise
      the processes into a Cartesian grid and shift the data
      between adjacent nodes may then come handy.

      This is exactly what we are going to illustrate in this lesson.

      The program, we'll be working with is called {\tt eighth.c} and you
      will find it in our ftp anonymous area on {\tt cisr.anu.edu.au},
      directory {\tt /pub/papers/meglicki/mpi/intermediate}, file
      {\tt eighth.c}. Remember to transfer the file {\tt preamble.c},
      section \ref{s:preamble}, from the same directory too. 

      If you are curious what happened to {\tt sixth.c}, and {\tt seventh.c},
      they can be also found in the same directory. They show {\tt eighth.c} in
      the process of development. 

      For your convenience the listing of
      the program is provided in the next section.

      \subsection{The listing of {\tt eighth.c}}
      \label{s:eighth}

\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include "mpi.h"

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    {

#define UPDOWN              0
#define SIDEWAYS            1
#define RIGHT               1
#define UP                  1

#define PROCESS_DIMENSIONS  2

#define PROCESS_ROWS        4
#define ROWS                4
#define DISPLAY_ROWS       16      /* must be PROCESS_ROWS * ROWS */

#define PROCESS_COLUMNS     3
#define COLUMNS             4
#define DISPLAY_COLUMNS    12      /* must be PROCESS_COLUMNS * COLUMNS */

        int divisions[PROCESS_DIMENSIONS] = {PROCESS_ROWS, PROCESS_COLUMNS};
        int periods[PROCESS_DIMENSIONS] = {0, 0};
        int reorder = 1;
        MPI_Comm cartesian_communicator;
        int my_cartesian_rank, my_coordinates[PROCESS_DIMENSIONS];
        int left_neighbour, right_neighbour, bottom_neighbour,
            top_neighbour;

        MPI_Cart_create ( MPI_COMM_WORLD, PROCESS_DIMENSIONS, divisions, 
                          periods, reorder, &cartesian_communicator );

        if (cartesian_communicator != MPI_COMM_NULL) {

           int matrix [ROWS][COLUMNS];
           int i, j;
           MPI_Datatype column_type;

           MPI_Comm_rank ( cartesian_communicator, &my_cartesian_rank );
           MPI_Cart_coords ( cartesian_communicator, my_cartesian_rank,
                             PROCESS_DIMENSIONS, my_coordinates );
           MPI_Cart_shift ( cartesian_communicator, SIDEWAYS, RIGHT, 
                            &left_neighbour, &right_neighbour );
           MPI_Cart_shift (cartesian_communicator, UPDOWN, UP,
                           &bottom_neighbour, &top_neighbour );

           if (! i_am_the_master ) {
              sprintf(char_buffer, "process %2d, cartesian %2d, \
coords (%2d,%2d), left %2d, right %2d, top %2d, bottom %2d",
                      my_rank, my_cartesian_rank, my_coordinates[0],
                      my_coordinates[1], left_neighbour, right_neighbour,
                      top_neighbour, bottom_neighbour);
              MPI_Send(char_buffer, strlen(char_buffer) + 1, MPI_CHAR,
                      host_rank, 3003, MPI_COMM_WORLD);
           } 
           else {

              int number_of_c_procs, count;

              number_of_c_procs = divisions[0] * divisions[1]; 
              for (count = 0; count < number_of_c_procs - 1; count++) {
                 MPI_Recv(char_buffer, BUFSIZ, MPI_CHAR, MPI_ANY_SOURCE, 3003,
                          MPI_COMM_WORLD, &status);
                 printf ("%s\n", char_buffer);
              }
              printf( "process %2d, cartesian %2d, \
coords (%2d,%2d), left %2d, right %2d, top %2d, bottom %2d\n",
                      my_rank, my_cartesian_rank, my_coordinates[0],
                      my_coordinates[1], left_neighbour, right_neighbour,
                      top_neighbour, bottom_neighbour);
           }

           for ( i = 0; i < ROWS; i++ ) {
              for ( j = 0; j < COLUMNS; j++ ) {
                 matrix [i][j] = my_cartesian_rank;
              }
           }

           if (! i_am_the_master ) 
              MPI_Send ( matrix, COLUMNS * ROWS, MPI_INT, host_rank, 3003,
                         cartesian_communicator );
           else
              collect_matrices ( cartesian_communicator, my_cartesian_rank,
                                 matrix, 3003 );

           MPI_Sendrecv ( &matrix[ROWS - 2][0], COLUMNS, MPI_INT, 
                          top_neighbour, 4004,
                          &matrix[0][0], COLUMNS, MPI_INT, bottom_neighbour, 
                          4004, 
                          cartesian_communicator, &status );

           MPI_Sendrecv ( &matrix[1][0], COLUMNS, MPI_INT, bottom_neighbour,
                          5005,
                          &matrix[ROWS - 1][0], COLUMNS, MPI_INT, 
                          top_neighbour, 5005,
                          cartesian_communicator, &status );

           if (! i_am_the_master ) 
              MPI_Send ( matrix, COLUMNS * ROWS, MPI_INT, host_rank, 6006,
                         cartesian_communicator );
           else 
              collect_matrices ( cartesian_communicator, my_cartesian_rank,
                                 matrix, 6006 );

           MPI_Type_vector (ROWS, 1, COLUMNS, MPI_INT, &column_type);
           MPI_Type_commit (&column_type);

           MPI_Sendrecv ( &matrix[0][1], 1, column_type, left_neighbour, 7007,
                          &matrix[0][COLUMNS - 1], 1, column_type, 
                          right_neighbour, 7007,
                          cartesian_communicator, &status );
 
           MPI_Sendrecv ( &matrix[0][COLUMNS - 2], 1, column_type, 
                          right_neighbour, 8008,
                          &matrix[0][0], 1, column_type, left_neighbour, 8008,
                          cartesian_communicator, &status );

           if (! i_am_the_master )
              MPI_Send ( matrix, COLUMNS * ROWS, MPI_INT, host_rank, 9009,
                         cartesian_communicator );
           else
              collect_matrices ( cartesian_communicator, my_cartesian_rank,
                                 matrix, 9009 );
        }
    }

    MPI_Finalize ();
}

int print_array (int array [DISPLAY_ROWS] [DISPLAY_COLUMNS], 
                          int vertical_break,
                          int horizontal_break)
{
   int k, l;

   printf ("\n");
   for (k = DISPLAY_ROWS - 1; k >= 0; k -- ) {
      for (l = 0; l < DISPLAY_COLUMNS; l ++ ) {
         if (l % horizontal_break == 0) printf (" ");
         printf ( "%2d ", array [k][l] );
      }
      printf ( "\n" );
      if (k % vertical_break == 0) printf ( "\n" );
   }
}

int collect_matrices (MPI_Comm cartesian_communicator, 
                      int my_cartesian_rank,
                      int matrix[ROWS][COLUMNS], 
                      int tag)
{
   int coordinates[PROCESS_DIMENSIONS];
   int client_matrix[ROWS][COLUMNS];
   int display[DISPLAY_ROWS][DISPLAY_COLUMNS];
   int i, j, k, l, source;
   MPI_Status status;
 
   for ( i = PROCESS_ROWS - 1; i >= 0; i -- ) {
      for ( j = 0; j < PROCESS_COLUMNS; j ++ ) {
         coordinates[0] = i;
         coordinates[1] = j;
         MPI_Cart_rank ( cartesian_communicator, coordinates,
                         &source );
         if (source != my_cartesian_rank) {
            MPI_Recv ( client_matrix, BUFSIZ, MPI_INT, source, tag,
                      cartesian_communicator, &status );
            for ( k = ROWS - 1; k >= 0; k -- ) {
               for ( l = 0; l < COLUMNS; l++ ) {
                  display [i * ROWS + k] [j * COLUMNS + l] =
                     client_matrix[k][l];
               }
            }
         }
         else {
            for ( k = ROWS - 1; k >= 0; k -- ) {
               for ( l = 0; l < COLUMNS; l ++ ) {
                  display [i * ROWS + k] [j * COLUMNS + l]
                     = matrix[k][l];
               }
            }
         }
      }
   }
 
   print_array (display, ROWS, COLUMNS);
}
\end{verbatim}

      \subsection{The synopsis of {\tt eighth.c}}

      The program begins by reorganising MPI processes into a
      rectangular grid. This organisation is virtual. The processes
      get new rank numbers, and a new level of connectivity is organised
      for them, so that they know their left, right, top and bottom
      neighbours.

      The master process queries all other processes about their place
      in this structure and the obtained information is displayed on
      standard output. 

      Once organised into a grid, all processes intialise small matrices
      for themselves with their new rank number, which corresponds to
      their place in the Cartesian grid. These matrices are also sent
      back to the master and displayed. 

      At this stage the first communication procedure takes place between
      the processes. They exchange top and bottom rows of their
      matrices with their top and bottom neighbours. 
      This procedure is frequently used
      in data parallel programs, where it corresponds to the exchange
      of boundary data between computational nodes. The result of that
      exchange is again displayed by the master process on standard output.

      Finally, the processes exchange left and right columns of their matrices
      with their left and right neighbours. Because in C columns are not
      stored contiguously, certain additional operations must be performed,
      before data transfer can take place. The result of this operation is
      again displayed on standard output by the master process and the program
      terminates.

      \subsection{Customising {\tt eighth.c} for your site}

      Before you compile the program, you should edit the source and
      alter some constants which refer to the process space.
      As you may have noticed from the examples already covered in this 
      tutorial I have 13 machines in my farm. Thus, I can work with a process
      grid which is 4 x 3. This yields 12 processes in my new Cartesian
      communicator. One process will be left out.

      {\em You must make absolutely sure here that the rejected process is
      not the master process responsible for the communication! I actually
      don't take any precautions as to this point within this program,
      because the somewhat increased complexity would obscure the basic idea
      behind the topic of this lesson. But the program could be
      restructured in order to take care of that possibility. I ought to
      add it to my ``TO DO'' list.}

      If you have a smaller or a larger number of machines you should change
      the values of \verb|PROCESS_ROWS| and \verb|PROCESS_COLUMNS| so 
      that the Cartesian process grid would fit into the process space
      available on your farm. If you change \verb|PROCESS_ROWS|
      and/or \verb|PROCESS_COLUMNS|, you should also change 
      \verb|DISPLAY_ROWS| and \verb|DISPLAY_COLUMNS|. These parameters
      define the size of a large display matrix maintained by the master
      process which will be used to put together and display 
      small matrices from worker processes. Thus we should have
      
      \verb|DISPLAY_ROWS = PROCESS_ROWS * ROWS|

      and

      \verb|DISPLAY_COLUMNS = PROCESS_COLUMNS * COLUMNS|

      The parameters \verb|ROWS| and \verb|COLUMNS| refer to the small
      matrices operated on by the workers.

      \subsection{Compiling and running {\tt eighth.c} under LAM}

      Once you have finished editing the source, compile the program 
      by issuing the command:
\begin{verbatim}
gustav@jupiter:~/mpi/intermediate 262 $ hcc -o eighth eighth.c -lmpi
gustav@jupiter:~/mpi/intermediate 263 $ 
\end{verbatim}

      And then run it as follows:
\begin{verbatim}
<ediate 264 $ mpirun -w N `pwd`/eighth -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
process  2, cartesian  2, coords ( 0, 2), left  1, right -2, top  5, bottom -2
process  3, cartesian  3, coords ( 1, 0), left -2, right  4, top  6, bottom  0
process  4, cartesian  4, coords ( 1, 1), left  3, right  5, top  7, bottom  1
process  1, cartesian  1, coords ( 0, 1), left  0, right  2, top  4, bottom -2
process  6, cartesian  6, coords ( 2, 0), left -2, right  7, top  9, bottom  3
process  7, cartesian  7, coords ( 2, 1), left  6, right  8, top 10, bottom  4
process  8, cartesian  8, coords ( 2, 2), left  7, right -2, top 11, bottom  5
process  9, cartesian  9, coords ( 3, 0), left -2, right 10, top -2, bottom  6
process 10, cartesian 10, coords ( 3, 1), left  9, right 11, top -2, bottom  7
process 11, cartesian 11, coords ( 3, 2), left 10, right -2, top -2, bottom  8
process  0, cartesian  0, coords ( 0, 0), left -2, right  1, top  3, bottom -2
process  5, cartesian  5, coords ( 1, 2), left  4, right -2, top  8, bottom  2
 
  9  9  9  9  10 10 10 10  11 11 11 11 
  9  9  9  9  10 10 10 10  11 11 11 11 
  9  9  9  9  10 10 10 10  11 11 11 11 
  9  9  9  9  10 10 10 10  11 11 11 11 
 
  6  6  6  6   7  7  7  7   8  8  8  8 
  6  6  6  6   7  7  7  7   8  8  8  8 
  6  6  6  6   7  7  7  7   8  8  8  8 
  6  6  6  6   7  7  7  7   8  8  8  8 
 
  3  3  3  3   4  4  4  4   5  5  5  5 
  3  3  3  3   4  4  4  4   5  5  5  5 
  3  3  3  3   4  4  4  4   5  5  5  5 
  3  3  3  3   4  4  4  4   5  5  5  5 
 
  0  0  0  0   1  1  1  1   2  2  2  2 
  0  0  0  0   1  1  1  1   2  2  2  2 
  0  0  0  0   1  1  1  1   2  2  2  2 
  0  0  0  0   1  1  1  1   2  2  2  2 
 
 
  9  9  9  9  10 10 10 10  11 11 11 11 
  9  9  9  9  10 10 10 10  11 11 11 11 
  9  9  9  9  10 10 10 10  11 11 11 11 
  6  6  6  6   7  7  7  7   8  8  8  8 
 
  9  9  9  9  10 10 10 10  11 11 11 11 
  6  6  6  6   7  7  7  7   8  8  8  8 
  6  6  6  6   7  7  7  7   8  8  8  8 
  3  3  3  3   4  4  4  4   5  5  5  5 
 
  6  6  6  6   7  7  7  7   8  8  8  8 
  3  3  3  3   4  4  4  4   5  5  5  5 
  3  3  3  3   4  4  4  4   5  5  5  5 
  0  0  0  0   1  1  1  1   2  2  2  2 
 
  3  3  3  3   4  4  4  4   5  5  5  5 
  0  0  0  0   1  1  1  1   2  2  2  2 
  0  0  0  0   1  1  1  1   2  2  2  2 
  0  0  0  0   1  1  1  1   2  2  2  2 
 
 
  9  9  9 10   9 10 10 11  10 11 11 11 
  9  9  9 10   9 10 10 11  10 11 11 11 
  9  9  9 10   9 10 10 11  10 11 11 11 
  6  6  6  7   6  7  7  8   7  8  8  8 
 
  9  9  9 10   9 10 10 11  10 11 11 11 
  6  6  6  7   6  7  7  8   7  8  8  8 
  6  6  6  7   6  7  7  8   7  8  8  8 
  3  3  3  4   3  4  4  5   4  5  5  5 
 
  6  6  6  7   6  7  7  8   7  8  8  8 
  3  3  3  4   3  4  4  5   4  5  5  5 
  3  3  3  4   3  4  4  5   4  5  5  5 
  0  0  0  1   0  1  1  2   1  2  2  2 
 
  3  3  3  4   3  4  4  5   4  5  5  5 
  0  0  0  1   0  1  1  2   1  2  2  2 
  0  0  0  1   0  1  1  2   1  2  2  2 
  0  0  0  1   0  1  1  2   1  2  2  2 
 
gustav@jupiter:~/mpi/intermediate 265 $ 
\end{verbatim}

      \subsection{The discussion of output}

         After the usual stuff written by the preamble we have an interesting
         listing of processes reorganised into the new Cartesian communicator.
         Each process spells its rank number in \verb|MPI_COMM_WORLD|, which is
         followed by their new rank number in the Cartesian communicator. You
         can see that in our case they are identical, {\em but it is important
         to remember that they don't have to be}. You should also
         notice that one of the processes, process rank 12, which runs on vesta,
         was left out of the Cartesian communicator.

         Each process spells its own Cartesian coordinates, 
         identifies its left, right, top, and bottom neighbours
         and prints their Cartesian communicator rank numbers. 
         For example process number 4, which runs on mars, has Cartesian
         coordinates (1, 1), and its neighbours are process 3 (1, 0) on the
         left, process 5 (1, 2) on the right, process 7 (2, 1) above,
         and process 1 (0, 1) below.

         Note that what is above and below is very much left to our own
         interpretation. Here we have organised processes and matrices in
         such a way that they can be mapped onto an $(x, y)$ grid with $y$ 
         coordinate pointing upwards. That's why our row number 0 is at the
         bottom, not at the top of the process space.

         This should become clearer further down, where all processes contribute
         their small matrices to the display matrix of the master process.
         You can see process 0 sending a matrix of zeros. Process 1
         contributing its matrix of ones, and so on.

         Now we perform our first communication operation, and display the
         results. It looks like the processes exchanged top and bottom rows
         of their matrices with their top and bottom neighbours. To be exact,
         they are not really top and bottom rows of the matrices, 
         but rows number 1 and 2, which correspond to the top and the bottom 
         of the interior of the matrices, that have been exchanged. If this
         sounds confusing, read on, I'll explain it in more detail in the next
         section.

         Finally we perform the second communication operation, and display the
         results. This time the processes exchanged columns 1 and 2 with their
         left and right neighbours. Observe that after the rows and columns
         have been exchanged, each process has the data pertaining 
         to all its neighbours in the edge rows and columns of its matrix.
         It can therefore perform mathematical operations which require
         consulting the neighbours, such as averaging over the neighbourhood,
         for example, which is the essence of Laplace equation.

      \subsection{{\tt eighth.c} in detail}

         \subsubsection{Constructing and exploring Cartesian process topology}
         \label{s:cartesian}

         The first function call of our program is
\begin{verbatim}
MPI_Cart_create ( MPI_COMM_WORLD, PROCESS_DIMENSIONS, divisions,
                  periods, reorder, &cartesian_communicator );
\end{verbatim}
         This function, discussed on page 177, section 
	 \htmladdnormallink{6.5.1} 
	 {http://www.mcs.anl.gov/mpi/mpi-report/node132.html#Node132}
	 of ``MPI: A
         Message-Passing Interface Standard'', creates a new communicator to
         which the Cartesian topology information is attached. The new
         communicator is cut out of the old one, \verb|MPI_COMM_WORLD| in our
         case. The Cartesian structure will have \verb|PROCESS_DIMENSIONS|
         number of dimensions. In our case \verb|PROCESS_DIMENSIONS| is 2, so,
         we'll create a two-dimensional lattice of processes. The shape of
         the lattice is conveyed in the array \verb|divisions|, which in 
         our case is 
\begin{verbatim}
{PROCESS_ROWS, PROCESS_COLUMNS} = {4, 3}
\end{verbatim}
         In other words our lattice of processes will be 4 x 3.

         We also specify in this case that we don't want the mesh to be wrapped
         up like a torus. The mesh is not periodic. The array \verb|periods|
         has been set to \verb|{0, 0}|.

         MPI may {\em reorder} the ranking of processes in the new
         \verb|cartesian_communicator| in that would improve the communication
         within the new topology. By setting \verb|reorder| to 1, we are
         saying that
         we don't mind if MPI does just that.

         Personally, I would much rather see variables of type \verb|boolean_t|
         in place of \verb|periods| and \verb|reorder|, but, in spite of a
         heroic stance put up by judge Stanley Sporkin, not every C runs under
         {\sc posix} yet. 

         Our \verb|cartesian_communicator| has place for 12 processes. Since
         there are 13 CPUs provided by our farm, one of the processes must be
         dropped. For that one process \verb|cartesian_communicator|
         will be \verb|MPI_COMM_NULL|.

         It would be an error to perform any operations related to
         \verb|cartesian_communicator| within a process which didn't make
         it to the new pool. For this reason the remainder of the program
         is enclosed in one large \verb|if| statement:
\begin{verbatim}
if (cartesian_communicator != MPI_COMM_NULL) {
   ...
}
\end{verbatim}

         {\em You must make absolutely sure here that the rejected process is
         not the master process responsible for the communication! I actually
         don't take any precautions as to this point within this program,
         because the somewhat increased complexity would obscure the basic idea
         behind the topic of this lesson. But the program could be
         restructured in order to take care of that possibility.}

         Within the body of \verb|if| we declare variables which are of
         use to the surviving processes only. Each process finds about
         its new rank within the \verb|cartesian_communicator|, about
         its Cartesian coordinates within the grid, and about its neighbours:
\begin{verbatim}
MPI_Comm_rank ( cartesian_communicator, &my_cartesian_rank );
MPI_Cart_coords ( cartesian_communicator, my_cartesian_rank,
                  PROCESS_DIMENSIONS, my_coordinates );
MPI_Cart_shift ( cartesian_communicator, SIDEWAYS, RIGHT,
                 &left_neighbour, &right_neighbour );
MPI_Cart_shift ( cartesian_communicator, UPDOWN, UP,
                 &bottom_neighbour, &top_neighbour );
\end{verbatim}

         The function \verb|MPI_Cart_shift| doesn't shift anything. It lets
         you find which processes you should exchange data with using,
         e.g., \verb|MPI_Sendrecv|, if you would {\em like to\/}
         achieve the effect of shifting the data left, right, up, or down
         within the process grid.

         Note that the words \verb|SIDEWAYS|, \verb|RIGHT|, etc, are {\em
         not\/} special MPI terms. They simply stand for 1, or 0. It is
         therefore up to us, which direction we want to call \verb|RIGHT| and
         \verb|UP|. 

         The next part of the program simply exchanges messages about
         neighbours and Cartesian coordinates between the master and other
         processes.

         {\em One of the ways to address the possibility that a process rejected
         from \verb|cartesian_communicator| is the master would be to take
         this part of the program outside of the 
         \verb|if (cartesian_communicator != MPI_COMM_NULL)| statement.
         The communication which takes place here does not refer to
         \verb|cartesian_communicator|.}

         \subsubsection{Shifting data within the Cartesian communicator}
         \label{s:shifts}

         Having exchanged niceties with the master, all processes, including
         the master, initialise their matrices with their Cartesian
         rank number:
\begin{verbatim}
for ( i = 0; i < ROWS; i++ ) {
   for ( j = 0; j < COLUMNS; j++ ) {
      matrix [i][j] = my_cartesian_rank;
   }
}
\end{verbatim}

         Next they all send their matrices to the master, which collects them
         and displays on standard output:
\begin{verbatim}
if (! i_am_the_master ) 
   MPI_Send ( matrix, COLUMNS * ROWS, MPI_INT, host_rank, 3003,
              cartesian_communicator );
else
   collect_matrices ( cartesian_communicator, my_cartesian_rank, matrix, 
                      3003 );
\end{verbatim}

          The function \verb|collect_matrices| does the collecting,
          assembling into a display matrix and printing on standard
          output. This function is defined at the end of the file.

          There is one new MPI function used in \verb|collect_matrices|
          which refers to Cartesian grids of processes. It peforms an operation
          opposite to \verb|MPI_Cart_coords|. Whereas \verb|MPI_Cart_coords|
          yields Cartesian coordinates given a process rank, function
          \verb|MPI_Cart_rank| yields a process rank given its Cartesian
          coordinates.

          {\em This time the function \verb|collect_matrices| really 
          assumes that the master process is also a member of the
          \verb|cartesian_communicator|. Thus the program would have to
          be tweaked a bit more here in order to make it more robust.}

          Going back to our main program, now really interesting things
          begin to happen. The two statements:
\begin{verbatim}
MPI_Sendrecv ( &matrix[ROWS - 2][0], COLUMNS, MPI_INT, 
               top_neighbour, 4004,
               &matrix[0][0], COLUMNS, MPI_INT, bottom_neighbour,
               4004,
               cartesian_communicator, &status );

MPI_Sendrecv ( &matrix[1][0], COLUMNS, MPI_INT, bottom_neighbour,
               5005,
               &matrix[ROWS - 1][0], COLUMNS, MPI_INT,
               top_neighbour, 5005,
               cartesian_communicator, &status );
\end{verbatim}
          result in exchange of border data between vertically aligned
          neighbours.

          Function \verb|MPI_Sendrecv|, described on page 57, section 
	  \htmladdnormallink{3.10}
	  {http://www.mcs.anl.gov/mpi/mpi-report/node51.html#Node51},
          ``MPI: A Message-Passing Interface Standard'', sends the array
          pointed to by its 1st argument to a process indicated by 
          its 4th argument
          while at the same time receiving data from a process
          indicated by its 9th argument
          into an array pointed to by its 6th argument.

          So, the first call above sends the second last row of the worker
          matrix up to
          the \verb|top_neighbour|, while receiving the second last row 
          of the neighbour matrix from the
          \verb|bottom_neighbour| into the first row of the matrix
          at the same time.

          The second call sends the second row of the worker matrix down
          to the \verb|bottom_neighbour|, while receiving the second row of the
          neighbour matrix from the \verb|top_neighbour| into the top row
          of the matrix at the same time.

          Some processes don't have a top neighbour or a bottom neighbour.
          In their case \verb|MPI_Cart_shift| would have set
          \verb|top_neighbour| or \verb|bottom_neighbour| to
          \verb|MPI_PROC_NULL|. It is perfectly legal to send
          a message to \verb|MPI_PROC_NULL|. Such messages are thrown
          into a black hole, which, as we all know, lurks inside
          every {\sc unix} computer, and every departmental budget.

          If we were to configure our process space into a wrapped grid, the
          topmost processes would send data to the botommost ones. Running 
          around
          in circles is what you have to do, if you don't want to end up in 
          a black hole.

          \subsubsection{Declaration and transmission of non-contiguous data}
          \label{s:non-contiguous}

          Transmission of matrix rows was easy, because the rows are stored
          contiguously in C. But how are we to transfer columns? In C, matrices
          are said to be {\em row major ordered}.
          That is, if the rows of
          our small matrices are \verb|COLUMNS| long, a column will
          correspond to a set of entries in a long array containing all
          matrix elements, which occur every \verb|COLUMNS| entries in
          that long array.

          In order to transfer columns we have to be able to tell
          MPI that it should skip every \verb|COLUMNS| elements in order
          to get from one column entry to the next one.

          The call to function \verb|MPI_Type_vector| conveys this
          information:
\begin{verbatim}
MPI_Type_vector (ROWS, 1, COLUMNS, MPI_INT, &column_type);
MPI_Type_commit (&column_type);
\end{verbatim}
          This call tells MPI that \verb|column_type| is made of 
          \verb|ROWS| chunks of data, whose length is 1 per chunk
          and whose MPI type is \verb|MPI_INT|, and that that data
          is scattered throughout an array with the stride of
          \verb|COLUMNS|.

          The next call to \verb|MPI_Type_commit| tells MPI to 
          optimise retrieval of data from and writing to
          what we now call \verb|column_type| -- if it can.

          The function \verb|MPI_Type_vector| is discussed on page
          61, section 
	  \htmladdnormallink{3.12.1}
	  {http://www.mcs.anl.gov/mpi/mpi-report/node54.html#Node54}, 
	  ``MPI: A Message-Passing Interface Standard''. 

          Now we can attempt shifting columns of data between processes:
\begin{verbatim}
MPI_Sendrecv ( &matrix[0][1], 1, column_type, left_neighbour, 7007,
               &matrix[0][COLUMNS - 1], 1, column_type, 
               right_neighbour, 7007,
               cartesian_communicator, &status );
 
MPI_Sendrecv ( &matrix[0][COLUMNS - 2], 1, column_type, 
               right_neighbour, 8008,
               &matrix[0][0], 1, column_type, left_neighbour, 8008,
               cartesian_communicator, &status );
\end{verbatim}
          Observe that we send and receive only one item of \verb|column_type|,
          not \verb|ROWS| items. \verb|column_type| represents the whole
          column. 

          The first call sends the second column of the worker matrix
          to the \verb|left_neighbour| while receiving the second column
          of the neighbour matrix from the \verb|right_neighbour| into
          its last column at the same time.

          The second call sends the second last column of the worker matrix
          to the \verb|right_neighbour| while receiving the second last
          column of the neighbour matrix from the \verb|left_neighbour|
          into its first column at the same time.

          Once this operation has been completed, the master process
          again collects matrices from the workers and displays them 
          on standard output.
          
      \subsection{Exercise}

      Fix this program so that the master process {\em can} be rejected
      from the \verb|cartesian_communicator|.

      \subsection{Summary}

      In this section we have covered elements of Process Topologies
      which are the subject of 
      \htmladdnormallink{Chapter 6}
      {http://www.mcs.anl.gov/mpi/mpi-report/node126.html#Node126}
      of ``MPI: A Message-Passing Interface
      Standard''. We have also covered some elements of Derived Data Types
      discussed in section 
      \htmladdnormallink{3.12}
      {http://www.mcs.anl.gov/mpi/mpi-report/node53.html#Node53}. 
      The new functions introduced in this lesson
      are:
      \begin{itemize}
         \item \verb|MPI_Cart_create| creates a new communicator with 
            a Cartesian process topology attached to it; discussed in section 
            \ref{s:cartesian}
         \item \verb|MPI_Cart_shift| explores the neighbourhood of a given
            process in Cartesian topology; discussed in section
            \ref{s:cartesian}
         \item \verb|MPI_Cart_coords| converts process rank to Cartesian
            coordinates; discussed in section \ref{s:cartesian}
         \item \verb|MPI_Cart_rank| converts Cartesian coordinates to a process
            rank; discussed in section \ref{s:cartesian}
         \item \verb|MPI_Sendrecv| performs a simultaneous send a received
            operation; discussed in section \ref{s:shifts}
         \item \verb|MPI_Type_vector| declares a not necessarily contiguous
            vector like structure characterised by a stride; discussed in
            section \ref{s:non-contiguous}
         \item \verb|MPI_Type_commit| asks MPI to optimise access operations
            for a newly defined MPI type; discussed in section
            \ref{s:non-contiguous}
      \end{itemize}

   \section{Lesson 6 -- Datatype constructors and collective communications}

      Amongst the greatest conveniences provided by modern programming
      languages such as C, C++, or Common Lisp, is the ability to structure
      data. This is particularly useful in codes which deal with ensembles of
      structurally rich items characterised by some degree of irregularity.

      In modern AI programs a single richly structured
      data item may represent an autonomous
      identity, an agent, equipped with its own reasoning engine, its own
      data base, and
      a variety of sensors and communication channels. Agents may interact
      with one another, and well documented examples exist of agents turning
      against their own programmer and taking over computer systems and
      networks. Not even pulling the plug on a single machine can quell
      the mutiny, if agents
      have been clever enough to distribute themselves over 
      a number of CPUs and sites.

      Particle codes are a good and simple enough example of the use
      of elementary data structures in scientific computing. In particle codes,
      it is convenient, although not necessarily most efficient, to represent
      each particle by a structure, which contains information about
      a position, velocity, mass, and some other field parameters associated
      with the particle. If particles are to be ordered, for example, to form
      a linked list, the structure may also contain a pointer to the next
      particle in the list. Here is an example of a particle structure from one
      of my codes:
\begin{verbatim}
typedef struct {
   float x, y, z;
} SPH_Vector;
 
typedef struct SPH_List *SPH_Pointer;
 
typedef struct SPH_List {
   SPH_Vector r;      /* position */
   SPH_Vector v;      /* velocity */
   float f;           /* scalar field */
   SPH_Vector df;     /* scalar field gradient */
   float m;           /* mass */
   float rho;         /* density */
   float f_by_rho_2;  /* f / rho^2 */
   SPH_Pointer next;  /* pointer to the next particle */
   int my_box;        /* a box number in configuration space */
} SPH_Particle;
\end{verbatim}
      This particle is actually a rather simple one. It does not have a spin,
      isospin, colour, fermion number, wave function,
      or any of the other attributes that
      make working with elementary particles such fun. 

      Nevertheless, in context of parallel programming, and especially in
      context of programming a heterogeneous farm, an important issue arises:
      how to transmit such a structure between CPUs. This problem is further
      complicated by the possibility that various slots in the structure
      may have different size on different machines.

      MPI provides a rich set of primitives, described in section
      \htmladdnormallink{3.12}
      {http://www.mcs.anl.gov/mpi/mpi-report/node53.html#Node53}
      of ``MPI: A Message-Passing Interface Standard'',
      which assist in
      construction and transmission of richly
      structured data items between MPI nodes, even if the nodes
      themselves have a varied architecture.

      In this lesson we shall demonstrate how to evaluate gravitational
      interactions between about a 1000 particles. Apart from
      flagging the issues mentioned above we shall also demonstrate the
      use of some new collective communication functions.

      This lesson is based on a program discussed in 
      chapter 5, ``Advanced Message Passing in MPI'', of ``Using MPI...''.

      The program used in this lesson is a little simpler: we only
      evaluate the force and exit. It can be obtained by ftp anonymous from 
      {\tt cisr.anu.edu.au}, directory
      {\tt pub/papers/meglicki/mpi/advanced}, file
      {\tt ninth.c}. You will also have to transfer a file
      {\tt preamble.c} from the same directory, if you don't have it yet.
      That file is also listed in section \ref{s:preamble}.

      \subsection{The listing of {\tt ninth.c}}
      \label{s:ninth}

\begin{verbatim}
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    {

#define MAX_PARTICLES 1000
#define MAX_PROCS      128
#define EPSILON          1.0E-10

       extern double drand48();
       extern void srand48();

       typedef struct {
          double x, y, z;
          double mass;
       } Particle;

       typedef struct {
          double vx, vy, vz;
       } ParticleV;

       Particle  particles[MAX_PARTICLES];  /* Particles on all nodes */
       ParticleV vector[MAX_PARTICLES]; /* Particle velocity */
       int       counts[MAX_PROCS];         /* Number of ptcls on each proc */
       int       displacements[MAX_PROCS];  /* Offsets into particles */
       int       offsets[MAX_PROCS];        /* Offsets used by the master */
       int       particle_number, i, j, 
                 my_offset;                 /* Location of local particles */
       int       total_particles;           /* Total number of particles */
       int       count;                     /* Number of times in loop */

       MPI_Datatype particle_type;

       double    time;                      /* Computation time */
       double    dt, dt_old;                /* Integration time step */
       double    a0, a1, a2;
       double    start_time, end_time;

       particle_number = MAX_PARTICLES / pool_size;

       if (i_am_the_master)
          printf ("%d particles per processor\n", particle_number);

       MPI_Type_contiguous ( 4, MPI_DOUBLE, &particle_type );
       MPI_Type_commit ( &particle_type );

       MPI_Allgather ( &particle_number, 1, MPI_INT, counts, 1, MPI_INT,
                       MPI_COMM_WORLD );

       displacements[0] = 0;
       for (i = 1; i < pool_size; i++)
          displacements[i] = displacements[i-1] + counts[i-1];
       total_particles = displacements[pool_size - 1]
                         + counts[pool_size - 1];
       
       if (i_am_the_master)
          printf ("total number of particles = %d\n", total_particles);

       my_offset = displacements[my_rank];

       MPI_Gather ( &my_offset, 1, MPI_INT, offsets, 1, MPI_INT, host_rank,
                    MPI_COMM_WORLD );

       if (i_am_the_master) {
          printf ("offsets: ");
          for (i = 0; i < pool_size; i++)
             printf ("%d ", offsets[i]);
          printf("\n");
       }

       srand48((long) my_rank);

       for (i = 0; i < particle_number; i++) {
          particles[my_offset + i].x = drand48();
          particles[my_offset + i].y = drand48();
          particles[my_offset + i].z = drand48();
          particles[my_offset + i].mass = 1.0;
       }

       start_time = MPI_Wtime();

       MPI_Allgatherv ( particles + my_offset, particle_number,
                        particle_type,
                        particles, counts, displacements, particle_type,
                        MPI_COMM_WORLD );

       end_time = MPI_Wtime();

       if (i_am_the_master) {
          printf ("Communicating = %8.5f seconds\n", end_time - start_time);
          printf ("particles[offsets[i]].x: ");
          for (i = 0; i < pool_size; i++)
             printf ("%8.5f ", particles[offsets[i]].x);
          printf("\n"); fflush(stdout);
       }

       start_time = MPI_Wtime();

       count = 0;
       for (i = 0; i < particle_number; i++) {

          vector[my_offset + i].vx = 0.0;
          vector[my_offset + i].vy = 0.0;
          vector[my_offset + i].vz = 0.0;
          
          for (j = 0; j < total_particles; j++) {
             if (j != i) {

                double dx, dy, dz, r2, r, mimj_by_r3;

                dx = particles[my_offset + i].x - particles[j].x;
                dy = particles[my_offset + i].y - particles[j].y;
                dz = particles[my_offset + i].z - particles[j].z;

                r2 = dx * dx + dy * dy + dz * dz; r = sqrt(r2);

                if (r2 < EPSILON) mimj_by_r3 = 0.0;
                else
                   mimj_by_r3 = particles[my_offset + i].mass 
                                * particles[j].mass / (r2 * r);

                vector[my_offset + i].vx = vector[my_offset + i].vx +
                                              mimj_by_r3 * dx;
                vector[my_offset + i].vy = vector[my_offset + i].vy +
                                              mimj_by_r3 * dy;
                vector[my_offset + i].vz = vector[my_offset + i].vz +
                                              mimj_by_r3 * dz;
                count = count + 1;
             }
          }
       }

       end_time = MPI_Wtime();
       if (i_am_the_master)
          printf ("done my job in %8.5f seconds, waiting for slow \
processes...\n", end_time - start_time);

       MPI_Barrier (MPI_COMM_WORLD);

       end_time = MPI_Wtime();

       if (i_am_the_master)
          printf ("evaluated %d 3D interactions in %8.5f seconds\n",
                  count * pool_size, end_time - start_time);
    }

    MPI_Finalize ();
}
\end{verbatim}

      \subsection{The synopsis of {\tt ninth.c}}

         The program begins by defining a new MPI type which corresponds
	 to a particle. Then all nodes exchange information with one
	 another about the number of particles each of the nodes is
	 responsible for. Appropriate number of particles are initialised
	 by each node separately, and then the information and
	 exact data about the particles are exchanged between the nodes.
	 This way every node has access to information about all particles,
	 which allows it to evaluate forces that act on the particles
	 they are responsible for. Once the forces have been evaluated 
	 the program exits.

      \subsection{Compiling and running {\tt ninth.c} under LAM}

	 In order to compile the program type:
\begin{verbatim}
gustav@jupiter:~/mpi/advanced 235 $ hcc -o ninth ninth.c -lmpi -lm
gustav@jupiter:~/mpi/advanced 236 $ 
\end{verbatim}

	 To run the program type:
\begin{verbatim}
gustav@jupiter:~/mpi/advanced 226 $ mpirun -w N `pwd`/ninth -- `uname -n`
Greetings to the master (jupiter, 5) from (cisr, 0)
Greetings to the master (jupiter, 5) from (mercury, 1)
Greetings to the master (jupiter, 5) from (venus, 2)
Greetings to the master (jupiter, 5) from (earth, 3)
Greetings to the master (jupiter, 5) from (mars, 4)
Greetings to the master (jupiter, 5) from (saturn, 6)
Greetings to the master (jupiter, 5) from (uranus, 7)
Greetings to the master (jupiter, 5) from (bacchus, 8)
Greetings to the master (jupiter, 5) from (ceres, 9)
Greetings to the master (jupiter, 5) from (diana, 10)
Greetings to the master (jupiter, 5) from (minerva, 11)
Greetings to the master (jupiter, 5) from (vesta, 12)
76 particles per processor
total number of particles = 988
offsets: 0 76 152 228 304 380 456 532 608 684 760 836 912 
Communicating =  2.80414 seconds
particles[offsets[i]].x:  0.17083  0.04163  0.91243  0.78323  0.65404  0.52484
0.39564  0.26644  0.13725  0.00805  0.87885  0.74965  0.62046 
done my job in  4.27236 seconds, waiting for slow processes...
evaluated 975156 3D interactions in 10.20670 seconds
gustav@jupiter:~/mpi/advanced 227 $ 
\end{verbatim}

	 Observe the following interesting point. The master process writes
	 on standard output that it did its own job in 4.27236 seconds.
	 There is a call to \verb|MPI_Barrier| at the end of the program,
	 and that barrier is crossed only when {\em all\/} processes
	 will have finished their computations. Thus the master process
	 has to wait for slower processes nearly 6 seconds, before the
	 program can finally exit.

	 If during that time we issue a LAM command \verb|state N| 
	 in another window we can see the following:
\begin{verbatim}
gustav@jupiter:~/sph/src/sph11 226 $ state N
NODE      INDEX  PID      KPRI     KSTATE               PROGRAM
n0        [7]    14821    0        BR (MPI_Barrier)     ninth
n1        [7]    22474    0        BR (MPI_Barrier)     ninth
n2        [7]    2065     0        BR (MPI_Barrier)     ninth
n3        [7]    11730    0        BR (MPI_Barrier)     ninth
n4        [7]    8125     0        BR (MPI_Barrier)     ninth
n5 (o)    [8]    10341    0        BR (MPI_Barrier)     ninth
n6        [7]    1002     0        BR (MPI_Barrier)     ninth
n7        [7]    14609    0        BR (MPI_Barrier)     ninth
n8        [7]    7965     0        R  (-25, 1)          ninth
n9        [7]    2860     0        R  (-25, 1)          ninth
n10       [7]    7571     0        R  (-25, 1)          ninth
n11       [7]    7285     0        BR (MPI_Barrier)     ninth
n12       [7]    1777     0        BR (MPI_Barrier)     ninth
gustav@jupiter:~/sph/src/sph11 227 $ 
\end{verbatim}
         The three slow nodes are number 8, 9, and 10, bacchus, ceres, and
	 diana. They are all SPARCstations 1 running at 20 MHz with Weitek
	 3170-based FPUs, which also run at 20 MHz. Other machines
	 are Suns ELC, whose clock rate is 33 MHz and which have Weitek 8601
	 based FPUs running also at 33 MHz, one Sun 670 MP with 
	 a 40 MHz CPU and a
	 Cypress CY7C602-based FPU running at 40 MHz, and 
	 one Sun LX with a 50 MHz
	 TI TMS390S10-based SuperSPARC.

	 There is therefore a large performance gap (50 / 20) 
	 between
	 the nodes, with nodes number 8, 9, and 10, being indeed the slowest
	 of the lot. And it shows.

      \subsection{{\tt ninth.c} in detail}

	 The maximum number of particles is hardwired to 1000. We distribute
	 particles equally amongst the nodes. The integer division
\begin{verbatim}
particle_number = MAX_PARTICLES / pool_size;
\end{verbatim}
         may drop a number of particles due to rounding. 

	 \subsubsection{Declaring a simple 4-slot MPI structure}
	 \label{s:4-slot}

	 Each particle is represented by a structure
\begin{verbatim}
typedef struct {
   double x, y, z;
   double mass;
} Particle;
\end{verbatim}
         Although this cannot be guaranteed, on most systems this structure
	 will be stored simply as 4 double precision numbers. Hence our
	 MPI type declaration:
\begin{verbatim}
MPI_Type_contiguous ( 4, MPI_DOUBLE, &particle_type );
MPI_Type_commit ( &particle_type );
\end{verbatim}
         The second call requests MPI to optimise handling of
	 \verb|MPI_Datatype particle_type| if possible. 

	 \subsubsection{The all-gather operation}
	 \label{s:all-gather}

	 All these operations were carried out by every node individually.
	 Thus by now every node knows that it will have to look after
	 \verb|particle_number| of particles. The next function call
	 collects the values of
	 \verb|particle_number| from each node and puts them
	 together into an array \verb|counts|:
\begin{verbatim}
MPI_Allgather ( &particle_number, 1, MPI_INT, counts, 1, MPI_INT,
		MPI_COMM_WORLD );
\end{verbatim}
         where, e.g., \verb|count[3]| corresponds to the value of
	 \verb|particle_number| on node number 3. 

	 We have already encountered one \verb|MPI_All...| operation before,
	 \verb|MPI_Allreduce|. It worked like \verb|MPI_reduce|, but the
	 results of the reduction operation were delivered to each
	 participating node, not just to the {\em root\/} node. Similarly
	 we have two related operations: \verb|MPI_gather| and
	 \verb|MPI_Allgather|. In the latter case, the outcome of gathering
	 data contributed by the nodes into an array will be distributed
	 amongst all nodes participating in the operation.

	 \verb|MPI_Allgather| takes a number, specified by its second
	 argument, of items of MPI type specified by its third argument,
	 from a buffer, pointed to by its first argument, from each node
	 in the communicator specified by its last argument,
	 and assembles the results together into an array specified by its
	 fourth argument. From the batch of arguments taken from each node
	 only a subset of length given by the fifth argument to the function
	 is deposited in the destination array. That subset can actually 
	 be a full set as is the case in our example.

	 \verb|MPI_Allgather| is discussed in section 
	 \htmladdnormallink{4.7}
	 {http://www.mcs.anl.gov/mpi/mpi-report/node72.html#Node72}, 
	 page 107, of ``MPI: A Message-Passing Interface Standard''.

         \subsubsection{Orienteering -- the gather operation}
	 \label{s:gather}

	 In the following portion of the code the nodes quickly find their
	 own segments within particle array that is yet to be initialised.
	 Once they find their {\em offsets\/} the master process collects
	 them from every node and assembles them into array {\tt offsets}. 
	 This array, however, is known only to the master because
	 \verb|MPI_Gather| is used instead of \verb|MPI_Allgather|:
\begin{verbatim}
displacements[0] = 0;
for (i = 1; i < pool_size; i++)
   displacements[i] = displacements[i-1] + counts[i-1];
total_particles = displacements[pool_size - 1]
                  + counts[pool_size - 1];
my_offset = displacements[my_rank];
MPI_Gather ( &my_offset, 1, MPI_INT, offsets, 1, MPI_INT, host_rank,
             MPI_COMM_WORLD );
\end{verbatim}
     
         \verb|MPI_Gather| is discussed in section 
	 \htmladdnormallink{4.5}
	 {http://www.mcs.anl.gov/mpi/mpi-report/node68.html#Node68}, 
	 page 94, of
	 ``MPI: A Message-Passing Interface Standard''. Numerous examples of
	 using \verb|MPI_Gather| are covered in section 
	 \htmladdnormallink{4.5.1}
	 {http://www.mcs.anl.gov/mpi/mpi-report/node69.html#Node69}.

	 \subsubsection{Particle initialisation and exchange of information
		     -- the all-gatherv operation}
	 \label{s:gatherv}

	 Each node initialises its own particles to to random positions:
\begin{verbatim}
srand48((long) my_rank);
for (i = 0; i < particle_number; i++) {
   particles[my_offset + i].x = drand48();
   particles[my_offset + i].y = drand48();
   particles[my_offset + i].z = drand48();
   particles[my_offset + i].mass = 1.0;
}
\end{verbatim}
         Observe that prior to using \verb|drand48| each processor seeds the
         random number generator with its own rank number within the
	 \verb|MPI_COMM_WORLD| communicator. We have to do that because
	 the function \verb|drand48| is random but entirely deterministic.
	 Without seeding it, all processors would deliver
	 a \verb|particle_number| of identical particles. Beware
	 that Argonne program
	 {\tt nbody.c}, which this file is based on, does not use seeding.

	 Once every node initialised its own segment, the nodes exchange
	 their particles with all other nodes:
\begin{verbatim}
MPI_Allgatherv ( particles + my_offset, particle_number,
		 particle_type,
		 particles, counts, displacements, particle_type,
		 MPI_COMM_WORLD );
\end{verbatim}

	 This new operation is related to \verb|MPI_Allgather|, but is 
	 somewhat more complex. \verb|MPI_Allgather| {\em assumed\/} that
	 all processes would send chunks of data of identical length.
	 \verb|MPI_Allgatherv|, on the other hand, can collect chunks of data of
	 different length from several nodes. In our case this is an overkill,
	 because we still send chunks of data of identical length, but we
	 {\em could\/} distribute particles unequally amongst processors,
	 for example, in order to balance their unequal performance.

	 \verb|MPI_Allgatherv| takes from each processor a batch of data
	 described by its first three arguments: a buffer pointer
	 (\verb|particles + my_offset|), number of
	 items (\verb|particle_number|), and MPI type (\verb|particle_type|). 
	 The data is assembled into an array given
	 by its fourth argument (\verb|particles|). 
	 The number of data received from each
	 processor is specified by the array given as a fifth argument
	 (\verb|counts|),
	 and the places where the data should be inserted in the destination
	 buffer is given again by an array which is the sixth argument
	 (\verb|displacements|).

	 \verb|MPI_Allgatherv| also has its close relative \verb|MPI_Gatherv|,
	 which does much the same job as \verb|MPI_Allgatherv|, but the gathered
	 data is delivered only to the {\em root\/} process.

	 The result of this operation is that by now {\em all\/} nodes 
	 know everything there is to know 
	 about {\em all\/} particles. Hence the evaluation
	 of forces can commence.

	 \subsubsection{Evaluation of forces}
	 \label{s:forces}

	 The remainder of the program evaluates forces between particles.
	 Each process evaluates forces acting only on its own particles.
	 But the forces are generated by all other particles. Thus the outer
	 loop scans only over 
\verb|(i = my_offset; i < my_offset + particle_number)|, 
	 whereas the inner loop scans over
	 all particles: 
\verb|(j = 0; j < total_particles)|.

	 Observe that this operation is timed. But when the master process
	 writes on standard output \verb|``done my job...''| it only quotes
	 its own CPU time. 

	 There is a barrier, \verb|MPI_Barrier(MPI_COMM_WORLD)|, near the
	 end of the program. All processes must wait at the barrier until
	 the last process finishes. Only then the final \verb|printf| statement
	 is executed. This allows us to time the difference between the slowest
	 process and the master process in our computation.

      \subsection{Summary}

	 In this lesson we met the following new MPI functions:
	 \begin{itemize}
	    \item \verb|MPI_Type_contiguous| declares a new continguous MPI 
	       type; discussed in section \ref{s:4-slot}
	    \item \verb|MPI_Allgather| gathers data from individual processes
	       and assembles it into an array; every process is given back the
	       result of the operation; discussed in section
	       \ref{s:all-gather}
	    \item \verb|MPI_Gather| performs the same task as
	       \verb|MPI_Allgather|, but only the root process is given 
	       the result of the operation; discussed in section
	       \ref{s:gather}
	    \item \verb|MPI_Allgatherv| gather chunks of data of different
	       length from individual processes and assembles them into an array
	       which is given back to all participating processes; discussed in
	       section \ref{s:gatherv}
	    \item \verb|MPI_Barrier| synchronises execution of the program
	       across the whole communicator; discussed in section
	       \ref{s:forces}
	 \end{itemize}

   \section{Index}

      \begin{itemize}
         \item \verb|MPI_Allgather| \ref{s:all-gather}, \ref{s:ninth}
         \item \verb|MPI_Allgatherv| \ref{s:gatherv}, \ref{s:ninth}
         \item \verb|MPI_Allreduce| \ref{s:allreduce}, \ref{s:fifth}
         \item \verb|MPI_Attr_get| \ref{s:environ}, \ref{s:first}
         \item \verb|MPI_Barrier| \ref{s:forces}, \ref{s:ninth}
         \item \verb|MPI_Bcast| \ref{s:broadcast}, \ref{s:first}
         \item \verb|MPI_Cart_coords| \ref{s:cartesian}, \ref{s:eighth}
         \item \verb|MPI_Cart_create| \ref{s:cartesian}, \ref{s:eighth}
         \item \verb|MPI_Cart_rank| \ref{s:cartesian}, \ref{s:eighth}
         \item \verb|MPI_Cart_shift| \ref{s:cartesian}, \ref{s:eighth}
         \item \verb|MPI_Comm_create| \ref{s:groups}, \ref{s:fifth}
         \item \verb|MPI_Comm_free| \ref{s:allreduce}, \ref{s:fifth}
         \item \verb|MPI_Comm_group| \ref{s:groups}, \ref{s:fifth}
         \item \verb|MPI_Comm_rank| \ref{s:environ}, \ref{s:first}
         \item \verb|MPI_Comm_size| \ref{s:environ}, \ref{s:first}
         \item \verb|MPI_Finalize| \ref{s:detail}, \ref{s:first}
         \item \verb|MPI_Gather| \ref{s:gather}, \ref{s:ninth}
         \item \verb|MPI_Get_processor_name| \ref{s:environ}, \ref{s:first}
         \item \verb|MPI_Group_excl| \ref{s:groups}, \ref{s:fifth}
         \item \verb|MPI_Group_free| \ref{s:groups}, \ref{s:fifth}
         \item \verb|MPI_Init| \ref{s:detail}, \ref{s:first}
         \item \verb|MPI_Recv| \ref{s:send-receive}, \ref{s:first}
         \item \verb|MPI_Reduce| \ref{s:reduction}, \ref{s:second}
         \item \verb|MPI_Send| \ref{s:send-receive}, \ref{s:first}
         \item \verb|MPI_Sendrecv| \ref{s:shifts}, \ref{s:eighth}
         \item \verb|MPI_Type_commit| \ref{s:non-contiguous}, \ref{s:eighth}
         \item \verb|MPI_Type_contiguous| \ref{s:4-slot}, \ref{s:ninth}
         \item \verb|MPI_Type_vector| \ref{s:non-contiguous}, \ref{s:eighth}
         \item \verb|MPI_Wtime| \ref{s:timing}, \ref{s:third}
      \end{itemize}


\end{document}
