!
! *****************************************************************
! *                                                               *
! *    Copyright (c) Digital Equipment Corporation, 1995          *
! *                                                               *
! *   All Rights Reserved.  Unpublished rights  reserved  under   *
! *   the copyright laws of the United States.                    *
! *                                                               *
! *   The software contained on this media  is  proprietary  to   *
! *   and  embodies  the  confidential  technology  of  Digital   *
! *   Equipment Corporation.  Possession, use,  duplication  or   *
! *   dissemination of the software and media is authorized only  *
! *   pursuant to a valid written license from Digital Equipment  *
! *   Corporation.                                                *
! *                                                               *
! *   RESTRICTED RIGHTS LEGEND   Use, duplication, or disclosure  *
! *   by the U.S. Government is subject to restrictions  as  set  *
! *   forth in Subparagraph (c)(1)(ii)  of  DFARS  252.227-7013,  *
! *   or  in  FAR 52.227-19, as applicable.                       *
! *                                                               *
! *****************************************************************
!
!
! HISTORY
!
! File:  lu.f90
!
! This program times the standard LU-decomposition kernel on a 1000x1000
! array distributed (*,CYCLIC).

program lu

parameter(n = 1000)

real, dimension(n, n) :: a
integer i, j, k
double precision start_time, finish_time

!HPF$ distribute(*, cyclic) :: a

! Arbitrary initialization of the array a:
forall (i = 1:n, j = 1:n)  a(i, j) = (i+j)*(i+j)*(i+j)

print 100, n, n
100   format ('Computing ', I4, 'x', I4, ' LU decomposition...')

call timer(start_time)

do k = 1, n-1
  a(k+1:n, k) = a(k+1:n, k)/a(k, k)
  forall (i = k+1:n, j = k+1:n) a(i, j) = a(i, j) - a(i, k)*a(k, j)
end do
      
call timer(finish_time)

print 101, finish_time - start_time
101 format('Elapsed time in seconds:  ',f10.3)

end


subroutine timer(t)

integer val(8)
double precision t

call date_and_time(values=val)
t = dble(val(8))*1d-3 + dble(val(7)) +          &
    dble(val(6))*60d0 + dble(val(5))*3600d0

end subroutine
