|| || || || || || || ||

---

PxDistribution.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2003-2004, University of Amsterdam, The Netherlands.
 *  All rights reserved.
 *
 *  Author(s):
 *  Frank Seinstra <fjseins@wins.uva.nl>
 */

#ifndef __PxDistribution_h_
#define __PxDistribution_h_


#include "Core/Array/Pattern/ArrayFunc.h"
#include "Core/Array/Pattern/PxArrayFunc.h"
#include "Core/Array/Pattern/PxPartition.h"
#include "Core/Array/Pattern/PxSystem.h"
#include "mpi.h"

#include "PxxStream.h"

namespace Impala
{
namespace Core
{
namespace Array
{
namespace Pattern
{


/*** Enumerations *****************************************************/

enum tags   { INIT_TAG, BCAST_TAG, RDUCE_TAG,   // message tags
              SCAT_TAG, GATH_TAG, RDIST_TAG, BOR_TAG };
enum dists  { PX_OFT, PX_SBT, PX_MPI };                 // distribution types
enum dirs   { TO, FROM };                       // direction types
enum planes { YZ_PLANE, XZ_PLANE, XY_PLANE };   // plane types


/*** Global variables *************************************************/

static int  _myCPU;     /* Local CPU number */
static int  _nrCPUs;    /* Total nr of CPUs */
static int  _logCPUs;   /* Shortest CPU->CPU path in hypercube */
static int  _maxCPUs;   /* Max nr of CPUs in complete hypercube */

int  dp = 0;

/**********************************************************************
 *** Distribution initialization & helper functions                 ***
 **********************************************************************/

inline void
PxInitDistribution(int aw, int ah, int ad,
                   int xcpus, int ycpus, int zcpus)
{
    _myCPU   = PxMyCPU();
    _nrCPUs  = PxNrCPUs();
    _logCPUs = log((double)_nrCPUs)/log(2.0);
    _maxCPUs = pow((double)2, _logCPUs);
    if (_maxCPUs < _nrCPUs) {
        _logCPUs++;
        _maxCPUs *= 2;
    }
    PxInitPartition(aw, ah, ad, xcpus, ycpus, zcpus);
}


template<class ArrayT>
static inline void
PxLclArrayCopy(ArrayT* loc, int dir, ArrayT* blk, int root, bool bdata)
{
    /*** Copy local array (loc) TO/FROM 'head' of block array (blk) ***/

    typedef typename ArrayT::StorType storT;
    storT* lPtr = (bdata) ? ArrayPB(loc) : ArrayCPB(loc);
    storT* bPtr = (bdata) ? ArrayPB(blk) : ArrayCPB(blk);
    int eSize = ArrayT::ElemSize();

    if (_myCPU == root) {
        bPtr = ArrayPB(blk) + eSize*(PxLclStart(_myCPU,
               ArrayBW(blk), ArrayBH(blk), ArrayBD(blk)));
        if (bdata) {
            bPtr -= (ArrayW(blk)*ArrayH(blk)*ArrayBD(blk) +
                ArrayW(blk)*ArrayBH(blk) + ArrayBW(blk)) * eSize;
        }
    }

    int width  = (bdata) ? ArrayW(loc) : ArrayCW(loc);
    int height = (bdata) ? ArrayH(loc) : ArrayCH(loc);
    int depth  = (bdata) ? ArrayD(loc) : ArrayCD(loc);

    for (int z=0; z<depth; z++) {
        for (int y=0; y<height; y++) {
            storT* locPtr = lPtr + (ArrayW(loc)*ArrayH(loc)*z +
                                    ArrayW(loc)*y) * eSize;
            storT* blkPtr = bPtr + (ArrayW(blk)*ArrayH(blk)*z +
                                    ArrayW(blk)*y) * eSize;
            for (int x=0; x<width*eSize; x++) {
                (dir == TO) ? (*blkPtr++ = *locPtr++) :
                              (*locPtr++ = *blkPtr++);
            }
        }
    }
}


static inline int
PxResponsibilityRange(int cpuNr)
{
    /*** Calculate the number of CPUs that 'cpuNr' is responsible ***/
    /*** for in hypercube (SBT) collective communication.         ***/

    int range = 1;

    for (int i=0; i<_logCPUs; i++) {
        if ((cpuNr & range) == range) {
            return range;
        }
        range <<= 1;
    }
    return range;
}


static inline int
PxLineRange(int dimension, int firstIndex, int *order)
{
    /*** Calculate number of 'lines' that will be send/received in ***/
    /*** SBT collective communication. The CPU at 'firstIndex' in  ***/
    /*** the SBTarray 'order' is responsible for a range of CPUs,  ***/
    /*** and will send and receive data for that range.            ***/

    int minNrLines = 0;
    int overflow = 0;

    switch (dimension) {
        case 1: if (PxLogPartXcpus() == 1) {    // get nr of 'columns'
                    return PxFullWidth();
                }
                minNrLines = PxMinLclWidth();
                overflow = PxOverflowX();
                break;
        case 2: if (PxLogPartYcpus() == 1) {    // get nr of 'rows'
                    return PxFullHeight();
                }
                minNrLines = PxMinLclHeight();
                overflow = PxOverflowY();
                break;
        case 3: if (PxLogPartZcpus() == 1) {    // get nr of 'pages'
                    return PxFullDepth();
                }
                minNrLines = PxMinLclDepth();
                overflow = PxOverflowZ();
                break;
        default:
                return 0;
    }

    int lastIndex = firstIndex + PxResponsibilityRange(firstIndex) - 1;
    int nrIndexes = lastIndex - firstIndex + 1;
    int firstCPU = order[firstIndex];
    if (firstIndex == 0) {
        firstCPU = (firstCPU / 4) * 4;
    }
    int lastCPU = firstCPU + nrIndexes - 1;
    if (lastCPU >= _nrCPUs) {
        lastCPU = _nrCPUs - 1;
    }

    /*** Normalize CPU numbers to grid indices in correct dimension ***/

    switch (dimension) {
        case 1: firstCPU = PxGridIndexX(firstCPU);
                lastCPU  = PxGridIndexX(lastCPU);
                break;
        case 2: firstCPU = PxGridIndexY(firstCPU);
                lastCPU  = PxGridIndexY(lastCPU);
                break;
        case 3: firstCPU = PxGridIndexZ(firstCPU);
                lastCPU  = PxGridIndexZ(lastCPU);
                break;
    }
    int nrBufCPUs = lastCPU - firstCPU + 1;
    if (firstCPU >= overflow) {
        return (minNrLines * nrBufCPUs);
    } else if (lastCPU < overflow) {
        return ((minNrLines + 1) * nrBufCPUs);
    } else {
        return (minNrLines * nrBufCPUs + overflow - firstCPU);
    }
}


static inline void
PxGetSBTorder(int root, int orderLength, int *order, int *myIndex)
{
    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    if (root >= _nrCPUs) {
        return;
    }
    *myIndex    = 0;
    order[0]    = root;
    int rootTmp = root;
    int index   = orderLength / 2;

    for (int i=0; i<_logCPUs; i++) {
        if (rootTmp < index) {
            for (int j=0; j<index; j++) {
                int cpu = index+j + root-rootTmp;
                if (cpu == _myCPU) {
                    *myIndex = index+j;
                }
                order[index+j] = cpu;
            }
        } else {
            for (int j=0; j<index; j++) {
                int cpu = j + root-rootTmp;
                if (cpu == _myCPU) {
                    *myIndex = index+j;
                }
                order[index+j] = cpu;
            }
            rootTmp -= index;
        }
        index /= 2;
    }
}


/**********************************************************************
 *** Communication initialization                                   ***
 **********************************************************************/

inline void
PxInitCommunication()
{
    char        dummy;
    MPI_Status  stat;

    /*** Send up and receive down ***/

    if (PxMyUpCPU() != -1) {
        MPI_Send(&dummy, 1, MPI_BYTE,
                        PxMyUpCPU(), INIT_TAG, MPI_COMM_WORLD);
    }
    if (PxMyDownCPU() != -1) {
        MPI_Recv(&dummy, 1, MPI_BYTE,
                        PxMyDownCPU(), INIT_TAG, MPI_COMM_WORLD, &stat);

    /*** Send down and receive up ***/

        MPI_Send(&dummy, 1, MPI_BYTE,
                        PxMyDownCPU(), INIT_TAG, MPI_COMM_WORLD);
    }
    if (PxMyUpCPU() != -1) {
        MPI_Recv(&dummy, 1, MPI_BYTE,
                        PxMyUpCPU(), INIT_TAG, MPI_COMM_WORLD, &stat);
    }
}



/**********************************************************************
 *** Broadcast Value                                                ***
 **********************************************************************/

template<class ArithT>
static inline ArithT
PxBcastValueOFT(ArithT val, int root)
{
    /*** Broadcast value using One-level Flat Tree ***/

    MPI_Datatype elem;
    MPI_Status stat;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    if (_myCPU == root) {               // send value to all other CPUs
        for (int partner=0; partner<_nrCPUs; partner++) {
            if (partner != _myCPU) {
                MPI_Send(&val,
                         1, elem, partner, BCAST_TAG, MPI_COMM_WORLD);
            }
        }
    } else {                            // receive value from root
        MPI_Recv(&val, 1, elem, root, BCAST_TAG, MPI_COMM_WORLD, &stat);
    }
    MPI_Type_free(&elem);
    return val;
}


template<class ArithT>
static inline ArithT
PxBcastValueSBT(ArithT val, int root)
{
    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    int *order = new int[_maxCPUs], myIndex;
    PxGetSBTorder(root, _maxCPUs, order, &myIndex);

    /*** Broadcast value using Spanning Binomial Tree ***/

    MPI_Datatype elem;
    MPI_Status stat;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    int mask = 1 << _logCPUs-1;
    for (int i=0; i<_logCPUs; i++) {
        int partnerIndex = myIndex ^ mask;
        int partner = order[partnerIndex];
        if ((myIndex % mask == 0) && (partner < _nrCPUs)) {
            if (myIndex < partnerIndex) {
                MPI_Send(&val, 1, elem,
                         partner, BCAST_TAG, MPI_COMM_WORLD);
            } else {
                MPI_Recv(&val, 1, elem,
                         partner, BCAST_TAG, MPI_COMM_WORLD, &stat);
            }
        }
        mask >>= 1;
    }
    MPI_Type_free(&elem);
    delete order;
    return val;
}


template<class ArithT>
static inline ArithT
PxBcastValueMPI(ArithT val, int root)
{
    /*** MPI-Bcast. May be faster than OFT/SBT broadcast provided   ***/
    /*** above as specific communication subsystem capabilities may ***/
    /*** be incorporated in the implementation.                     ***/

    MPI_Datatype elem;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);
    MPI_Bcast(&val, 1, elem, root, MPI_COMM_WORLD);
    MPI_Type_free(&elem);
    return val;
}


template<class ArithT>
inline ArithT
PxBcastValue(ArithT val, int root=0, int dtype=PX_MPI)
{
    /*** root  : rank of broadcast root node                        ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MPI)      ***/
    
    switch (dtype) {
        case PX_OFT : return PxBcastValueOFT(val, root); break;
        case PX_SBT : return PxBcastValueSBT(val, root); break;
        case PX_MPI :
        default         : return PxBcastValueMPI(val, root);
    }
}


/**********************************************************************
 *** Reduce Value to Root                                           ***
 **********************************************************************/

template<class ArithT, class RedOpT>
static inline void
PxReduce(ArithT *in, ArithT *inout, int *len, MPI_Datatype *dptr)
{
    /*** This is the user defined function as MPI likes to have it ***/

    RedOpT op;
    op.DoIt(*inout, *in);
}


template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToRootOFT(ArithT val, RedOpT redOp, int root)
{
    /*** Reduce value to root using One-level Flat Tree ***/

    ArithT result = ArithT(val);
    MPI_Datatype elem;
    MPI_Status stat;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    if (_myCPU == root) {           // receive value from all other CPUs
        for (int partner=0; partner<_nrCPUs; partner++) {
            if (partner != _myCPU) {
                ArithT recvVal;
                MPI_Recv(&recvVal, 1, elem,
                         partner, RDUCE_TAG, MPI_COMM_WORLD, &stat);
                redOp.DoIt(result, recvVal);
            }
        }
    } else {                        // send value to root
        MPI_Send(&val, 1, elem, root, RDUCE_TAG, MPI_COMM_WORLD);
    }
    MPI_Type_free(&elem);
    return result;
}


template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToRootSBT(ArithT val, RedOpT redOp, int root)
{
    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    int *order = new int[_maxCPUs], myIndex;
    PxGetSBTorder(root, _maxCPUs, order, &myIndex);

    /*** Reduce value to root using Spanning Binomial Tree ***/

    ArithT result = ArithT(val);
    MPI_Datatype elem;
    MPI_Status stat;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    int mask = 1;
    for (int i=0; i<_logCPUs; i++) {
        int partnerIndex = myIndex ^ mask;
        int partner = order[partnerIndex];
        if ((myIndex % mask == 0) && (partner < _nrCPUs)) {
            if (myIndex > partnerIndex) {
                MPI_Send(&result, 1, elem,
                         partner, RDUCE_TAG, MPI_COMM_WORLD);
            } else {
                ArithT recvVal;
                MPI_Recv(&recvVal, 1, elem,
                         partner, RDUCE_TAG, MPI_COMM_WORLD, &stat);
                redOp.DoIt(result, recvVal);
            }
        }
        mask <<= 1;
    }
    MPI_Type_free(&elem);
    delete order;
    return result;
}


template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToRootMPI(ArithT val, RedOpT redOp, int root)
{
    /*** MPI-based global reduction NOT using any of the predefined ***/
    /*** reduce operations supported by MPI. At all times, a handle ***/
    /*** to a correctly instantiated (Horus) reduction operation is ***/
    /*** passed as a parameter to MPI_Reduce(...).                  ***/

    int commute = 1;
    ArithT result = RedOpT::NeutralElement();
    MPI_Op operation;
    MPI_Datatype elem;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    /*** Make instantiation of PxReduce, bind to MPIop handle & run ***/

    void (*reduceOp)(ArithT*, ArithT*, int*, MPI_Datatype*) =
                                            PxReduce<ArithT, RedOpT>;
    MPI_Op_create((void (*)(void*, void*, int*, MPI_Datatype*))reduceOp,
                                            commute, &operation);
    MPI_Reduce(&val, &result, 1, elem, operation, root, MPI_COMM_WORLD);
    MPI_Op_free(&operation);
    MPI_Type_free(&elem);
    return result;
}


template<class ArithT, class RedOpT>
inline ArithT
PxReduceValueToRoot(ArithT val, RedOpT redOp, int root=0, int dtype=PX_MPI)
{
    /*** redOp : Horus reduction operation                          ***/
    /*** root  : rank of reduction root                             ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MPI)      ***/
    
    switch (dtype) {
        case PX_OFT : return PxReduceValueToRootOFT(val, redOp, root);
                                          break;
        case PX_SBT : return PxReduceValueToRootSBT(val, redOp, root);
                                          break;
        case PX_MPI :
        default         : return PxReduceValueToRootMPI(val, redOp, root);
    }
}


/**********************************************************************
 *** Reduce Value to All                                            ***
 **********************************************************************/

template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToAllOFT(ArithT val, RedOpT redOp)
{
    /*** All-to-one (root) and one-to-all (broadcast) using OFT ***/

    return PxBcastValueOFT(PxReduceValueToRootOFT(val, redOp, 0), 0);
}


template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToAllSBT(ArithT val, RedOpT redOp)
{
    /*** All-to-one (root) and one-to-all (broadcast) using SBT ***/

    return PxBcastValueSBT(PxReduceValueToRootSBT(val, redOp, 0), 0);
}


template<class ArithT, class RedOpT>
static inline ArithT
PxReduceValueToAllMPI(ArithT val, RedOpT redOp)
{
    /*** MPI-based global Allreduce NOT using any of the predefined ***/
    /*** reduce operations supported by MPI. At all times, a handle ***/
    /*** to a correctly instantiated (Horus) reduction operation is ***/
    /*** passed as a parameter to MPI_Allreduce(...).               ***/

    int commute = 1;
    ArithT result = RedOpT::NeutralElement();
    MPI_Op operation;
    MPI_Datatype elem;
    MPI_Type_contiguous(sizeof(ArithT), MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    /*** Make instantiation of PxReduce, bind to MPIop handle & run ***/

    void (*reduceOp)(ArithT*, ArithT*, int*, MPI_Datatype*) =
                                            PxReduce<ArithT, RedOpT>;
    MPI_Op_create((void (*)(void*, void*, int*, MPI_Datatype*))reduceOp,
                                            commute, &operation);
    MPI_Allreduce(&val, &result, 1, elem, operation, MPI_COMM_WORLD);
    MPI_Op_free(&operation);
    MPI_Type_free(&elem);
    return result;
}


template<class ArithT, class RedOpT>
inline ArithT
PxReduceValueToAll(ArithT val, RedOpT redOp, int dtype=PX_MPI)
{
    /*** redOp : Horus reduction operation                          ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MPI)      ***/
    
    switch (dtype) {
        case PX_OFT : return PxReduceValueToAllOFT(val, redOp); break;
        case PX_SBT : return PxReduceValueToAllSBT(val, redOp); break;
        case PX_MPI :
        default         : return PxReduceValueToAllMPI(val, redOp);
    }
}


/**********************************************************************
 *** Scatter Array Data                                             ***
 **********************************************************************/

template<class ArrayT>
static inline void
PxScatterArrayOFT(ArrayT* glob, ArrayT* loc, int root, bool bdata)
{
    /*** Scatter array data using One-level Flat Tree ***/

    int eSize = ArrayT::ElemSize();
    int sSize = eSize*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);

    if (_myCPU == root) {           // send array data to all other CPUs
        for (int partner=0; partner<_nrCPUs; partner++) {
            if (partner != _myCPU) {
                int xSize  = PxLclWidth(partner);
                int ySize  = PxLclHeight(partner);
                int zSize  = PxLclDepth(partner);
                int offset = PxLclStart(partner, ArrayBW(glob),
                             ArrayBH(glob), ArrayBD(glob));
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                    offset -=
                    (ArrayW(glob)*ArrayH(glob)*ArrayBD(glob) +
                     ArrayW(glob)*ArrayBH(glob) + ArrayBW(glob));
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize,
                                    ArrayW(glob), elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                    ArrayW(glob)*ArrayH(glob)*sSize,
                                    blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Send(ArrayPB(glob) + offset*eSize, 1, blk3d,
                                    partner, SCAT_TAG, MPI_COMM_WORLD);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            }
        }
        PxLclArrayCopy(loc, FROM, glob, root, bdata);

    } else {                        // receive array data from root
        int xSize = PxLclWidth(_myCPU);
        int ySize = PxLclHeight(_myCPU);
        int zSize = PxLclDepth(_myCPU);
        if (bdata) {
            xSize += 2*ArrayBW(glob);
            ySize += 2*ArrayBH(glob);
            zSize += 2*ArrayBD(glob);
        }
        if (xSize != 0 && ySize != 0 && zSize != 0) {
            MPI_Type_vector(ySize, xSize, ArrayW(loc), elem, &blk2d);
            MPI_Type_hvector(zSize, 1,
                    ArrayW(loc)*ArrayH(loc)*sSize, blk2d, &blk3d);
            MPI_Type_commit(&blk3d);
            MPI_Recv((bdata) ? ArrayPB(loc) : ArrayCPB(loc), 1,
                    blk3d, root, SCAT_TAG, MPI_COMM_WORLD, &stat);
            MPI_Type_free(&blk3d);
            MPI_Type_free(&blk2d);
        }
    }
    MPI_Type_free(&elem);
}


template<class ArrayT>
static inline void
PxScatterArraySBT(ArrayT* glob, ArrayT* loc, int root, bool bdata)
{
    /*** SBT (Spanning Binomial Tree) or hypercube scatter. Its  ***/
    /*** use is restricted slightly for implementation reasons.  ***/
    /*** NOTE: should work correctly, except:                    ***/
    /***                                                         ***/
    /*** - IF ("zCPUs" == 1) AND ("yCPUs" > 1) THEN              ***/
    /***   "xCPUs" MUST be a power of 2 (so: 1, 2, 4, 8, etc...) ***/
    /*** - IF ("zCPUs" > 1) THEN                                 ***/
    /***   "xCPUs" AND "yCPUs" both MUST be a power of 2         ***/

    int     *order = new int[_maxCPUs], myIndex, bufW, bufH, bufD;
    ArrayT  *buf;

    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    PxGetSBTorder(root, _maxCPUs, order, &myIndex);

    /*** Create temporary array buffer (including borders) ***/

    if (_myCPU == root) {
        bufW = ArrayW(glob);
        bufH = ArrayH(glob);
        bufD = ArrayD(glob);
        buf  = glob;
    } else {
        bufW = PxLineRange(1, myIndex, order);
        bufH = PxLineRange(2, myIndex, order);
        bufD = PxLineRange(3, myIndex, order);
        buf  = PxArrayCreate<ArrayT>(bufW, bufH, bufD, ArrayBW(glob),
                                     ArrayBH(glob), ArrayBD(glob));
        bufW += 2*ArrayBW(glob);
        bufH += 2*ArrayBH(glob);
        bufD += 2*ArrayBD(glob);
    }

    /*** Scatter array data using Spanning Binomial Tree ***/

    int eSize = ArrayT::ElemSize();
    int sSize = eSize*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);

    int mask = 1 << (_logCPUs-1);
    for (int i=0; i<_logCPUs; i++) {
        int partnerIndex = myIndex ^ mask;
        int partner = order[partnerIndex];
        if ((myIndex % mask == 0) && (partner < _nrCPUs)) {
            if (myIndex < partnerIndex) {   // send data to SBT child
                int xSize  = PxLineRange(1, partnerIndex, order);
                int ySize  = PxLineRange(2, partnerIndex, order);
                int zSize  = PxLineRange(3, partnerIndex, order);
                int offset = PxLclStart(partner, ArrayBW(glob),
                                    ArrayBH(glob), ArrayBD(glob));
                if (_myCPU != root) {
                    offset -= (PxLclStart(_myCPU, ArrayBW(glob),
                               ArrayBH(glob), ArrayBD(glob)) -
                              (bufW*bufH*ArrayBD(glob) +
                               bufW*ArrayBH(glob) + ArrayBW(glob)));
                }
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                    offset -= (bufW*bufH*ArrayBD(glob) +
                               bufW*ArrayBH(glob) + ArrayBW(glob));
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize, bufW, elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                     bufW*bufH*sSize, blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Send(ArrayPB(buf) + offset*eSize, 1, blk3d,
                             partner, SCAT_TAG, MPI_COMM_WORLD);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            } else {                    // receive data from SBT parent
                int xSize = PxLineRange(1, myIndex, order);
                int ySize = PxLineRange(2, myIndex, order);
                int zSize = PxLineRange(3, myIndex, order);
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize, bufW, elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                     bufW*bufH*sSize, blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Recv((bdata) ? ArrayPB(buf) : ArrayCPB(buf),
                            1, blk3d, partner, SCAT_TAG, MPI_COMM_WORLD,
                            &stat);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            }
        }
        mask >>= 1;
    }
    MPI_Type_free(&elem);
    PxLclArrayCopy(loc, FROM, buf, root, bdata);
    if (_myCPU != root) {
        delete buf;
    }
    delete order;
}


template<class ArrayT>
static inline void
PxScatterArrayMPI(ArrayT* glob, ArrayT* loc, int root, int bdata=true)
{
    /*** MPI-Scatter. May be faster than OFT/SBT scatter provided   ***/
    /*** above as specific communication subsystem capabilities may ***/
    /*** be incorporated in the implementation.                     ***/
    /***                                                            ***/
    /*** NOTE: This function can only be used if:                   ***/
    /***                                                            ***/
    /*** Borders must be sent as well OR all border sizes are 0 AND ***/
    /*** xcpus>=1 AND ArrayH(glob)=1 AND ArrayD(glob)=1, OR     ***/
    /*** xcpus=1 AND ycpus>=1 AND ArrayD(glob)=1, OR              ***/
    /*** xcpus=1 AND ycpus=1 AND zcpus>=1.                          ***/
    /***                                                            ***/
    /*** This is due to the lack of possibility to define an array  ***/
    /*** of derived datatypes at the root. In principle, we could   ***/
    /*** temporarily re-arrange the layout of the array to still be ***/
    /*** able to use the standard MPI functionality. However, this  ***/
    /*** requires a costly copy operation which is not implemented. ***/

    int *counts, *disps;

    /*** Determine size and position of each local array ***/

    if (_myCPU == root) {
        counts = new int[_nrCPUs];
        disps  = new int[_nrCPUs];
        int slice = 2*((ArrayW(glob)*ArrayH(glob)*ArrayBD(glob)) +
                       (ArrayW(glob)*ArrayBH(glob)));
        for (int i=0; i<_nrCPUs; i++) {
            counts[i] = (PxLclWidth(i)  + 2*ArrayBW(glob)) *
                        (PxLclHeight(i) + 2*ArrayBH(glob)) *
                        (PxLclDepth(i)  + 2*ArrayBD(glob));
            disps[i] = (i==0) ? 0 : disps[i-1] + counts[i-1] - slice; 
        }
    }

    /*** Scatter array data ***/

    int sSize = ArrayT::ElemSize()*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_commit(&elem);
    MPI_Scatterv(ArrayPB(glob), counts, disps, elem, ArrayPB(loc),
                 ArrayW(loc)*ArrayH(loc)*ArrayD(loc), elem, root,
                 MPI_COMM_WORLD);
    MPI_Type_free(&elem);

    if (_myCPU == root) {
        delete counts;
        delete disps;
    }
}


template<class ArrayT>
inline void
PxScatterArray(ArrayT* glob, ArrayT* loc,
               int root=0, int dtype=PX_SBT, bool bdata=false)
{
    /*** root  : rank of sending node                               ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MPI)      ***/
    /*** bdata : indicates whether border data is scattered as well ***/
    
if (dp) PX_COUT << "Scatter: " << "NrCPUs = " << PxNrCPUs() << PX_ENDL;

    switch (dtype) {
        case PX_MPI : PxScatterArrayMPI(glob, loc, root, true ); break;
        case PX_OFT : PxScatterArrayOFT(glob, loc, root, bdata); break;
        case PX_SBT :
        default         : PxScatterArraySBT(glob, loc, root, bdata);
    }
}


/**********************************************************************
 *** Gather Array Data                                              ***
 **********************************************************************/

template<class ArrayT>
static inline void
PxGatherArrayOFT(ArrayT* glob, ArrayT* loc, int root, bool bdata)
{
    /*** Gather array data using One-level Flat Tree ***/

    int eSize = ArrayT::ElemSize();
    int sSize = eSize*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);

    if (_myCPU == root) {   // receive array data from all other CPUs
        PxLclArrayCopy(loc, TO, glob, root, bdata);
        for (int partner=0; partner<_nrCPUs; partner++) {
            if (partner != _myCPU) {
                int xSize  = PxLclWidth(partner);
                int ySize  = PxLclHeight(partner);
                int zSize  = PxLclDepth(partner);
                int offset = PxLclStart(partner,ArrayBW(glob),
                             ArrayBH(glob), ArrayBD(glob));
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                    offset -=
                    (ArrayW(glob)*ArrayH(glob)*ArrayBD(glob) +
                     ArrayW(glob)*ArrayBH(glob) + ArrayBW(glob));
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize,
                                    ArrayW(glob), elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                    ArrayW(glob)*ArrayH(glob)*sSize,
                                    blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Recv(ArrayPB(glob) + offset*eSize, 1, blk3d,
                             partner, GATH_TAG, MPI_COMM_WORLD, &stat);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            }
        }
    } else {                // send array data to root
        int xSize = PxLclWidth(_myCPU);
        int ySize = PxLclHeight(_myCPU);
        int zSize = PxLclDepth(_myCPU);
        if (bdata) {
            xSize += 2*ArrayBW(glob);
            ySize += 2*ArrayBH(glob);
            zSize += 2*ArrayBD(glob);
        }
        if (xSize != 0 && ySize != 0 && zSize != 0) {
            MPI_Type_vector(ySize, xSize, ArrayW(loc), elem, &blk2d);
            MPI_Type_hvector(zSize, 1,
                    ArrayW(loc)*ArrayH(loc)*sSize, blk2d, &blk3d);
            MPI_Type_commit(&blk3d);
            MPI_Send((bdata) ? ArrayPB(loc) : ArrayCPB(loc),
                     1, blk3d, root, GATH_TAG, MPI_COMM_WORLD);
            MPI_Type_free(&blk3d);
            MPI_Type_free(&blk2d);
        }
    }
    MPI_Type_free(&elem);
}


template<class ArrayT>
static inline void
PxGatherArraySBT(ArrayT* glob, ArrayT* loc, int root, bool bdata)
{
    /*** SBT (Spanning Binomial Tree) or hypercube gather. Its   ***/
    /*** use is restricted slightly for implementation reasons.  ***/
    /*** NOTE: should work correctly, except:                    ***/
    /***                                                         ***/
    /*** - IF ("zCPUs" == 1) AND ("yCPUs" > 1) THEN              ***/
    /***   "xCPUs" MUST be a power of 2 (so: 1, 2, 4, 8, etc...) ***/
    /*** - IF ("zCPUs" > 1) THEN                                 ***/
    /***   "xCPUs" AND "yCPUs" both MUST be a power of 2         ***/

    int     *order = new int[_maxCPUs], myIndex, bufW, bufH, bufD;
    ArrayT  *buf;

    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    PxGetSBTorder(root, _maxCPUs, order, &myIndex);

    /*** Create temporary array buffer (including borders) ***/

    if (_myCPU == root) {
        bufW = ArrayW(glob);
        bufH = ArrayH(glob);
        bufD = ArrayD(glob);
        buf  = glob;
    } else {
        bufW = PxLineRange(1, myIndex, order);
        bufH = PxLineRange(2, myIndex, order);
        bufD = PxLineRange(3, myIndex, order);
        buf  = PxArrayCreate<ArrayT>(bufW, bufH, bufD, ArrayBW(glob),
                                     ArrayBH(glob), ArrayBD(glob));
        bufW += 2*ArrayBW(glob);
        bufH += 2*ArrayBH(glob);
        bufD += 2*ArrayBD(glob);
    }
    PxLclArrayCopy(loc, TO, buf, root, bdata);

    /*** Gather array data using Spanning Binomial Tree ***/

    int eSize = ArrayT::ElemSize();
    int sSize = eSize*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);

    int mask = 1;
    for (int i=0; i<_logCPUs; i++) {
        int partnerIndex = myIndex ^ mask;
        int partner = order[partnerIndex];
        if ((myIndex % mask == 0) && (partner < _nrCPUs)) {
            if (myIndex > partnerIndex) {   // send data to SBT parent
                int xSize = PxLineRange(1, myIndex, order);
                int ySize = PxLineRange(2, myIndex, order);
                int zSize = PxLineRange(3, myIndex, order);
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize, bufW, elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                     bufW*bufH*sSize, blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Send((bdata) ? ArrayPB(buf) : ArrayCPB(buf),
                           1, blk3d, partner, GATH_TAG, MPI_COMM_WORLD);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            } else {                    // receive data from SBT child
                int xSize  = PxLineRange(1, partnerIndex, order);
                int ySize  = PxLineRange(2, partnerIndex, order);
                int zSize  = PxLineRange(3, partnerIndex, order);
                int offset = PxLclStart(partner, ArrayBW(glob),
                                    ArrayBH(glob), ArrayBD(glob));
                if (_myCPU != root) {
                    offset -= (PxLclStart(_myCPU, ArrayBW(glob),
                               ArrayBH(glob), ArrayBD(glob)) -
                              (bufW*bufH*ArrayBD(glob) +
                               bufW*ArrayBH(glob) + ArrayBW(glob)));
                }
                if (bdata) {
                    xSize += 2*ArrayBW(glob);
                    ySize += 2*ArrayBH(glob);
                    zSize += 2*ArrayBD(glob);
                    offset -= (bufW*bufH*ArrayBD(glob) +
                               bufW*ArrayBH(glob) + ArrayBW(glob));
                }
                if (xSize != 0 && ySize != 0 && zSize != 0) {
                    MPI_Type_vector(ySize, xSize, bufW, elem, &blk2d);
                    MPI_Type_hvector(zSize, 1,
                                     bufW*bufH*sSize, blk2d, &blk3d);
                    MPI_Type_commit(&blk3d);
                    MPI_Recv(ArrayPB(buf) + offset*eSize, 1, blk3d,
                             partner, GATH_TAG, MPI_COMM_WORLD, &stat);
                    MPI_Type_free(&blk3d);
                    MPI_Type_free(&blk2d);
                }
            }
        }
        mask <<= 1;
    }
    MPI_Type_free(&elem);
    if (_myCPU != root) {
        delete buf;
    }
    delete order;
}


template<class ArrayT>
static inline void
PxGatherArrayMPI(ArrayT* glob, ArrayT* loc, int root, bool bdata=true)
{
    /*** MPI-Gather. May be faster than OFT/SBT scatter provided    ***/
    /*** above as specific communication subsystem capabilities may ***/
    /*** be incorporated in the implementation.                     ***/
    /***                                                            ***/
    /*** NOTE: This function can only be used if:                   ***/
    /***                                                            ***/
    /*** Borders must be sent as well OR all border sizes are 0 AND ***/
    /*** xcpus>=1 AND ArrayH(glob)=1 AND ArrayD(glob)=1, OR     ***/
    /*** xcpus=1 AND ycpus>=1 AND ArrayD(glob)=1, OR              ***/
    /*** xcpus=1 AND ycpus=1 AND zcpus>=1.                          ***/
    /***                                                            ***/
    /*** This is due to the lack of possibility to define an array  ***/
    /*** of derived datatypes at the root. In principle, we could   ***/
    /*** temporarily re-arrange the layout of the array to still be ***/
    /*** able to use the standard MPI functionality. However, this  ***/
    /*** requires a costly copy operation which is not implemented. ***/

    int *counts, *disps;

    /*** Determine size and position of each local array ***/

    if (_myCPU == root) {
        counts = new int[_nrCPUs];
        disps  = new int[_nrCPUs];
        int slice = 2*((ArrayW(glob)*ArrayH(glob)*ArrayBD(glob)) +
                       (ArrayW(glob)*ArrayBH(glob)));
        for (int i=0; i<_nrCPUs; i++) {
            counts[i] = (PxLclWidth(i)  + 2*ArrayBW(glob)) *
                        (PxLclHeight(i) + 2*ArrayBH(glob)) *
                        (PxLclDepth(i)  + 2*ArrayBD(glob));
            disps[i] = (i==0) ? 0 : disps[i-1] + counts[i-1] - slice; 
        }
    }

    /*** Scatter array data ***/

    int sSize = ArrayT::ElemSize()*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_commit(&elem);
    MPI_Gatherv(ArrayPB(loc),
                ArrayW(loc)*ArrayH(loc)*ArrayD(loc), elem,
                ArrayPB(glob), counts, disps, elem, root,
                MPI_COMM_WORLD);
    MPI_Type_free(&elem);

    if (_myCPU == root) {
        delete counts;
        delete disps;
    }
}


template<class ArrayT>
inline void
PxGatherArray(ArrayT* glob, ArrayT* loc,
              int root=0, int dtype=PX_SBT, bool bdata=false)
{
    /*** root  : rank of receiving node                             ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MP)       ***/
    /*** bdata : indicates whether border data is gathered as well  ***/
    
if (dp) PX_COUT << "Gather..." << PX_ENDL;

    switch (dtype) {
        case PX_MPI : PxGatherArrayMPI(glob, loc, root, true ); break;
        case PX_OFT : PxGatherArrayOFT(glob, loc, root, bdata); break;
        case PX_SBT :
        default         : PxGatherArraySBT(glob, loc, root, bdata);
    }
}


/**********************************************************************
 *** Broadcast Array Data                                           ***
 **********************************************************************/

template<class ArrayT>
static inline void
PxBcastArrayOFT(ArrayT* a, int root, bool bdata)
{
    /*** Broadcast array data using One-level Flat Tree ***/

    int sSize = ArrayT::ElemSize()*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    if (!bdata) {
        MPI_Type_vector(ArrayCH(a),
                        ArrayCW(a), ArrayW(a), elem, &blk2d);
        MPI_Type_hvector(ArrayCD(a), 1,
                        ArrayW(a)*ArrayH(a)*sSize, blk2d, &blk3d);
        MPI_Type_commit(&blk3d);
    }
    if (_myCPU == root) {       // send array data to all other CPUs
        for (int partner=0; partner<_nrCPUs; partner++) {
            if (partner != _myCPU) {
                if (!bdata) {
                    MPI_Send(ArrayCPB(a), 1,
                             blk3d, partner, BCAST_TAG, MPI_COMM_WORLD);
                } else {
                    MPI_Send(ArrayPB(a),
                             ArrayW(a)*ArrayH(a)*ArrayD(a),
                             elem, partner, BCAST_TAG, MPI_COMM_WORLD);
                }
            }
        }
    } else {                    // receive array data from root
        if (!bdata) {
            MPI_Recv(ArrayCPB(a), 1, blk3d, root,
                     BCAST_TAG, MPI_COMM_WORLD, &stat);
        } else {
            MPI_Recv(ArrayPB(a), ArrayW(a)*ArrayH(a)*ArrayD(a),
                     elem, root, BCAST_TAG, MPI_COMM_WORLD, &stat);
        }
    }
    if (!bdata) {
        MPI_Type_free(&blk3d);
        MPI_Type_free(&blk2d);
    }
    MPI_Type_free(&elem);
}


template<class ArrayT>
static inline void
PxBcastArraySBT(ArrayT* a, int root, bool bdata)
{
    /*** Determine ordering of CPUs in Spanning Binomial Tree ***/

    int *order = new int[_maxCPUs], myIndex;
    PxGetSBTorder(root, _maxCPUs, order, &myIndex);

    /*** Broadcast array data using Spanning Binomial Tree ***/

    int sSize = ArrayT::ElemSize()*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    if (!bdata) {
        MPI_Type_vector(ArrayCH(a),
                        ArrayCW(a), ArrayW(a), elem, &blk2d);
        MPI_Type_hvector(ArrayCD(a), 1,
                        ArrayW(a)*ArrayH(a)*sSize, blk2d, &blk3d);
        MPI_Type_commit(&blk3d);
    }
    int mask = 1 << (_logCPUs-1);
    for (int i=0; i<_logCPUs; i++) {
        int partnerIndex = myIndex ^ mask;
        int partner = order[partnerIndex];
        if ((myIndex % mask == 0) && (partner < _nrCPUs)) {
            if (myIndex < partnerIndex) {   // send data to SBT child
                if (!bdata) {
                    MPI_Send(ArrayCPB(a), 1, blk3d,
                             partner, BCAST_TAG, MPI_COMM_WORLD);
                } else {
                    MPI_Send(ArrayPB(a),
                             ArrayW(a)*ArrayH(a)*ArrayD(a), elem,
                             partner, BCAST_TAG, MPI_COMM_WORLD);
                }
            } else {                    // receive data from SBT parent
                if (!bdata) {
                    MPI_Recv(ArrayCPB(a), 1, blk3d, partner,
                             BCAST_TAG, MPI_COMM_WORLD, &stat);
                } else {
                    MPI_Recv(ArrayPB(a),
                             ArrayW(a)*ArrayH(a)*ArrayD(a), elem,
                             partner, BCAST_TAG, MPI_COMM_WORLD, &stat);
                }
            }
        }
        mask >>= 1;
    }
    if (!bdata) {
        MPI_Type_free(&blk3d);
        MPI_Type_free(&blk2d);
    }
    MPI_Type_free(&elem);
    delete order;
}


template<class ArrayT>
static inline void
PxBcastArrayMPI(ArrayT* a, int root, bool bdata)
{
    /*** MPI-Bcast. May be faster than OFT/SBT broadcast provided ***/
    /*** above, as specific communication subsystem capabilities  ***/
    /*** may be incorporated in the implementation.               ***/

    int sSize = ArrayT::ElemSize()*sizeof(typename ArrayT::StorType);
    MPI_Datatype elem, blk2d, blk3d;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_commit(&elem);

    if (!bdata) {
        MPI_Type_vector(ArrayCH(a),
                        ArrayCW(a), ArrayW(a), elem, &blk2d);
        MPI_Type_hvector(ArrayCD(a), 1,
                         ArrayW(a)*ArrayH(a)*sSize, blk2d, &blk3d);
        MPI_Type_commit(&blk3d);
        MPI_Bcast(ArrayCPB(a), 1, blk3d, root, MPI_COMM_WORLD);
        MPI_Type_free(&blk3d);
        MPI_Type_free(&blk2d);
    } else {
        MPI_Bcast(ArrayPB(a), ArrayW(a)*ArrayH(a)*ArrayD(a),
                  elem, root, MPI_COMM_WORLD);
    }
    MPI_Type_free(&elem);
}


template<class ArrayT>
inline void
PxBcastArray(ArrayT* a, int root=0, int dtype=PX_MPI, bool bdata=false)
{
    /*** root  : rank of broadcast root node                        ***/
    /*** dtype : distribution type (PX_OFT, PX_SBT, or PX_MPI)      ***/
    /*** bdata : indicates whether border data is broadcast as well ***/
    
if (dp) PX_COUT << "Bcast: " << "NrCPUs = " << PxNrCPUs() << PX_ENDL;

    switch (dtype) {
        case PX_OFT : PxBcastArrayOFT(a, root, bdata); break;
        case PX_SBT : PxBcastArraySBT(a, root, bdata); break;
        case PX_MPI :
        default         : PxBcastArrayMPI(a, root, bdata);
    }
}


/**********************************************************************
 *** Redistribute Array Data                                        ***
 **********************************************************************/

template<class ArrayT>
inline void
PxRedistArray(ArrayT** a,
              int xcpus, int ycpus, int zcpus, bool bdata=false)
{
    /*** Complete redistribution of array data (incl. potential    ***/
    /*** borders) using the Direct Communication schema. This code ***/
    /*** is based on: 'Parallel Matrix Transpose Algorithms on     ***/
    /*** Distributed Memory Concurrent Computers' (Choi, Dongarra, ***/
    /*** Walker) Mathematical Sciences Section, Oak Ridge National ***/
    /*** Lab., and Dept. Computer Science, Univ. of Tennessee at   ***/
    /*** Knoxville, LAPACK Working Note 65 (CS-93-215), Oct 1993.  ***/

    int bw = ArrayBW(*a), bh = ArrayBH(*a), bd = ArrayBD(*a);

    /*** Save old local image characteristics ***/

    int* oldStartX = new int[_nrCPUs];
    int* oldStartY = new int[_nrCPUs];
    int* oldStartZ = new int[_nrCPUs];
    int* oldEndX = new int[_nrCPUs];
    int* oldEndY = new int[_nrCPUs];
    int* oldEndZ = new int[_nrCPUs];

    int i;
    for (i=0; i<_nrCPUs; i++) {
        oldStartX[i] = PxLclStartX(i);
        oldStartY[i] = PxLclStartY(i);
        oldStartZ[i] = PxLclStartZ(i);
        oldEndX[i]   = oldStartX[i] + PxLclWidth(i)  - 1;
        oldEndY[i]   = oldStartY[i] + PxLclHeight(i) - 1;
        oldEndZ[i]   = oldStartZ[i] + PxLclDepth(i)  - 1;
        if (bdata) {
            oldEndX[i] += 2*bw;
            oldEndY[i] += 2*bh;
            oldEndZ[i] += 2*bd;
        }
    }

    /*** Repartition according to new CPU grid ***/

    PxRePartition(xcpus, ycpus, zcpus);

    /*** Obtain new local image characteristics ***/

    int* newStartX = new int[_nrCPUs];
    int* newStartY = new int[_nrCPUs];
    int* newStartZ = new int[_nrCPUs];
    int* newEndX = new int[_nrCPUs];
    int* newEndY = new int[_nrCPUs];
    int* newEndZ = new int[_nrCPUs];

    for (i=0; i<_nrCPUs; i++) {
        newStartX[i] = PxLclStartX(i);
        newStartY[i] = PxLclStartY(i);
        newStartZ[i] = PxLclStartZ(i);
        newEndX[i]   = newStartX[i] + PxLclWidth(i)  - 1;
        newEndY[i]   = newStartY[i] + PxLclHeight(i) - 1;
        newEndZ[i]   = newStartZ[i] + PxLclDepth(i)  - 1;
        if (bdata) {
            newEndX[i] += 2*bw;
            newEndY[i] += 2*bh;
            newEndZ[i] += 2*bd;
        }
    }
    int newW = PxLclWidth(_myCPU);
    int newH = PxLclHeight(_myCPU);
    int newD = PxLclDepth(_myCPU);

    /*** Create new local array and temporary MPI buffer ***/

    ArrayT* nA = PxArrayCreate<ArrayT>(newW, newH, newD, bw, bh, bd);
    typedef typename ArrayT::StorType storT;
    int eSize = ArrayT::ElemSize();
    int sSize = eSize*sizeof(storT);
    int mpibuffsize = sizeof(MPI_BYTE) * sSize *
                      ArrayW(nA) * ArrayH(nA) * ArrayD(nA);
//  int mpibuffsize = sizeof(MPI_BYTE) * sSize * newW * newH * newD;
    char *mpibuff = new char[mpibuffsize];
    MPI_Buffer_attach(mpibuff, mpibuffsize);

    /*** Send to/receive from all other CPUs (direct communication) ***/

    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);

    for (i=0; i<_nrCPUs; i++) {

        int partner = (_nrCPUs - _myCPU + i) % _nrCPUs;

        /*** Send data to partner ***/

        int maxStartX = (oldStartX[_myCPU] > newStartX[partner]) ?
                         oldStartX[_myCPU] : newStartX[partner];
        int maxStartY = (oldStartY[_myCPU] > newStartY[partner]) ?
                         oldStartY[_myCPU] : newStartY[partner];
        int maxStartZ = (oldStartZ[_myCPU] > newStartZ[partner]) ?
                         oldStartZ[_myCPU] : newStartZ[partner];
        int minEndX   = (oldEndX[_myCPU] < newEndX[partner]) ?
                         oldEndX[_myCPU] : newEndX[partner];
        int minEndY   = (oldEndY[_myCPU] < newEndY[partner]) ?
                         oldEndY[_myCPU] : newEndY[partner];
        int minEndZ   = (oldEndZ[_myCPU] < newEndZ[partner]) ?
                         oldEndZ[_myCPU] : newEndZ[partner];
        int xSize = minEndX - maxStartX + 1;
        int ySize = minEndY - maxStartY + 1;
        int zSize = minEndZ - maxStartZ + 1;

        if (xSize > 0 && ySize > 0 && zSize > 0) {

            int offset = eSize *
            ((maxStartZ-oldStartZ[_myCPU]) * ArrayW(*a)*ArrayH(*a) +
             (maxStartY-oldStartY[_myCPU]) * ArrayW(*a) +
             (maxStartX-oldStartX[_myCPU]));

            MPI_Type_vector(ySize, xSize, ArrayW(*a), elem, &blk2d);
            MPI_Type_hvector(zSize, 1,
                        ArrayW(*a)*ArrayH(*a)*sSize, blk2d, &blk3d);
            MPI_Type_commit(&blk3d);
            MPI_Bsend((bdata) ?
                     ArrayPB(*a) + offset : ArrayCPB(*a) + offset,
                     1, blk3d, partner, RDIST_TAG, MPI_COMM_WORLD);
            MPI_Type_free(&blk3d);
            MPI_Type_free(&blk2d);
        }

        /*** Receive data from partner ***/

        maxStartX = (oldStartX[partner] > newStartX[_myCPU]) ?
                     oldStartX[partner] : newStartX[_myCPU];
        maxStartY = (oldStartY[partner] > newStartY[_myCPU]) ?
                     oldStartY[partner] : newStartY[_myCPU];
        maxStartZ = (oldStartZ[partner] > newStartZ[_myCPU]) ?
                     oldStartZ[partner] : newStartZ[_myCPU];
        minEndX   = (oldEndX[partner] < newEndX[_myCPU]) ?
                     oldEndX[partner] : newEndX[_myCPU];
        minEndY   = (oldEndY[partner] < newEndY[_myCPU]) ?
                     oldEndY[partner] : newEndY[_myCPU];
        minEndZ   = (oldEndZ[partner] < newEndZ[_myCPU]) ?
                     oldEndZ[partner] : newEndZ[_myCPU];
        xSize = minEndX - maxStartX + 1;
        ySize = minEndY - maxStartY + 1;
        zSize = minEndZ - maxStartZ + 1;

        if (xSize > 0 && ySize > 0 && zSize > 0) {

            int offset = eSize *
            ((maxStartZ-newStartZ[_myCPU]) * ArrayW(nA)*ArrayH(nA) +
             (maxStartY-newStartY[_myCPU]) * ArrayW(nA) +
             (maxStartX-newStartX[_myCPU]));

            MPI_Type_vector(ySize, xSize, ArrayW(nA), elem, &blk2d);
            MPI_Type_hvector(zSize, 1,
                        ArrayW(nA)*ArrayH(nA)*sSize, blk2d, &blk3d);
            MPI_Type_commit(&blk3d);
            MPI_Recv((bdata) ?
                    ArrayPB(nA) + offset : ArrayCPB(nA) + offset, 1,
                    blk3d, partner, RDIST_TAG, MPI_COMM_WORLD, &stat);
            MPI_Type_free(&blk3d);
            MPI_Type_free(&blk2d);
        }
    }

    /*** Clean up and make new local array permanent ***/

    delete oldStartX;
    delete oldStartY;
    delete oldStartZ;
    delete oldEndX;
    delete oldEndY;
    delete oldEndZ;
    delete newStartX;
    delete newStartY;
    delete newStartZ;
    delete newEndX;
    delete newEndY;
    delete newEndZ;
    MPI_Type_free(&elem);
    MPI_Buffer_detach(&mpibuff, &mpibuffsize);
    delete mpibuff;
    delete *a;
    *a = nA;
}


/**********************************************************************
 *** Array Border Data Exchange                                     ***
 **********************************************************************/

template<class ArrayT>
inline void
PxBorderExchange(ArrayT* a, int divide, int extent)
{
    int tw = ArrayW(a),  th = ArrayH(a),  td = ArrayD(a);
    int bw = ArrayBW(a), bh = ArrayBH(a), bd = ArrayBD(a);
    int cw = ArrayCW(a), ch = ArrayCH(a), cd = ArrayCD(a);

    typedef typename ArrayT::StorType storT;
    storT *ptr1, *ptr2, *ptr3, *ptr4;
    int part1, part2, xSize, ySize, zSize;
    int eSize = ArrayT::ElemSize();

    /*** Determine send/recv partners and buffer pointers & sizes ***/

    switch (divide) {
        case YZ_PLANE:                              // left <--> right
            part1 = PxMyLeftCPU();
            part2 = PxMyRightCPU();
            xSize = extent;
            ySize = ch;
            zSize = cd;
            ptr1 = ArrayPB(a) + (bd*tw*th+bh*tw+bw) * eSize;
            ptr2 = ArrayPB(a) + (bd*tw*th+bh*tw+bw+cw) * eSize;
            ptr3 = ArrayPB(a) + (bd*tw*th+bh*tw+bw+cw-xSize) * eSize;
            ptr4 = ArrayPB(a) + (bd*tw*th+bh*tw+bw-xSize) * eSize;
            break;

        case XZ_PLANE:                              // top <--> bottom
            part1 = PxMyUpCPU();
            part2 = PxMyDownCPU();
            xSize = tw;
            ySize = extent;
            zSize = cd;
            ptr1 = ArrayPB(a) + (bd*tw*th + bh*tw) * eSize;
            ptr2 = ArrayPB(a) + (bd*tw*th + (bh+ch)*tw) * eSize;
            ptr3 = ArrayPB(a) + (bd*tw*th + (bh+ch-ySize)*tw) * eSize;
            ptr4 = ArrayPB(a) + (bd*tw*th + (bh-ySize)*tw) * eSize;
            break;

        case XY_PLANE:                              // front <--> back
            part1 = PxMyFrontCPU();
            part2 = PxMyBackCPU();
            xSize = tw;
            ySize = th;
            zSize = extent;
            ptr1 = ArrayPB(a) + (bd*tw*th) * eSize;
            ptr2 = ArrayPB(a) + ((bd+cd)*tw*th) * eSize;
            ptr3 = ArrayPB(a) + ((bd+cd-zSize)*tw*th) * eSize;
            ptr4 = ArrayPB(a) + ((bd-zSize)*tw*th) * eSize;
            break;

        default:
            return;
    }

    if ((xSize == 0) || (ySize == 0) || (zSize == 0)) {
        return;
    }

    MPI_Datatype elem, blk2d, blk3d;
    MPI_Status stat;
    int sSize = eSize*sizeof(typename ArrayT::StorType);
    MPI_Type_contiguous(sSize, MPI_BYTE, &elem);
    MPI_Type_vector(ySize, xSize, tw, elem, &blk2d);
    MPI_Type_hvector(zSize, 1, tw*th*sSize, blk2d, &blk3d);
    MPI_Type_commit(&blk3d);

    /*** Send to first partner and receive from second partner ***/

    if (part1 != -1) {
        MPI_Send(ptr1, 1, blk3d, part1, BOR_TAG, MPI_COMM_WORLD);
    }
    if (part2 != -1) {
        MPI_Recv(ptr2, 1, blk3d, part2, BOR_TAG, MPI_COMM_WORLD, &stat);

    /*** Send to second partner and receive from first partner ***/

        MPI_Send(ptr3, 1, blk3d, part2, BOR_TAG, MPI_COMM_WORLD);
    }
    if (part1 != -1) {
        MPI_Recv(ptr4, 1, blk3d, part1, BOR_TAG, MPI_COMM_WORLD, &stat);
    }

    MPI_Type_free(&blk3d);
    MPI_Type_free(&blk2d);
    MPI_Type_free(&elem);
}


} // namespace Pattern
} // namespace Array
} // namespace Core
} // namespace Impala

#endif /* __PxDistribution_h_ */

---