utriang.f90.txt

来自「LU decomposition routines in matlab」· 文本 代码 · 共 109 行

program utriang
!
! compare the execution times for the solution of
! upper triangular systems using rthe column oriented
! version and the row oriented version of the backsolve.
!
! This program uses "second" to obtain the execution times
! can complie on SGI machines using
!  f90 utriang.f90

implicit none
integer :: i,ii,itry,j,n

real*4  :: ctime, rtime, start_time, end_time
real*4, external  :: second
real*8  :: A(1024,1024), b(1024), x(1024), b_save(1024)
real*8  :: error

n = 4

do ii = 3, 11

   n = 2*n
   
   ! create an upper triangular matrix
   A = 0.0
   x = 0.0 
   do j = 1, n
      b(j)=1.0
      do i = 1, j
         A(i,j)= sqrt((i+j)*1.0)
      end do
   end do

   ! save the right hand side
   
   b_save = b

   ! solve the upper triangular system using the COLUMN oriented version.
   ! we solve the system 10 times and average over the solution times 
   ! to 

   ctime = 0.
   
   do itry = 1, 10
   
      x = b_save
      start_time=second()

      do j = n, 1, -1
	   x(j) = x(j)/A(j,j)
           
           x(1:j-1) = x(1:j-1) - A(1:j-1,j)*x(j)
      end do
      end_time=second()
      ctime = ctime + (end_time - start_time)
   end do
  
   ctime = ctime / 10.
   
   ! check the solution
   b = b_save
   do j=1, n
      b(1:n) = b(1:n) - A(1:n,j)*x(j)
   end do

   if( sqrt(dot_product( b(1:n), b(1:n))) .gt. 1.e-10 ) then 
        print *, 'error=', sqrt(dot_product( b(1:n), b(1:n)))
   endif


   ! solve the upper triangular system using the ROW oriented version.
   ! we solve the system 10 times and average over the solution times 
   ! to 

   rtime = 0.
   
   do itry = 1, 10
   
      b = b_save
      
      start_time=second()

      x(n) = b(n)/ A(n,n)
      do i = n-1, 1, -1
	   x(i) = (b(i) - dot_product( A(i,i+1:n), x(i+1:n)))/ A(i,i)
      end do
      end_time=second()
      rtime = rtime + (end_time - start_time)
   end do
  
   rtime = rtime / 10.
  
   ! check the solution
   b = b_save
   do j=1, n
      b(1:n) = b(1:n) - A(1:n,j)*x(j)
   end do

   if( sqrt(dot_product( b(1:n), b(1:n))) .gt. 1.e-10 ) then 
        print *, 'error=', sqrt(dot_product( b(1:n), b(1:n)))
   endif
   
   print *, n, ctime, rtime
  
end do

end program utriang