program cannon_on_square_grid
　　include 'mpif.h'
　　integer myid, nproc, comm, ierr
　　integer dest, sour, left, right, upper, lower
　　integer maxlen, maxn, nps, n, status(mpi_status_size)
　　parameter (n=400,maxlen=n*n+1)
　　real*8 a(maxlen), b(maxlen), c(maxlen)
　　real*8 t1, t2, flops, perfs

  Initialize communicator.
　　call mpi_init( ierr )
　　call mpi_comm_dup( mpi_comm_world, comm, ierr )
　　call mpi_comm_rank( comm, rnyid, ierr )
　　call mpi_comm_size( comm, nproc, ierr )
　　nps = sqrt(nproc*1.0)
　　if(nps*nps. ne. nproc. or. mod(n, nps). ne. 0)goto 999
　　m=n/nps
 
 Initialize matrices.
　　myrow = mod(myid,nps)
　　mycol = myid/nps
　　do 10 j = 1, m
　　do 10 i = 1, m
　　　a(i+(j-1)*m) = (myrow*m+i)*n+mycol*m+j
　　　b(i+(j-1)*m) = (myrow*m+i)*n*n+mycol*m+j
　10 continue

　　t1=mpi_wtime()

 A1ignment of matrix A.
　　dest = myrow + mod(mycol-myrow+nps,nps)*nps
　　sour = myrow + mod(mycol+myrow,nps)*nps
　　call mpi_sendrecv_replace(a, m*m, mpi_double_precision, dest, 10, & sour, 10, comm, status, ierr)

 A1ignment of matrix B.
　　dest = mod(myrow-mycol + nps, nps) + mycol*nps
　　sour = mod(mycol + myrow, nps) + mycol*nps
　　call mpi_sendrecv_replace(b, m*m, mpi_double_precision, dest, 20, & sour, 20, comm, status, ierr)

 Initialize of matrix C.
　　do 20 j = 1, m
　　do 20 i = 1, m
　　c(i + (j-1)*m)=0.0
　20 continue

 Main procedure.
　　left = myrow + mod(mycol-1+nps,nps)*nps
　　upper = mod(myrow-1+nps,nps) + mycol*nps
　　right = myrow + mod(mycol+1,nps)*nps
　　lower = mod(myrow+1,nps) + mycol*nps
　　do 30 k=0, nps-1
　　　call dgemm('n', 'n', m, m, m, 1.0D0, a, m, b, m, 1.0D0, c, m)
　　　call mpi_sendrecv_replace(a, m*m, mpi_double_precision, left, & k+1, right, k+1, comm, status, ierr)
　　　call mpi_sendrecv_replace(b, m*m, mpi_double_precision, upper, & k+1+nps, lower, k+1+nps comm, status, ierr)
　30 continue

 Reverse alignment of matrix A.
　　sour = myrow + mod(mycol-myrow+nps, nps)*nps
　　dest = myrow + mod(mycol+myrow, nps)*nps
　　call mpi_sendrecv_replace(a, m*m, mpi_double_precision, dest, 100, ＆ sour, 100, comm, status, ierr)

 Reverse alignment of matrix B.
　　sour = mod(myrow-mycol+nps, nps) + mycol*nps
　　dest = mod(mycol+myrow, nps) + mycol*nps
　　call mpi_sendrecv_replace(b, m*m, mpi_double_precision, dest, 200, & sour, 200, comm, status, ierr)

　　t2 = mpi_wtime()
　　flops = 2.0d-6*n*n*n
　　perfs = flops/(t2-t1)
　　write(*,*)n, nprow, npcol, myid, perfs
　999 call mpi_comm_free(comm, ierr)
　　call mpi_finalize(ierr)
　　stop
　　end