Block Cyclic Distribution in MPI

Below is an annotated discussion of the main points in your question plus some strategies you can use to organize your data and communicate the sub-blocks more elegantly.

1. Why does my local array come out column-wise?

When you use

c
MPI_Type_create_darray(
    P*Q, rank, 2, gsizes, distribs, dargs, psizes,
    MPI_ORDER_FORTRAN, MPI_DOUBLE, &filetype
);

the data is read from file in Fortran order (column-major ordering). The DARRAY type in MPI will interpret “global” indexing in Fortran fashion, and then map the piece that belongs to each process into your local buffer.

Hence, the layout of each process’s local portion A is in column-major order within that buffer. The first elements you see in memory (e.g., 1, 13, 49, 61, …) reflect how Fortran indexing jumps down columns first.

If you want the internal storage of the local array in C row-major order, you can try changing the MPI_ORDER_FORTRAN to MPI_ORDER_C. However, be aware that DARRAY with MPI_ORDER_C can sometimes yield unexpected distributions if you are mixing it with code or data arranged in Fortran style. In practice, many people leave the DARRAY as is (Fortran ordering) and then do a local reorder later, or define a different MPI datatype for the final rearranged layout.

2. How do I get 2×2 blocks written consecutively within each process’s local array?

Strictly speaking, a single MPI_Type_create_darray call will not rearrange the sub-blocks in memory so that each of your 2×2 tiles appears contiguously in row-major order. The DARRAY logic is designed to pick out which elements each rank owns, but not necessarily to reorder them in memory the way you want.

Approach A: Two-step read

1. Step 1: Use MPI_Type_create_darray(...) to read in your chunk of the matrix. You end up with a local buffer A (size = (N*N)/(P*Q) elements) in Fortran-like ordering.
2. Step 2: Rearrange locally into your “tiling” layout.
   • You already know exactly which 2×2 blocks belong to you and in what order.
   • Create a small function that loops over your local 2D region in the correct increments and writes out the data in the order you want.
   • This is often easiest: just do a local copy from the DARRAY buffer A into a new buffer A_reordered that you control.

This two-step approach is by far the most common in practice because it is simpler to debug, even though it requires a local copy.

Approach B: A more advanced single-step read

You could define a “chain” of derived datatypes:

1. One describing the entire 12×12 in block-cyclic distribution,
2. Another describing how to reorder each rank’s sub-block into 2×2 tiles,
3. Possibly combined with a subarray or vector inside the file view.

But doing so is typically far more complicated than a local rearrangement in memory. You would need careful offset/stride calculations so that each sub-block arrives in a contiguous segment. Usually, HPC programmers accept the simpler approach A.

3. Distributing (broadcasting) sub-blocks among processes

You mentioned that, after reading data, you want to send certain block-rows or block-columns from one processor to another. For instance, in a 2×3 grid:

text
+-----+-----+-----+
|  1  |  2  |  3  |
+-----+-----+-----+
|  4  |  5  |  6  |
+-----+-----+-----+

• Processor 1 should send its block-row to processor 4 (i.e., the portion of the data that corresponds to blocks horizontally adjacent to P1).
• Processor 1 should also send its block-column to processors 2 and 3, etc.

You already discovered that creating 2D Cartesian sub-communicators via MPI_Cart_sub is the right step:

c
// 2D communicator “comm2d” of size 2 x 3
int dims[2]    = {P, Q};   
int periods[2] = {0, 0};   // usually no wrap-around
MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 0, &comm2d);

// Then create row-wise communicator (height=1)
int remain_dims[2] = {0, 1}; // keep dimension 1, drop dimension 0
MPI_Comm rowcomm;
MPI_Cart_sub(comm2d, remain_dims, &rowcomm);

// Then create column-wise communicator (width=1)
remain_dims[0] = 1; // keep dimension 0, drop dimension 1
remain_dims[1] = 0;
MPI_Comm colcomm;
MPI_Cart_sub(comm2d, remain_dims, &colcomm);

How do we actually send the subarray?

There are two main ways:

1. Local copy + MPI_Bcast (or MPI_Send/MPI_Recv)

   • Extract the 2×2 (or bigger) sub-block into a contiguous temporary buffer (e.g., block_buffer).
   • Call MPI_Bcast(block_buffer, block_size, MPI_DOUBLE, root, rowcomm) or a series of MPI_Send/MPI_Recv.
   • Everyone then has that sub-block in block_buffer.

2. Use a derived datatype describing that sub-block

   • Define a derived datatype that points into your local array A at the correct offsets so that a call like

     c
     MPI_Bcast(A, 1, block_subarray_type, root, rowcomm);

     actually sends (and receives) just those 2×2 elements.
   • The advantage is “zero-copy” at the sender. The disadvantage is that all processes (the sender and receivers) must define and commit the exact same derived type so that each ranks' MPI library understands how to pick out (and/or place) the elements.

Most people opt for the first approach (local copy into a small buffer) because it is much simpler and typically not a performance bottleneck, especially for relatively small 2×2 or 4×4 blocks. Trying to maintain a derived datatype for each possible sub-block or offset is more error-prone.

Can I Bcast a subarray directly without copying?

Yes, absolutely. The MPI standard allows you to define a subarray (or vector) datatype that describes where the elements live in your local array A. For example, you could do something like:

c
MPI_Datatype sub_block_t;
// Suppose your local array is 4 x 6 in Fortran order, and you want to broadcast
// the block-row consisting of the top 2 rows and 4 columns.

int sizes[2]    = {4, 6}; // full local "matrix" shape
int subsizes[2] = {2, 4}; // shape of the sub-block to broadcast
int starts[2]   = {0, 0}; // offset from top-left corner

MPI_Type_create_subarray(2, sizes, subsizes, starts,
                         MPI_ORDER_FORTRAN, MPI_DOUBLE, &sub_block_t);
MPI_Type_commit(&sub_block_t);

// Now you can broadcast:
MPI_Bcast(A, 1, sub_block_t, root, rowcomm);

• On the sending side, MPI will “walk” that sub_block_t in memory, sending only the 2×4 region.
• On the receiving side, the same call will deposit those elements into the corresponding “positions” in A (assuming A is allocated identically on all processes).

However, if your processes do not want to place the broadcasted data at the exact same indices in memory, you either:

1. Use the same derived type but also an “origin” offset (which can be tricky).
2. Or define an equivalent type on the receiving side that references a different region.
3. Or simply do a local copy on the receiver.

Hence the simplest path is almost always:

• Root: copy sub-block to a small contiguous buffer
• MPI_Bcast that buffer with a simple (count, MPI_DOUBLE) call
• Receiver: store that buffer wherever it likes

Given the small size of blocks (2×2, 4×4, etc.), this is usually faster to implement and less error-prone.

---

Summary and Recommendations

1. Accept the Fortran ordering in the DARRAY read (or set MPI_ORDER_C if you absolutely want row-major, but keep in mind that you can still end up with a non-intuitive layout in memory).
2. If you want each 2×2 block contiguous in your local buffer, do a post-read local rearrange. This is far simpler than trying to build a fancy DARRAY that also re-tiles memory.
3. For sending sub-blocks to other ranks (block-rows or block-columns):
   • Easiest method:
     1. Extract the sub-block into a small contiguous buffer on the sender,
     2. Use MPI_Bcast (or MPI_Send/MPI_Recv) of that contiguous buffer,
     3. Place it in the desired location on the receiver.
   • More advanced method:
     • Build a custom MPI derived datatype (e.g., MPI_Type_create_subarray or MPI_Type_vector) that describes your sub-block’s location in memory.
     • Call MPI_Bcast(A, 1, sub_block_type, root, comm) to send/receive in-place.
     • Make sure every process sees the same derived type description or an equivalent mapping for correct reception.

In practice, most codes keep the DARRAY approach for reading/writing files (since it nicely hides the complexity of global distribution) and then do local copies (or a simpler derived data type) for sub-block communication. Zero-copy solutions via complicated derived datatypes are possible but rarely worth the trouble unless you are chasing very specific performance bottlenecks.