# Example : Matrix Transpose

We will demonstrate the use of virtual topologies by way of a matrix transposition.
The matrix algebra for matrix transpose is:

The parallel algorithm is:

1. Select p and q such that the total number of processes, nprocs = p x q.
2. Partition the n x m matrix into a (blocked) p x q matrix whose elements are themselves matrices of size (n/p) x (m/q).
3. Perform a transpose on each of these sub-matrices. These are performed serially as the entire sub-matrix resides locally on a process. No inter-process communication is required.
4. Formally, the p x q matrix need be transposed to obtain the final result. However, in reality this step is often not necessary. If you need to access the element (or sub-matrix) “p,q” of the transposed matrix, all you need do is to access the element “q,p” which has already been transposed locally. Depending on what comes next in the calculation, unnecessary message passing may be avoided.
i,j (p)
aij

As an example, we define a 9 x 4 matrix and will use 6 processes. Next, map it into a 3 x 2 virtual cartesian grid, i.e., p=3, q=2. Coincidentally, each element of this cartesian grid is in turn a matrix of size 3 x 2, as illustrated in the figure below. For the physical grid, each square box represents one entry of the matrix. The pair of indices, “i,j”, on the first row gives the global cartesian coordinates while “(p)” is the process associated with the virtual grid allocated by calling `MPI_Cart_create` or `MPI_Comm_split`. aij on the second row is the value of the matrix element. On the right, the 3 x 2 virtual grid is depicted. Each box in this grid represents one process and contains one 3×2 submatrix. Finally, another communicator is created for the transposed virtual grid which has the dimensions of 2 x 3. The element at “1,0” of the transposed virtual grid, for instance, stores the value sent by the element at “0,1” of the virtual grid.

0,0 (0)
100
0,1 (0)
101
1,0 (0)
110
1,1 (0)
111
2,0 (0)
120
2,1 (0)
121
0,2 (1)
102
0,3 (1)
103
1,2 (1)
112
1,3 (1)
113
2,2 (1)
122
2,3 (1)
123

0,0 (0) 0,1 (1)
1,0 (2) 1,1 (3)
2,0 (4) 2,1 (5)
3,0 (2)
130
3,1 (2)
131
4,0 (2)
140
4,1 (2)
141
5,0 (2)
150
5,1 (2)
151
3,2 (3)
132
3,3 (3)
133
4,2 (3)
142
4,3 (3)
143
5,2 (3)
152
5,3 (3)
153
6,0 (4)
160
6,1 (4)
161
7,0 (4)
170
7,1 (4)
171
8,0 (4)
180
8,1 (4)
181
6,2 (5)
162
6,3 (5)
163
7,2 (5)
172
7,3 (5)
173
8,2 (5)
182
8,3 (5)
183

### Virtual Grid

0,0 (0) 0,1 (1)
1,0 (2) 1,1 (3)
2,0 (4) 2,1 (5)
``````      program matrix_transpose
implicit none
integer n, m, nv, nl, mv, ml, i, il, iv, j, jl, jv
integer p, ndim, reorder, ierr, grid_comm
integer master, me, Iam, source, dest, tag
parameter (n=9, m=8, nv=3, mv=2, nl=n/nv, ml=m/mv)
parameter (ndim=2, reorder=1)
integer a(nl,ml), at(ml, nl), b(m,n)
include "mpif.h"  !! This brings in pre-defined MPI constants, ...
integer dims(ndim), coord(ndim), req
logical period(ndim)
integer status(MPI_STATUS_SIZE)

data master/0/    !! 0 is defined as the master processor
data period/.false.,.false./  !! no cyclic boundary in either index
data tag/0/       !! a tag is not required in this case, set it to zero
data dest/0/      !! results are sent back to master

c**Starts MPI processes ...
call MPI_Init(ierr)                            !! starts MPI
call MPI_Comm_rank(MPI_COMM_WORLD, Iam, ierr)  !! get current process id
call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)    !! get number of processes

c**create cartesian topology for matrix
dims(1) = nv
dims(2) = mv
call MPI_Cart_create(MPI_COMM_WORLD, ndim, dims,
&       period, reorder, grid_comm, ierr)
call MPI_Comm_rank(grid_comm, me, ierr)
call MPI_Cart_coords(grid_comm, me, ndim, coord, ierr)
iv = coord(1)
jv = coord(2)

c**define local matrix according to virtual grid coordinates, (iv,jv)
do jl=1,ml
do il=1,nl
i = il + iv*nl
j = jl + jv*ml
a(il,jl) = i*10 + j
enddo
enddo

c**perform transpose on local matrix
do jl=1,ml
do il=1,nl
at(jl,il) = a(il,jl)
enddo
enddo

c**send "at" to Master for asembly and printing
call MPI_Isend(at, ml*nl, MPI_INTEGER, Master, tag,
&       grid_comm, req, ierr)

c**Master asembles all local transposes into final matrix and print
if(Iam .eq. Master)call asemble(at,ml,nl,grid_comm,b,m,n,p)

call MPI_Wait(req, status, ierr)       !! make sure all sends done

call MPI_Finalize(ierr)                !! let MPI finish up ...

end
subroutine asemble(at,ml,nl,comm,b,m,n,p)
implicit none
include "mpif.h"
integer ml, nl, comm, tag, source, m, n, ierr, p, ndim
integer iv, jv, i, j, il, jl, coord(2)
integer b(m,n), at(ml,nl), status(MPI_STATUS_SIZE)
data tag, ndim/0,2/

c**The Master asembles the final (transposed) matrix from local copies and print
do source=0,p-1
call MPI_Cart_coords(comm, source, ndim, coord, ierr)
call MPI_Recv(at,ml*nl, MPI_INTEGER, source, tag, comm,
&       status, ierr)
iv = coord(1)
jv = coord(2)
do jl=1,nl
j = jl + iv*nl       ! swap iv and jv for transpose
do il=1,ml
i = il + jv*ml
b(i,j) = at(il,jl)
enddo
enddo
enddo
write(*,'(9i5)')((b(i,j),j=1,n),i=1,m)

return
end``````