!Robert W. Numrich, Introduction to Co-Array Fortran
! http://www2.hpcl.gwu.edu/pgas09/tutorials/caf_tut.pdf
!
program matmul
implicit none
real, allocatable,dimension(:,:), codimension[:,:] :: a,b,c
integer :: i
integer :: j
integer :: k
integer :: l
integer,parameter :: n = 10
integer :: p
integer :: q
integer :: iAm
integer :: myP
integer :: myQ
p = num_images()
q = int(sqrt(float(p)))
iAm = this_image()
if (q*q /= p) then
if(iAm == 1) write (*,"('num_images must be square: p=',i5)") p
stop
end if
allocate(a(n,n)[q,*])
allocate(b(n,n)[q,*])
allocate(c(n,n)[q,*])
myP = this_image(c,1)
myQ = this_image(c,2)
a = 1.0
b = 1.0
c = 0.0
sync all
do i=1,n
do j=1,n
do k=1,n
do l=1,q
c(i,j) = c(i,j) + a(i,k)[myP, l]*b(k,j)[l,myQ]
end do
end do
end do
end do
if (any(c /= n*q)) write(*,"('error on image: ',2i5,e20.10)") myP, myQ, c(1,1)
write(*,"('check sum[',i5',',i5,']',e20.10)") myP, myQ, sum(c) - q*n**3
deallocate(a,b,c)
end program matmul
