Public Pastes
Guest

Angel de Vicente

By: a guest on Jul 20th, 2010  |  syntax: Fortran  |  size: 7.24 KB  |  hits: 84  |  expires: Never
download  |  raw  |  embed  |  report abuse
Copied
  1. PROGRAM FOX
  2.   IMPLICIT NONE
  3.   include "mpif.h"
  4.  
  5.   INTEGER :: procs, rank, error
  6.   INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
  7.  
  8.   INTEGER :: g_order, g_side, my_g_row, my_g_column, my_g_rank, comm, row_comm, col_comm, block_mpi_t
  9.   INTEGER :: rows
  10.   REAL, DIMENSION(:,:), ALLOCATABLE :: matrixA, matrixB, matrixC, localA, tempA, localB, localC
  11.  
  12.   CALL MPI_Init ( error )
  13.   CALL MPI_Comm_size ( MPI_COMM_WORLD, procs, error )
  14.   CALL MPI_Comm_rank ( MPI_COMM_WORLD, rank, error )
  15.  
  16.   CALL Read_Matrix(rank,error,rows,matrixA,matrixB,matrixC)
  17.   CALL Setup_grid(g_order,g_side,error,comm,my_g_rank,my_g_row,my_g_column,row_comm,col_comm)
  18.   CALL Distribute_matrixes(g_side,rows,block_mpi_t,error,matrixA,matrixB,localA,localB,localC,my_g_row,my_g_column)
  19.   CALL Perform_Fox_Algorithm(my_g_row,my_g_column,g_order,localA,tempA,localB,localC,row_comm,col_comm,status,error)
  20.   CALL Print_Last_Result(rank, my_g_rank, g_side, comm, procs, error, status, localC, matrixC, tempA)
  21.  
  22.   CALL MPI_Finalize (error)
  23.  
  24. CONTAINS
  25.  
  26.   SUBROUTINE Print_Last_Result(rank, my_g_rank, g_side, comm, procs, error, status, localC, matrixC, tempA)
  27.     INTEGER, INTENT(IN) :: rank, my_g_rank, g_side, comm, procs
  28.     INTEGER, INTENT(OUT) :: error
  29.     INTEGER, DIMENSION(:), INTENT(OUT) :: status
  30.     REAL, DIMENSION(:,:), INTENT(IN) :: localC
  31.     REAL, DIMENSION(:,:), INTENT(OUT) :: matrixC, tempA
  32.  
  33.     INTEGER :: grow,gcol,i
  34.     INTEGER,DIMENSION(2) :: coordinates
  35.  
  36.     IF (rank .EQ. 0 .AND. my_g_rank .NE. 0) PRINT*, "Houston, we have a problem with I/O"
  37.  
  38.     CALL MPI_SEND(localC,g_side*g_side,MPI_REAL,0,0,comm,error)
  39.  
  40.     IF (rank .EQ. 0) THEN
  41.        matrixC=0      
  42.        DO i=1,procs
  43.           CALL MPI_Recv(tempA,g_side*g_side,MPI_REAL,MPI_ANY_SOURCE,0,comm,status,error)
  44.           CALL MPI_Cart_coords(comm,status(MPI_Source),2,coordinates,error)
  45.           grow = coordinates(1)
  46.           gcol = coordinates(2)
  47.           matrixC(grow*g_side+1:(grow+1)*g_side,gcol*g_side+1:(gcol+1)*g_side) = tempA
  48.        END DO
  49.     PRINT*,"============ Matrix multiplication in parallel ==================="
  50.     CALL Print_Matrix(matrixC)
  51.     END IF
  52.   END SUBROUTINE Print_Last_Result
  53.  
  54.  
  55.   SUBROUTINE Print_Matrix(M)
  56.     REAL, DIMENSION(:,:) :: M
  57.     INTEGER :: side,row,column
  58.  
  59.     side = SIZE(M,1)
  60.  
  61.     DO row=1,side
  62.        DO column=1,side-1
  63.           WRITE(*,'(F10.3)',ADVANCE='NO'), M(row,column)
  64.        END DO
  65.        WRITE(*,'(F10.3)'), M(row,side)
  66.     END DO
  67.   END SUBROUTINE Print_Matrix
  68.  
  69.  
  70.   SUBROUTINE Perform_Fox_Algorithm(my_g_row,my_g_column,g_order,localA,tempA,localB,localC,row_comm,col_comm,status,error)
  71.     INTEGER, INTENT(IN) :: my_g_row,my_g_column,g_order,row_comm,col_comm
  72.     INTEGER, INTENT(OUT) :: error
  73.     INTEGER, DIMENSION(:), INTENT(OUT) :: status
  74.     REAL, DIMENSION(:,:), INTENT(IN) :: localA
  75.     REAL, DIMENSION(:,:), INTENT(INOUT) :: tempA,localB,localC
  76.  
  77.     INTEGER :: source,destination,stage,bcast_root
  78.  
  79.     localC = 0
  80.     source = MOD(my_g_row + 1,g_order)
  81.     destination = MOD(my_g_row + g_order - 1,g_order)
  82.  
  83.     DO stage=0,g_order-1
  84.        bcast_root = MOD(my_g_row + stage,g_order)
  85.  
  86.        IF (my_g_column .EQ. bcast_root) THEN
  87.           CALL MPI_BCAST(localA,g_side*g_side,MPI_REAL,bcast_root,row_comm,error)
  88.           CALL Matrix_Multiply(localA,localB,localC)
  89.        ELSE
  90.           CALL MPI_BCAST(tempA,g_side*g_side,MPI_REAL,bcast_root,row_comm,error)
  91.           CALL Matrix_Multiply(tempA,localB,localC)
  92.        END IF
  93.        CALL MPI_Sendrecv_replace(localB,g_side*g_side,MPI_REAL,destination,0,source,0,col_comm,status,error)
  94.     END DO
  95.   END SUBROUTINE Perform_Fox_Algorithm
  96.  
  97.  
  98.   SUBROUTINE Matrix_Multiply(A,B,C)
  99.     REAL, DIMENSION(:,:), INTENT(IN) :: A,B
  100.     REAL, DIMENSION(:,:), INTENT(OUT) :: C
  101.     INTEGER :: side,row,column
  102.  
  103.     side = SIZE(A,1)
  104.  
  105.     DO row=1,side
  106.        DO column=1,side
  107.           C(row,column) = C(row,column) + SUM(A(row,:)*B(:,column))
  108.        END DO
  109.     END DO
  110.   END SUBROUTINE Matrix_Multiply
  111.  
  112.  
  113.   SUBROUTINE Distribute_matrixes(g_side,rows,block_mpi_t,error,matrixA,matrixB,localA,localB,localC,my_g_row,my_g_column)
  114.     INTEGER, INTENT(IN) :: g_side,rows,my_g_row,my_g_column
  115.     INTEGER, INTENT(OUT) :: block_mpi_t,error
  116.     REAL, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: localA,localB,localC
  117.     REAL, DIMENSION(:,:), INTENT(INOUT) :: matrixA,matrixB
  118.  
  119.     CALL MPI_TYPE_VECTOR(g_side,g_side,rows,MPI_REAL,block_mpi_t,error)
  120.     CALL MPI_TYPE_COMMIT(block_mpi_t, error)
  121.  
  122.     !! We send the whole matrixA and matrixB. Should find an efficient way
  123.     !! of sending the matrixes directly to localA and localB
  124.  
  125.     CALL MPI_BCAST(matrixA,rows*rows,MPI_REAL,0,MPI_COMM_WORLD,error)
  126.     CALL MPI_BCAST(matrixB,rows*rows,MPI_REAL,0,MPI_COMM_WORLD,error)
  127.  
  128.     ALLOCATE(localA(g_side,g_side),tempA(g_side,g_side),localB(g_side,g_side),localC(g_side,g_side))
  129.     localA = matrixA(my_g_row*g_side+1:(my_g_row+1)*g_side,my_g_column*g_side+1:(my_g_column+1)*g_side)
  130.     localB = matrixB(my_g_row*g_side+1:(my_g_row+1)*g_side,my_g_column*g_side+1:(my_g_column+1)*g_side)
  131.   END SUBROUTINE Distribute_matrixes
  132.  
  133.  
  134.   SUBROUTINE Setup_grid(g_order,g_side,error,comm,my_g_rank,my_g_row,my_g_column,row_comm,col_comm)
  135.     INTEGER, INTENT(OUT) :: g_order, g_side,error,comm,my_g_rank,my_g_row,my_g_column,row_comm,col_comm
  136.  
  137.     INTEGER,DIMENSION(2) :: dimensions,coordinates
  138.     LOGICAL,DIMENSION(2) :: wrap_around,free_coords
  139.  
  140.     g_order = SQRT(REAL(procs))
  141.     g_side = rows / g_order
  142.  
  143.     dimensions = g_order
  144.     wrap_around = .TRUE.
  145.     CALL MPI_Cart_create(MPI_COMM_WORLD,2,dimensions,wrap_around,.TRUE.,comm,error)
  146.     CALL MPI_Comm_rank ( comm, my_g_rank, error )
  147.     CALL MPI_Cart_coords(comm,my_g_rank,2,coordinates,error)
  148.     my_g_row = coordinates(1)
  149.     my_g_column = coordinates(2)
  150.  
  151.     free_coords(1) = .FALSE.
  152.     free_coords(2) = .TRUE.
  153.     CALL MPI_Cart_sub(comm,free_coords,row_comm,error)
  154.  
  155.     free_coords(1) = .TRUE.
  156.     free_coords(2) = .FALSE.
  157.     CALL MPI_Cart_sub(comm,free_coords,col_comm,error)
  158.   END SUBROUTINE Setup_grid
  159.  
  160.  
  161.   SUBROUTINE Read_Matrix(rank,error,rows,matrixA,matrixB,matrixC)
  162. !!!
  163. !!!  Read the matrix in process 0, and generate the multiplication in Process 0 (for testing)
  164. !!!    
  165.     INTEGER, INTENT(IN) :: rank
  166.     INTEGER, INTENT(OUT) :: rows,error
  167.     INTEGER :: i
  168.     REAL, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: matrixA,matrixB,matrixC
  169.  
  170.     IF (rank .EQ. 0) THEN
  171.        READ*, rows
  172.        CALL MPI_BCAST(rows,1,MPI_REAL,0,MPI_COMM_WORLD,error)
  173.        ALLOCATE(matrixA(rows,rows),matrixB(rows,rows),matrixC(rows,rows))
  174.        DO i=1,rows
  175.           READ*, matrixA(i,:)
  176.        END DO
  177.        DO i=1,rows
  178.           READ*, matrixB(i,:)
  179.        END DO
  180.        
  181.        !! Calculate the matrix multiplication for comparison
  182.        matrixC = 0
  183.        CALL Matrix_Multiply(matrixA,matrixB,matrixC)
  184.  
  185.        PRINT*, "CALCULATED IN PROCESS 0"
  186.        PRINT*, "=======MATRIX A=========="
  187.        CALL Print_Matrix(matrixA)
  188.        PRINT*, "=======MATRIX B=========="
  189.        CALL Print_Matrix(matrixB)
  190.        PRINT*, "=======MATRIX C=========="
  191.        CALL Print_Matrix(matrixC)
  192.     ELSE
  193.        CALL MPI_BCAST(rows,1,MPI_REAL,0,MPI_COMM_WORLD,error)
  194.        ALLOCATE(matrixA(rows,rows),matrixB(rows,rows))
  195.     END IF
  196.   END SUBROUTINE Read_Matrix
  197.  
  198. END PROGRAM FOX
create a new version of this paste RAW Paste Data