c c c********************************************************************* c c Example 21 (mm_mult_cartesian.f) c c Objective : Matrix Matrix multiplication(Using Cartesian Topology) c c Input : Read files (mdata1.inp) for first input matrix c and (mdata2.inp) for second input matrix c c Output : Result of matrix matrix multiplication on Processor 0. c c Necessary Condition : Number of Processes should be less than c or equal to 8. Vector size for Vectors A and c B should be properly striped. that is Vector c size should be properly divisible by c Number of processes used. c c********************************************************************* program main include 'mpif.h' include 'cart.h' integer A_Bloc_MatrixSize, B_Bloc_MatrixSize integer NoofRows_A, NoofCols_A, NoofRows_B, NoofCols_B integer NoofRows_BlocA, NoofCols_BlocA integer NoofRows_BlocB, NoofCols_BlocB integer Matrix_Size(4) integer Local_Index, Global_Row_Index, Global_Col_Index integer Local_varb, irow, icol, jrow integer Local_iproc, Local_irow integer iproc, jproc, index, Proc_Id, Root integer Numprocs, MyRank, MyGlobalRank, root_p real Matrix_A(ROW_SIZE_A,COL_SIZE_A), $ Matrix_B(ROW_SIZE_B,COL_SIZE_B), $ Matrix_C(ROW_SIZE_A,COL_SIZE_B) real MatA_array(TOTAL_SIZE_A), MatB_array(TOTAL_SIZE_B), $ MatC_array(TOTAL_SIZE_C) real A_Bloc_Matrix(TOTAL_SIZE_A), $ B_Bloc_Matrix(TOTAL_SIZE_B), $ C_Bloc_Matrix(TOTAL_SIZE_C) real LocalMatrix_A(TOTAL_SIZE_A) real LocalMatrix_B(TOTAL_SIZE_B) real RowBlock_A(TOTAL_SIZE_A), ColBlock_B(TOTAL_SIZE_B) c For Setting Cartesian Topology c Periods - For Wraparound in each dimension. c Dimensions - Number of processors in each dimension. c Coordinates - processor Row and Column identification c Remain_dims - For row and column communicators. integer Periods(2) integer Dimensions(2) integer Coordinates(2) integer Remain_dims(2) integer Comm2D, Row_comm, Col_comm c .....Initialising call MPI_INIT (ierror) call MPI_COMM_SIZE(MPI_COMM_WORLD, Numprocs, ierror) call MPI_COMM_RANK(MPI_COMM_WORLD, MyGlobalRank, ierror) c ........Set up the MPI_COMM_WORLD and CARTESIAN TOPOLOGY Dimensions(1) = sqrt(dble(Numprocs)) Dimensions(2) = sqrt(dble(Numprocs)) root_p = sqrt(dble(Numprocs)) c ........Wraparound mesh in both dimensions. Periods(1) = 1 Periods(2) = 1 c ........Create Cartesian topology in two dimnesions and c Cartesian decomposition of the processes call MPI_DIMS_CREATE( Numprocs, 2, Dimensions, ierr ) call MPI_CART_CREATE(MPI_COMM_WORLD, 2, Dimensions, $ Periods, .TRUE., Comm2D, ierror) call MPI_COMM_RANK(Comm2D,MyRank,ierror) call MPI_CART_COORDS(Comm2D,MyRank,2,Coordinates,ierror) Row = Coordinates(1) Col = Coordinates(2) c Construction of row communicator and column communicators c (use cartesian row and columne machanism to get c Row and Col Communicators) Remain_dims(1) = 0 Remain_dims(2) = 1 c The output communicator represents the column c containing the process call MPI_CART_SUB(Comm2D, Remain_dims, Row_comm, ierror) Remain_dims(1) = 1 Remain_dims(2) = 0 c The output communicator represents the row containing c the process call MPI_CART_SUB(Comm2D, Remain_dims, Col_comm, ierror) c Reading Input Root = 0 if ( MyRank .eq. Root) then c ......Read the Matrix From Input file ..... open(unit=12, file = './data/mdata1.inp') read(12,*) NoofRows_A, NoofCols_A write(6,*) 'Input matrix - A ' do i = 1,NoofRows_A read(12,*) (Matrix_A(i,j),j=1,NoofCols_A) write(6,75) (Matrix_A(i,j),j=1,NoofCols_A) enddo open(unit=13, file = './data/mdata2.inp') read(13,*) NoofRows_B, NoofCols_B write(6,*) 'Input matrix - B' do i = 1,NoofRows_B read(13,*) (Matrix_B(i,j),j=1,NoofCols_B) write(6,75) (Matrix_B(i,j),j=1,NoofCols_B) enddo close(12) close(13) Matrix_Size(0) = NoofRows_A Matrix_Size(1) = NoofCols_A Matrix_Size(2) = NoofRows_B Matrix_Size(3) = NoofCols_B endif c ..........Send Matrix Size to all processors call MPI_BCAST (NoofRows_A, 1, MPI_INTEGER, Root, $ MPI_COMM_WORLD, ierror) call MPI_BCAST (NoofCols_A, 1, MPI_INTEGER, Root, $ MPI_COMM_WORLD, ierror) call MPI_BCAST (NoofRows_B, 1, MPI_INTEGER, Root, $ MPI_COMM_WORLD, ierror) call MPI_BCAST (NoofCols_B, 1, MPI_INTEGER, Root, $ MPI_COMM_WORLD, ierror) c ........ print*,'MyRank',MyRank,'A= ',NoofRows_A,NoofCols_A if(NoofCols_A .ne. NoofRows_B) then if(MyRank .eq. Root) $ print*," Matrix Dimension is in correct for Multiplication" goto 100 endif c rem_Rows_A = mod(NoofRows_A, Dimensions(1)) c rem_Cols_A = mod(NoofCols_A, Dimensions(1)) c c rem_Rows_B = mod(NoofRows_B, Dimensions(1)) c rem_Cols_B = mod(NoofCols_B, Dimensions(1)) c if((mod(NoofRows_A, Dimensions(1)) .ne. 0) .or. $ (mod(NoofCols_A, Dimensions(1)) .ne. 0) .or. $ (mod(NoofRows_B, Dimensions(1)) .ne. 0) .or. $ (mod(NoofCols_B, Dimensions(1)) .ne. 0)) then if(MyRank .eq. Root) $ print*,"Matrices can't be divided among processors equally" goto 100 endif NoofRows_BlocA = NoofRows_A / Dimensions(1) NoofCols_BlocA = NoofCols_A / Dimensions(1) NoofRows_BlocB = NoofRows_B / Dimensions(1) NoofCols_BlocB = NoofCols_B / Dimensions(1) A_Bloc_MatrixSize = NoofRows_BlocA * NoofCols_BlocA B_Bloc_MatrixSize = NoofRows_BlocB * NoofCols_BlocB c print*,'MyRank',MyRank,NoofRows_BlocA, NoofCols_BlocA c ....Memory allocating for Bloc Matrices c ....Rearrange the input matrices in one dimensional c arrays by approriate order if(MyRank .eq. Root) then c ...... Convert Matrix into one-dimensional array c according to Checkerboard Local_Index = 1 do iproc = 0, root_p-1 do jproc = 0, root_p-1 do irow = 1, NoofRows_BlocA Global_Row_Index = iproc * NoofRows_BlocA + irow do icol = 1, NoofCols_BlocA Global_Col_Index = jproc * NoofCols_BlocA + icol MatA_array(Local_Index) = $ Matrix_A(Global_Row_Index, Global_Col_Index) Local_Index = Local_Index + 1 enddo enddo enddo enddo Local_Index = 1 do iproc = 0, root_p-1 do jproc = 0, root_p-1 do irow = 1, NoofRows_BlocB Global_Row_Index = iproc * NoofRows_BlocB + irow do icol = 1, NoofCols_BlocB Global_Col_Index = jproc * NoofCols_BlocB + icol MatB_array(Local_Index) = $ Matrix_B(Global_Row_Index, Global_Col_Index) Local_Index = Local_Index + 1 enddo enddo enddo enddo endif c Scatter the Data to all processes by MPI_SCATTER call MPI_SCATTER(MatA_array, A_Bloc_MatrixSize, MPI_REAL, $ A_Bloc_Matrix, A_Bloc_MatrixSize, MPI_REAL, $ Root, Comm2D, ierror) call MPI_SCATTER(MatB_array, B_Bloc_MatrixSize, MPI_REAL, $ B_Bloc_Matrix, B_Bloc_MatrixSize, MPI_REAL, $ Root, Comm2D, ierror) c .........Broadcasting Data amongst rows and columns c call MPI_ALLGATHER(A_Bloc_Matrix, A_Bloc_MatrixSize, $ MPI_REAL, RowBlock_A, A_Bloc_MatrixSize, MPI_REAL, $ Row_comm, ierror) call MPI_ALLGATHER (B_Bloc_Matrix, B_Bloc_MatrixSize, $ MPI_REAL, ColBlock_B, B_Bloc_MatrixSize, MPI_REAL, $ Col_comm, ierror) call MPI_BARRIER(Comm2D, ierror) c ..Rearrange Data to do matrix multiplication for RowBlock_A index = 1 do irow = 0, NoofRows_BlocA-1 Local_irow = irow*NoofCols_BlocA do iproc = 0, root_p-1 Local_iproc = iproc * A_Bloc_MatrixSize do icol = 1, NoofCols_BlocA Local_varb = icol + Local_irow + Local_iproc LocalMatrix_A(index) = RowBlock_A(Local_varb) index = index + 1 enddo enddo enddo c ... Rearrange Data to do matrix multiplication for ColBlock_B index = 1 do icol = 1, NoofCols_BlocB do iproc = 0, root_p-1 Local_iproc = iproc*B_Bloc_MatrixSize do irow = 0, NoofRows_BlocB - 1 Local_varb = icol + Local_iproc + irow*NoofCols_BlocB LocalMatrix_B(index) = ColBlock_B(Local_varb) index = index + 1 enddo enddo enddo C ........Multiply LocalBlock matrices to get C_Bloc_matrix index = 1 do irow = 0, NoofRows_BlocA-1 do icol = 0, NoofCols_BlocB-1 C_Bloc_Matrix(index) = 0.0 do jrow = 1, NoofRows_B C_Bloc_Matrix(index) = C_Bloc_Matrix(index) + $ LocalMatrix_A(irow*NoofCols_A + jrow) * $ LocalMatrix_B(icol*NoofRows_B+jrow) enddo index = index + 1 enddo enddo c ........Gather output block matrices at processor 0 call MPI_GATHER(C_Bloc_Matrix, NoofRows_BlocA*NoofCols_BlocB, $ MPI_REAL, MatC_array, NoofRows_BlocA * NoofCols_BlocB, $ MPI_REAL, Root, Comm2D, ierror) c ...Rearranging the output matrix in a array by approriate order if (MyRank .eq. Root) then do iproc = 0, root_p-1 do jproc = 0, root_p-1 Proc_Id = iproc * root_p + jproc do irow = 1, NoofRows_BlocA Global_Row_Index = iproc * NoofRows_BlocA + irow do icol = 1, NoofCols_BlocB Global_Col_Index = jproc * NoofCols_BlocB + icol Local_Index = (Proc_Id*NoofRows_BlocA*NoofCols_BlocB) $ + ((irow-1) * NoofCols_BlocB) + icol Matrix_C(Global_Row_Index,Global_Col_Index) = $ MatC_array(Local_Index) enddo enddo enddo enddo write(6,*) 'Output matrix ' do i = 1,NoofRows_A write(6,75) (Matrix_C(i,j),j=1,NoofCols_B) enddo 75 format(8(2x,f8.3)) endif 100 call MPI_FINALIZE(ierror) stop end