Subdividing Scalapack grid

Question

I am trying to calculate the eigenspectrum for a number of large matrices using Scalapack, but rather than distribute each matrix across all 32 processes. I'd rather distribute each matrix across 4 processes and calculate 8 matrices in parallel. I know how to subdivide an MPI Grid using MPI_Comm_split, but I it seems Scalapack doesn't take custom communicators. Instead it seems to use a BLACS grid rooted in PVM.

How can I implement this subdivision in Scalapack?

Answer 1

This is done by BLACS and grid setup.

The reference functions are

• BLACS_GRIDINIT( ICONTXT, ORDER, NPROW, NPCOL )
• BLACS_GRIDMAP( ICONTXT, USERMAP, LDUMAP, NPROW, NPCOL )

The documentation for these routines states:

These routines take the available processes, and assign, or map, them into a BLACS process grid.

Each BLACS grid is contained in a context (its own message passing universe), so that it does not interfere with distributed operations which occur within other grids/contexts.

These grid creation routines may be called repeatedly in order to define additional contexts/grids.

Which means that you can create 8 different grids and pass each ICONTXT to scalapack routines for each matrix.

Both of them get an IN/OUT argument

ICONTXT

(input/output) INTEGER

On input, an integer handle indicating the system context to be used in creating the BLACS context. The user may obtain a default system context via a call to BLACS_GET. On output, the integer handle to the created BLACS context.

You can use these contexts recursively on the same manner.

Answer 2

I implemented @ztik suggestion and this is the result I came up with. It seems to work:

program main
    use mpi
    implicit none
    integer              :: ierr, me, nProcs, color, i,j,k, my_comm, dims(2), global_contxt
    integer              :: cnt, n_colors, map(2,2)
    integer, allocatable :: contxts(:)
    integer, parameter   :: group_size = 4

    call MPI_Init(ierr)

    call MPI_Comm_size(MPI_COMM_WORLD, nProcs, ierr)
    call MPI_Comm_rank(MPI_COMM_WORLD, me, ierr)

    color = get_color(group_size)
    n_colors =  nProcs / group_size
    allocate(contxts(n_colors))

    dims = calc_2d_dim(group_size)

    call BLACS_GET(0, 0, global_contxt)
    if(me == 0) write (*,*) global_contxt
    contxts = global_contxt

    do k = 1,n_colors 
        ! shift for each context
        cnt = group_size * (k-1)
        if(me==0) write (*,*) "##############", cnt

        ! create map
        do i=1,2
            do j=1,2
                map(i,j) = cnt
                cnt = cnt + 1
            enddo
        enddo

        call BLACS_GRIDMAP(contxts(k), map, 2, 2, 2)
        do i = 0,nProcs
            if(i == me) then
                write (*,*) me, contxts(k)
            endif
            call MPI_Barrier(MPI_COMM_WORLD, ierr)
        enddo
    enddo

    call MPI_Finalize(ierr)
contains
    function get_color(group_size) result(color)
        implicit none
        integer, intent(in)     :: group_size
        integer                 :: me, nProcs, color, ierr, i, cnt

        call MPI_Comm_size(MPI_COMM_WORLD, nProcs, ierr)
        call MPI_Comm_rank(MPI_COMM_WORLD, me, ierr)

        if(mod(nProcs, group_size) /= 0) then
            write (*,*) "Nprocs not divisable by group_size", mod(nProcs, group_size)
            call MPI_Abort(MPI_COMM_WORLD, 0, ierr)
        endif

        color = 0
        do i = 1,me
            if(mod(i, group_size) == 0) then
                color = color + 1
            endif
        enddo
    end function get_color

    function calc_2d_dim(sz) result(dim)
        implicit none
        integer, intent(in)    :: sz
        integer                :: dim(2), cand

        cand = nint(sqrt(real(sz)))

        do while(mod(sz, cand) /= 0)
            cand = cand - 1
        enddo
        dim(1) = sz/cand
        dim(2) = cand
    end function calc_2d_dim

end program main