[Fortran] 纯文本查看 复制代码
program main
use f95_precision
use lapack95
implicit none
!---------------------------------------------------------------------------
integer(4) :: info
integer(4), parameter :: m = 4
integer(4) :: i,j
real(8) :: a(m,m),aa(m,m),b(m)
integer(4) :: p(m)
real(8) :: cb1,ce1,cb2,ce2
!---------------------------------------------------------------------------
data ((a(i,j),j=1,m),i=1,m) &
/ 1.8 , 2.88, 2.05, -0.89, &
5.25, -2.95, -0.95, -3.8 , &
1.58, -2.69, -2.9 , -1.04, &
-1.11, -0.66, -0.59, 0.8 /
data b /9.52, 24.35, 0.77, -6.22/
!---------------------------------------------------------------------------
! output the A, b
write(6,*) 'A ='
do i = 1,m
write(6,'(1x,<m-1>(f5.2,2x),f5.2)') (a(i,j),j=1,m)
end do
write(6,*) 'b ='
do i = 1,m
write(6,'(1x,f5.2)') b(i)
end do
!---------------------------------------------------------------------------
! guass solver
write(6,*) '=============================='
write(6,*) 'User defined solver:'
aa = a
call cpu_time(cb1)
call solve(m,aa,b)
call cpu_time(ce1)
write(6,*) 'x ='
do i = 1,m
write(6,'(1x,f5.2)') b(i)
end do
!
! b = matmul(a,b)
! write(6,*) 'Ax ='
! do i = 1,m
! write(6,'(1x,f5.2)') b(i)
! end do
write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce1-cb1
!---------------------------------------------------------------------------
! mkl solver
write(6,*) '=============================='
write(6,*) 'MKL solver:'
call cpu_time(cb2)
! f77 solution
! call dgetrf(4, 4, a, 4, p, info )
! call dgetrs('N', m, 1, a, m, p, b, m, info)
! f95 solution
call getrf(a,p)
call getrs(a,p,b)
call cpu_time(ce2)
write(6,*) 'x ='
do i = 1,m
write(6,'(1x,f5.2)') b(i)
end do
! !
! b = matmul(a,b)
! write(6,*) 'Ax ='
! do i = 1,m
! write(6,'(1x,f5.2)') b(i)
! end do
write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce2-cb2
!---------------------------------------------------------------------------
continue
end program
[Fortran] 纯文本查看 复制代码
subroutine solve(n,a,b)
!
implicit none
!
integer(kind=4) :: i,j,k,n
real(kind=8) :: a(n,n),b(n)
real(kind=8) :: x
!
do i = 1,n
continue
do j = i+1,n
x = a(j,i)/a(i,i)
a(j,i) = x
do k = i+1,n
a(j,k) = a(j,k) - x*a(i,k)
end do
end do
end do
!
do i = 2,n
do j = 1,i-1
b(i) = b(i) - b(j)*a(i,j)
end do
end do
!
b(n) = b(n)/a(n,n)
do i = n-1,1,-1
do j = n,i+1,-1
b(i) = b(i) - b(j)*a(i,j)
end do
b(i) = b(i)/a(i,i)
end do
!
return
!
end subroutine