!-- From http://www.cisl.ucar.edu/docs/blueice/examples/matmul_omp.f.html

program matmul_omp
 
   implicit none
   integer, parameter :: nx = 1024, &
                         nm = nx, &
                         ny = nx

   integer, dimension(nx, nm) :: a
   integer, dimension(nm, ny) :: b
   integer, dimension(nx, ny) :: m
   
   integer :: i,j,n
 
!--Initialize the Matrix arrays

!$omp parallel do default(shared) private(n, i)
   do n=1, nm
   do i=1, nx
      a(i,n) = i+n
   end do
   end do
 
!$omp parallel do default(shared) private(n, j)
   do j=1, ny
   do n=1, nm
      b(n,j) = j-n
   end do
   end do
 
!$omp parallel do default(shared) private(i, j)
   do j=1, ny
   do i=1, nx
      m(i,j) = 0
   end do
   end do
 
!--Matrix-Matrix Multiplication

!$omp parallel do default(shared) private(i, j, n)
   do j=1,ny
   do n=1,nm
   do i=1,nx
      m(i,j) = m(i,j) + a(i,n)*b(n,j)
   end do
   end do
   end do

end program matmul_omp