Subject: Example attempt...

From: [email-masked] (Tom Henderson)
Date: Wed, 7 Jul 93 16:59:01 MDT

Hi all,

At the last meeting I promised to try and find some "more real" program 
examples and "re-write" them using MPI contexts.  Here's one...  It is the 
first application I thought of when I first heard the "group" concept.  I 
thought it would be "real simple" but it's gotten a bit long.  At the very 
least I'm learning a lot about the current draft!  :-)

ISAR EXAMPLE

This example is taken from an old ISAR (Inverse Synthetic Aperture Radar) 
translational motion compensation and display program that I originally wrote 
in Occam long, long ago.  The program converts raw radar data into a human-
readable sequence of animated images and displays them in "real time" on a 
monitor like this:  

            +--------+        +------------+        +---------+
raw         |        |        |            |        |         |
radar --->  | SEARCH |  --->  | COMPENSATE |  --->  | DISPLAY |  --->  monitor
data        |        |        |  (FFT's)   |        |         |
            +--------+        +------------+        +---------+

The program divides naturally into three tasks:  search, compensate, and 
display.  The search task uses an iterative algorithm to determine the current 
motion of the target relative to the radar set (distance, velocity, and 
acceleration).  The compensation task modifies phase in the raw data according 
to the current motion and uses FFT's to convert from frequency space to 
"image" space.  The display task does various scaling, thresholding, 
centering, etc. and displays the resulting image on a monitor.  Each time 
new radar data arrives, a new image is generated.  

The search task must pass current raw radar data and current motion parameters 
to the compensation task.  The compensation task must pass the raw image to 
the display task.  There is no feedback loop in the program, this is a pure 
pipeline.  

To be scalable, the program must be able to choose how many processes are 
allocated to each task at run time.  The easiest way to do this is to assign a 
separate sub-group to each task.  The sub-groups can then be populated at run 
time in some "sensible" way to achieve load balance.  

MY ATTEMPT AT AN MPI VERSION

The "pseudo code" below is an attempt to set up the sub-groups and get them 
communicating.  Processing details are mostly left hidden inside functions.  
I've left a few details exposed so the "flavor" of the original code is still 
there.  Some details:  

    1)    I have avoided using publish/subscribe because we have already seen 
          an example of this in the last meeting.  Actually, this code is 
          simpler if MPI_COMM_PUBLISH_SUBSCRIBE() is used.  If there is any 
          interest, I'll send out the version that uses publish/subscribe 
          later.  

    2)    Communication in the pipeline is as simple as possible.  I have not 
          used double or triple buffering to improve performance.  If I figure 
          the "simple" example out, I'll add optimizations.  

    3)    The compensation sub-group calls routines from a third-party FFT 
          library.  I put this detail in because I'm not quite sure how it all 
          works yet.  For now, I'm assuming that this library can behave like 
          the MPI collective communication calls and that it requires an 
          initialization call (INIT_FFT_LIB() in the code below).  

    4)    When I have had a choice of MPI functions to use, I have attempted 
          to use the "highest level" function.  

(It helps to look at the code when reading the following description...)  

The code starts with the usual #includes (  :-)  ) and #defines.  The include 
file for the FFT library is fftLib.h.  Three keys are #defined for sub-group 
creation using MPI_COLL_SUBGROUP().  MAX_FLAT_GROUP_SIZE is a number pulled 
out of thin air as maximum size (in bytes) of a "flattened" group.  

Two structures are typedef'ed.  The isarParameters structure contains all the 
information required by the compensation sub-group to compensate raw radar 
data.  The most interesting fields are the motion fields (r, v, a).  These are 
filled in for each image by the search sub-group.  A decompParameters 
structure contains information about how data is decomposed in a sub-group.  
It is used to determine how data is re-organized when arrays are passed from 
one sub-group to the next.  

After MPI_INIT(), INIT_FFT_LIB(), etc., function getGroupSizes() is called to 
determine how many processes will be assigned to each sub-group.  This 
function also sets the membership key for each process.  MPI_COLL_SUBGROUP() 
is then called using the membership key to create all three sub-groups.  

Every process is then made aware of each sub-group.  Process 0 in each 
sub-group sends a flattened group description to process 0 in the "ALL" 
group.  The flattened groups are then broadcast to all processes and 
unflattened.  

Communicators for each sub-group are then created using MPI_COMM_MAKE().  
Intergroup communicators for each of the communicating subgroups are then 
created using MPI_COMM_MERGE().  

Decomposition parameters are then set up for each sub-group using function 
makeDecompParams().  (MPI topology functions would probably live inside 
makeDecompParams().)    After "other" set up stuff, the main loops for each 
sub-group begin.  

In the main loop for the search sub-group, function addNewData() gets new raw 
radar data and stores it in array radarData.  Function findNewMotion() 
determines target motion and updates the (r, v, a) fields of isarParams.  
Function sendCurrentData() sends radarData to the compensate sub-group (with 
appropriate re-organization according to data decompositions described by 
searchDecomp and compensateDecomp).  Finally, function sendMotion broadcasts 
the (r, v, a) fields of isarParams to every process in the compensate 
sub-group.  

In the main loop for the compensate sub-group, function recvCurrentData() gets 
updated raw radar data from the search sub-group.  Function recvMotion() then 
gets the current (r, v, a) fields of isarParams from the search sub-group.  
Function compensateData() uses (r, v, a) to apply phase corrections to the 
radar data.  It then calls the FFT routines from the third-party library to 
create the raw image.  Function sendCurrentData() passes the raw image to the 
display sub-group.  

In the main loop for the display sub-group, function recvCurrentData() gets 
the raw image from the compensate sub-group.  Function scaleThresholdCenter() 
produces the human-readable image.  Function displayImage stuffs the image 
into the monitor.  

SOME QUESTIONS

1)    Am I doing the context/sub-group creation in the "best" way (given that 
      I did not use publish/subscribe)?  Is this even correct?  

2)    Assume that the third-party FFT library makes collective communication 
      calls.  Will it work as written?  Do I need to pass any "extra" contexts 
      or communicators into compensateData() to make things work?  Do I really 
      need INIT_FFT_LIB()?  

3)    Does MPI_COMM_MERGE() create new contexts?  If not, is it possible to 
      get confused between communications that use (for example) 
      compensateComm and sendComm?  This question comes up when I want to add 
      the double/triple buffering.  

Please fire away...  If there is any interest, I'll pass on the 
publish/subscribe version and maybe even work on a double/triple buffering 
version...  

I hope I don't find too many stupid mistakes when I read this tomorrow!  :-)

Tom Henderson
NOAA Forecast Systems Laboratory

------------------------------ BEGIN "PSEUDO" CODE ---------------------------

/* #includes */
#include <mpi.h>       /* :-) */
#include <fftLib.h>    /* For FFT library. */
...

/* #defines */

/* Membership keys for sub-group creation. */
#define SEARCH_GROUP_KEY      0
#define COMPENSATE_GROUP_KEY  1
#define DISPLAY_GROUP_KEY     2

/* Buffer size for flattened groups. */
#define MAX_FLAT_GROUP_SIZE   4096                            /* whatever... */

...

/* typedefs */

typedef struct isarParameters                             /* ISAR parameters */
    {
    int num_bursts;            /* number of bursts (rows) in each ISAR frame */
    int num_pulses;         /* number of pulses (columns) in each ISAR frame */
    float r;                                              /* target range: m */
    float v;                                         /* target velocity: m/s */
    float a;                                   /* target acceleration: m/s^2 */
    float prf;                                  /* pulse repetition rate: Hz */
    float deltaf;                                 /* frequency step size: Hz */
    float f0;                                         /* start frequency: Hz */
    } isarParameters;

typedef struct decompParameters             /* data decomposition parameters */
    {
    int numRows;                               /* number of rows in 2D array */
    int numCols;                            /* number of columns in 2D array */
    int myNumRows;                        /* number of rows in my 2D "chunk" */
    int myNumCols;                     /* number of columns in my 2D "chunk" */
    int myStartRow;                            /* start row of my 2D "chunk" */
    int myStartCol;                         /* start column of my 2D "chunk" */
    int numProcRows;                  /* number of logical rows of processes */
    int numProcCols;               /* number of logical columns of processes */
    int myProcRow;                                 /* my logical process row */
    int myProcCol;                              /* my logical process column */
    } decompParameters;

main()
    {
    group searchGroup, compensateGroup, displayGroup, myGroup;
    comm searchComm, compensateComm, displayComm, sendComm, recvComm, tmpComm;
    float **radarData;
    char **image;
    isarParameters isarParams;
    decompParameters searchDecomp, compensateDecomp, displayDecomp;
    int membershipKey, myRank, myRankInAll, length;
    int allGroupSize, searchGroupSize, compensateGroupSize, displayGroupSize;
    char myFlatGroup[MAX_FLAT_GROUP_SIZE], 
         searchFlatGroup[MAX_FLAT_GROUP_SIZE], 
         compensateFlatGroup[MAX_FLAT_GROUP_SIZE], 
         displayFlatGroup[MAX_FLAT_GROUP_SIZE];
    comm_handle searchHandle, compensateHandle, displayHandle;
    return_handle status;

/* Generic setup. */
    MPI_INIT();
    MPI_COMM_SIZE(MPI_COMM_ALL, &allGroupSize);
    MPI_COMM_RANK(MPI_COMM_ALL, &myRankInAll);

/* Initialize library called inside function compensateData(). */
    INIT_FFT_LIB(...);

/* Other setup stuff (malloc() 2D arrays, etc...) */
    ...

/* Build sub-groups. */
    getGroupSizes(allGroupSize, &searchGroupSize, &compensateGroupSize, 
                  &displayGroupSize, &membershipKey, ...);
    MPI_COLL_SUBGROUP(MPI_COMM_ALL, membershipKey, myGroup);

/* Process 0 in each sub-group sends flattened group description to process */
/* 0 in the "ALL" group.  Flattened groups are then broadcast to all */
/* processes and unflattened.  Beware that process 0 in the "ALL" group may */
/* also be process 0 in a sub-group. */
    /* Process 0 in "ALL" posts receives for flattened groups. */
    if (myRankInAll == 0)
        {
        MPI_IRECVC(searchHandle, searchFlatGroup, MAX_FLAT_GROUP_SIZE, 
                   MPI_BYTE, MPI_DONTCARE, SEARCH_GROUP_KEY, MPI_COMM_ALL);
        MPI_IRECVC(compensateHandle, compensateFlatGroup, MAX_FLAT_GROUP_SIZE, 
                   MPI_BYTE, MPI_DONTCARE, COMPENSATE_GROUP_KEY, MPI_COMM_ALL);
        MPI_IRECVC(displayHandle, displayFlatGroup, MAX_FLAT_GROUP_SIZE, 
                   MPI_BYTE, MPI_DONTCARE, DISPLAY_GROUP_KEY, MPI_COMM_ALL);
        }
    /* Process 0 in sub-group flattens group and sends to process 0 in */
    /* "ALL". */
    MPI_GROUP_RANK(myGroup, &myRank);
    if (myRank == 0)
        {
        MPI_GROUP_FLATTEN(myGroup, MAX_FLAT_GROUP_SIZE, myFlatGroup, &length);
        if (length > MAX_FLAT_GROUP_SIZE)
            {
            handleError(...);
            }
        MPI_SENDC(myFlatGroup, MAX_FLAT_GROUP_SIZE, MPI_BYTE, 0, 
                  membershipKey, MPI_COMM_ALL);
        }
    /* Process 0 in "ALL" waits for completion of receives and then */
    /* broadcasts flat groups to all processes. */
    if (myRankInAll == 0)
        {
        MPI_WAIT(searchHandle, status);
        MPI_WAIT(compensateHandle, status);
        MPI_WAIT(displayHandle, status);
        }
    MPI_BCASTC(searchFlatGroup, MAX_FLAT_GROUP_SIZE, MPI_BYTE, 
               MPI_COMM_ALL, 0);
    MPI_BCASTC(compensateFlatGroup, MAX_FLAT_GROUP_SIZE, MPI_BYTE, 
               MPI_COMM_ALL, 0);
    MPI_BCASTC(displayFlatGroup, MAX_FLAT_GROUP_SIZE, MPI_BYTE, 
               MPI_COMM_ALL, 0);
    /* Finally, unflatten all group descriptions. */
    MPI_GROUP_UNFLATTEN(MAX_FLAT_GROUP_SIZE, searchFlatGroup, searchGroup);
    MPI_GROUP_UNFLATTEN(MAX_FLAT_GROUP_SIZE, compensateFlatGroup, 
                        compensateGroup);
    MPI_GROUP_UNFLATTEN(MAX_FLAT_GROUP_SIZE, displayFlatGroup, displayGroup);

/* Build communicators. */
    MPI_COMM_MAKE(MPI_COMM_ALL, searchGroup, searchComm);
    MPI_COMM_MAKE(MPI_COMM_ALL, compensateGroup, compensateComm);
    MPI_COMM_MAKE(MPI_COMM_ALL, displayGroup, displayComm);

/* Build intergroup communicators. */
    if (membershipKey == SEARCH_GROUP_KEY)
        {
        MPI_COMM_MERGE(searchComm, compensateComm, sendComm, tmpComm);
        }
    else if (membershipKey == COMPENSATE_GROUP_KEY)
        {
        MPI_COMM_MERGE(compensateComm, searchComm, tmpComm, recvComm);
        MPI_COMM_MERGE(compensateComm, displayGroup, sendComm, tmpComm);
        }
    else if (membershipKey == DISPLAY_GROUP_KEY)
        {
        MPI_COMM_MERGE(displayGroup, compensateComm, tmpComm, recvComm);
        }
    else
        {
        handleError(...);
        }
    MPI_COMM_FREE(tmpComm);

/* Data decomposition set-up code for all groups... */
    makeDecompParams(searchGroupSize, searchDecomp, ...);
    makeDecompParams(compensateGroupSize, compensateDecomp, ...);
    makeDecompParams(displayGroupSize, displayDecomp, ...);

/* Other stuff... */
    ...

/* MAIN LOOPS for each sub-group... */
    /* MAIN LOOP for "search" group... */
    if (membershipKey == SEARCH_GROUP_KEY)
        {
        searchDone = 0;
        while (searchDone == 0)
            {
            ...
            addNewData(searchComm, radarData, searchDecomp, ...);
            findNewMotion(searchComm, radarData, isarParams, searchDecomp, ...);
            sendCurrentData(sendComm, radarData, searchDecomp, 
                            compensateDecomp);
            sendMotion(sendComm, isarParams);
            ...
            }  /* end of searchGroup while */
        }  /* end of searchGroup if */
    /* MAIN LOOP for "compensate" group... */
    else if (membershipKey == COMPENSATE_GROUP_KEY)
        {
        compensateDone = 0;
        while (compensateDone == 0)
            {
            ...
            recvCurrentData(recvComm, radarData, compensateDecomp, 
                            searchDecomp);
            recvMotion(recvComm, isarParams);
            compensateData(compensateComm, radarData, isarParams, 
                           compensateDecomp, ...);
            sendCurrentData(sendComm, radarData, compensateDecomp, 
                            displayDecomp);
            ...
            }  /* end of compensateGroup while */
        }  /* end of compensateGroup if */
    /* MAIN LOOP for "display" group... */
    else if (membershipKey == DISPLAY_GROUP_KEY)
        {
        displayDone = 0;
        while (displayDone == 0)
            {
            ...
            recvCurrentData(recvComm, radarData, displayDecomp, 
                            compensateDecomp);
            scaleThresholdCenter(displayComm, radarData, image, 
                                 displayDecomp, ...);
            displayImage(displayComm, image, displayDecomp);
            ...
            }  /* end of displayGroup while */
        }  /* end of displayGroup if */
     else
        {
        handleError(...);
        }

/* Cleanup... */
    ...

    }  /* end of main() */