!----------------------------------------------------------------!
! The current code has three purpose:                            !
! (1) how to parallalize nested loops.                           !
! (2) matrix-matrix product.                                     !
! (3) algorithm's importance.                                    !
!----------------------------------------------------------------!
! Xiaoxu Guan at LSU HPC.                                        !
! March 17, 2016                                                 !
!----------------------------------------------------------------!

      program nested_loops_v2
      implicit none
      integer            :: i,j,k,nosize
      integer, parameter :: idp = 8, nsize = 6000
      real(kind=idp)     :: temp,time_start,time_end,omp_get_wtime
      real(kind=idp), dimension(1:nsize,1:nsize) :: a,b,c


      nosize = nsize
! initialize arrays: a, b, and c:
      do i = 1, nosize
      do j = 1, nosize
         a(j,i) = SIN( dble(i-j) )
         b(j,i) = COS( dble(i+j) )
      end do
      end do

      c = 0.0_idp


! c = a x b 

      time_start = omp_get_wtime()

! only parallelize the outer loop in jki:
!$omp parallel do default(none) private(i,j,k) shared(nosize,a,b,c)
      do j = 1, nosize
      do k = 1, nosize
         do i = 1, nosize
            c(i,j) = c(i,j) + a(i,k) * b(k,j)
         end do
      end do
      end do
!$omp end parallel do

      time_end = omp_get_wtime()

      write(*,*) 'elapsed time (seconds) = ', time_end - time_start
      write(*,*)

!----------------------------------------------------------------------------!




      end program nested_loops_v2

