/home/marwa/GE/ReportsN3000/NoOptNoVec/luOmp.f90

Line% of fetchesSource
1  
!-- OpenMp Do version
2  
program lu
3  
  !uncomment when using ifort for rand()
4  
  !use IFPORT
5  
  implicit none
6  
  integer, parameter :: DP=kind(0.0D0)
7  
8  
!--  Variables
9  
  integer :: i,j,k,nthr,id,myseed=7654321,n
10  
  real(kind=DP), dimension(:,:), ALLOCATABLE :: A, U, B
11  
  real(kind=DP), dimension(:), ALLOCATABLE ::LHS, X
12  
  integer, dimension(:), ALLOCATABLE :: pvt
13  
  real(kind=DP) :: timer,error,walltime,ipvt, resd,swap, mx, mn
14  
  integer :: OMP_GET_MAX_THREADS
15  
  
16  
  nthr=OMP_GET_MAX_THREADS()
17  
  timer=walltime()
18  
19  
!-- Read n and Allocate A
20  
   !switch to variable n when do sp
21  
  write (*,*) 'Enter n: '
22  
  read *,n
23  
  ALLOCATE(A(n,n+1), LHS(n), U(n,n), B(n,n), pvt(n),X(n))
24  
  i=rand(myseed) 
25  
!-- Initiate matrix
26  
 do j=1,n
27  
     do i=1,n
28  
        A(i,j)= rand(0)
29  
        B(i,j)=A(i,j)
30  
     end do
31  
     pvt(j)=j
32  
end do
33  
 LHS(1:n)=0.0
34  
 A(1:n,n+1)=1
35  
36  
 timer=walltime()
37  
!-- LU-factorize
38  
  do k=1,n-1
39  
   !----Find Pivot k
40  
    do i=k+1,n
41  
      if ( abs(A(i,k))> abs(A(pvt(k),k)) ) then
42  
         pvt(k)=i;
43  
      end if
44  
     end do
45  
46  
  !---------------end Find Pivot
47  
  if (pvt(k)>k) then
48  
   swap=A(k,k); A(k,k)=A(pvt(k),k); A(pvt(k),k)=swap
49  
  end if
50  
  do i=k+1,n
51  
    A(i,k)=A(i,k)/A(k,k)
52  
  end do
53  
54  
!$OMP PARALLEL DO private(i,j,swap)
55  
  do j=k+1,n+1
56  
   swap=A(pvt(k),j)
57  
   if (pvt(k)>k) then
58  
     A(pvt(k),j)=A(k,j); A(k,j)=swap
59  
   end if
60  
   do i=k+1,n
61 98.3%
     A(i,j)=A(i,j)-A(i,k)*swap
62  
   end do
63  
  end do
64  
   
65  
end do
66  
67  
  timer=walltime()-timer
68  
69  
  write(*,*) 'n = ',n,'   time = ',timer, 'nthreads= ', nthr 
70  
71  
  ! CHECK CORRECTNESS
72  
  
73  
  do j=1,n
74  
     U(j,j)=A(j,j)
75  
     do i=j+1,n
76  
        U(i,j)=0
77  
     end do
78  
     do i=1,j-1
79  
        U(i,j)=A(i,j)
80  
     end do
81  
  end do
82  
83  
  X(n)=A(n,n+1)/A(n,n)
84  
  do j=n,2, -1
85  
     do i=1,j-1
86  
        LHS(i)=LHS(i)+U(i,j)*X(j)
87  
     end do
88  
     i=j-1
89  
     X(i)=(A(i,n+1)-LHS(i))/A(i,i)
90  
  end do
91  
  
92  
  error=0.0
93  
  do i=1,n
94  
     resd=0.0
95  
     do j=1,n
96  
        resd=resd+B(i,j)*X(j) 
97  
     end do
98  
     error=error+abs(1-resd)   !err=abs(b-AX), b=1
99  
  end do
100  
101  
  write(*,*) 'ERROR: ',error
102  
  DEALLOCATE(A,B,LHS,X,U,pvt)
103  
104  
CONTAINS
105  
 subroutine PrintMat(A,n)
106  
    implicit none
107  
    REAL(kind=DP), dimension(:,:), INTENT (IN)   :: A
108  
    integer, INTENT (IN)   :: n
109  
    integer :: i
110  
    do i=1,n
111  
        print *, A(i,:)
112  
    end do
113  
    return
114  
  end subroutine PrintMat
115  
116  
end program lu
117  
118  

Copyright (c) 2006-2011 Rogue Wave Software, Inc. All Rights Reserved.
Patents pending.