Program gauss_elimination

parameter(n=400)
double precision, Allocatable,Dimension(:,:)  :: a
double precision, Allocatable,Dimension(:)      :: x
double precision time_begin,time_end

!HPF$ PROCESSORS PROC(NUMBER_OF_PROCESSORS())
!HPF$ ALIGN X(i) with A(i,*)
!HPF$ DISTRIBUTE (CYCLIC,*) ONTO PROC :: A 

interface
  Subroutine gauss_part_piv(a,x,n)
    double precision, Dimension(n,n+1)  :: a 
    double precision, Dimension(n)      :: x
    !HPF$ ALIGN X(i) with *A(i,*)
    !HPF$ DISTRIBUTE *(CYCLIC,*)  :: A 
  end subroutine
end interface

Allocate(A(n,n+1))
Allocate(X(n))
call  Generate_A_x(A,X,n)
call timer(time_begin)
call gauss_part_piv(a,x,n)
call timer(time_end)
print *,'Elapsed time for solving linear equation with dimension of'
print *,'coeefient matrix',n,'by ',n,'is '
print *,time_end-time_begin,' seconds'

end


!Subroutine for Solving linear equation A.X = B
! by Gauss elimination with partial pivoting

! Parameters
! A : Input,Coefficient matrix with dimension (n,n+1). Right hand side vector
!     of linear equation is stored into last column,(n+1)'th, of matrix A
! X : Unknown vector . This subroutine finds the unknown X values and
!     return them in X parameter
! n : Number of linear equation

Subroutine gauss_part_piv(a,x,n)
double precision, Dimension(n,n+1)  :: a
double precision, Dimension(n)      :: x
Integer :: i,j,k,temp
Real :: p(1),m,piv

!HPF$ ALIGN X(i) with *A(i,*)
!HPF$ DISTRIBUTE *(CYCLIC,*)  :: A 

!HPF INDEPENDENT,new(i,j,m)
Do i=1,n
   Do j=i+1,n
     m = a(j,i)/a(i,i)
     a(j,i:n+1) = a(j,i:n+1) - m*a(i,i:n+1)
   enddo
!   a(i,:) = a(i,:)/a(i,i)
enddo

!Set x(n)
! First calculate x(n) solving A(nrow(n),n+1) * X(n) = A(nrow(n),n)

x(n) = a(n,n+1)/a(n,n)

! Solve X(n-1),X(n-2),..X(1) by backward substitution

!HPF INDEPENDENT,new(i)
Do i=n-1,1,-1
  x(i) = (a(i,n+1)-sum(a(i,i+1:n)*x(i+1:n)))/a(i,i)
enddo
return

end

Subroutine Generate_A_x(A,X,n)
double precision A(n,n+1),X(n)
integer n,i


Call Random_number(A)
Call Random_number(X)
A = Aint(30*A)+1
X = Aint(30*X)+1
Do i=1,n
  A(i,n+1) = sum(A(i,1:n) * X)
enddo
return
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 timer


Subroutine print_mat(A,n)
Real A(n,n+1)
integer n,i
Do i=1,n
  print *,A(i,1:n+1)
enddo
end