LU Decomposition MPI

Question

This is a MPI code for LU Decomposition.

I have used the following strategy -

There is a master(rank 0) and others are slaves. The master sends rows to each slave. Since each slave might receive more than row, I store all the received rows in a buffer and then perform LU Decomposition on it. After doing that I send back the buffer to the master. The master does not do any computation. It just sends and receives.

for(i=0; i<n; i++)
    map[i] = i%(numProcs-1) + 1;

for(i=0; i<n-1; i++)
{
    if(rank == 0)
    {
        status = pivot(LU,i,n);

        for(j=0; j<n; j++)
            row1[j] = LU[n*i+j];
    }

    MPI_Bcast(&status, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(status == -1)
        return -1;

    MPI_Bcast(row1, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    int tag1 = 1, tag2 = 2, tag3 = 3, tag4 = 4;

    if(rank == 0)
    {
        int pno, start, index, l, rowsReceived = 0;
        MPI_Request req;
        MPI_Status stat;

        for(j=i+1; j<n; j++)
            MPI_Isend(&LU[n*j], n, MPI_DOUBLE, map[j], map[j], MPI_COMM_WORLD, &req);

        if(i>=n-(numProcs-1))
            cnt++;

        for(j=0; j<numProcs-1-cnt; j++)
        {
            MPI_Recv(&pno, 1, MPI_INT, MPI_ANY_SOURCE, tag2, MPI_COMM_WORLD, &stat);
            //printf("1. Recv from %d and j : %d and i : %d\n",pno,j,i);
            MPI_Recv(&rowsReceived, 1, MPI_INT, pno, tag3, MPI_COMM_WORLD, &stat);
            MPI_Recv(rowFinal, n*rowsReceived, MPI_DOUBLE, pno, tag4, MPI_COMM_WORLD, &stat);

            /* Will not go more than numProcs anyways */
            for(k=i+1; k<n; k++)
            {
                if(map[k] == pno)
                {
                    start = k;
                    break; 
                }
            }

            for(k=0; k<rowsReceived; k++)
            {
                index = start + k*(numProcs-1);

                for(l=0; l<n; l++)
                    LU[n*index+l] = rowFinal[n*k+l];
            }
        }
    }

    else
    {
        int rowsReceived = 0;
        MPI_Status stat, stats[3];
        MPI_Request reqs[3];

        for(j=i+1; j<n; j++)
            if(map[j] == rank)
                rowsReceived += 1;


        for(j=0; j<rowsReceived; j++)
        {
            MPI_Recv(&rowFinal[n*j], n, MPI_DOUBLE, 0, rank, MPI_COMM_WORLD, &stat);
        }

        for(j=0; j<rowsReceived; j++)
        {
            double factor = rowFinal[n*j+i]/row1[i];

            for(k=i+1; k<n; k++)
                rowFinal[n*j+k] -= (row1[k]*factor);

            rowFinal[n*j+i] = factor;
        }

        if(rowsReceived != 0)
        {
            //printf("Data sent from %d iteration : %d\n",rank,i);
            MPI_Isend(&rank, 1, MPI_INT, 0, tag2, MPI_COMM_WORLD, &reqs[0]);
            MPI_Isend(&rowsReceived, 1, MPI_INT, 0, tag3, MPI_COMM_WORLD, &reqs[1]);
            MPI_Isend(rowFinal, n*rowsReceived, MPI_DOUBLE, 0, tag4, MPI_COMM_WORLD, &reqs[2]);
        }
        //MPI_Waitall(3,reqs,stats);
    }
}

The problem that I am facing is that sometimes the program hangs. My guess is that the sends and receives are not being matched but I am not being able to figure out where the problem lies.

I ran test cases on matrices of size 1000x1000, 2000x2000, 3000x3000, 5000x5000 and 7000x7000. Presently the code hangs for 7000x7000. Could someone please help me out?

Things to note :-

1. map implements the mapping scheme, which row goes to which slave.

2. rowsReceived tells each slave the no of rows it will receive. I dont need to calculate that each and every time, but I will fix it later.

3. row1 is the buffer in which the active row will be stored.

4. rowFinal is the buffer of the rows being received and being modified.

5. cnt is not important and can be ignored. For that the check for rowReceived!=0 needs to be removed.

Answer 1

It looks like you are never completing your nonblocking operations. You have a bunch of calls to MPI_Isend and MPI_Irecv throughout the code, but you're never doing an MPI_Wait or MPI_Test (or one of the similar calls). Without that completion call, those nonblocking calls will never complete.