|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2011 |
|
|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2011 |
Some examples using the LAPACK routine dgesv for solving a linear system.
See the Slides from lectures from Lecture 11 for more information.
See Installing LAPACK and BLAS for hints on installing.
Here is the first example with static array allocation as in Fortran 77:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 | ! $CLASSHG/codes/lapack/random/randomsys1.f90
program randomsys1
implicit none
integer, parameter :: nmax=1000
real(kind=8), dimension(nmax) :: b, x
real(kind=8), dimension(nmax,nmax) :: a
real(kind=8) :: err
integer :: i, info, lda, ldb, nrhs, n
integer, dimension(nmax) :: ipiv
! initialize random number generator seed
! if you remove this, the same numbers will be generated each
! time you run this code.
call init_random_seed()
print *, "Input n ... "
read *, n
if (n<1 .or. n>nmax) then
print *, "n must be positive and no greater than ",nmax
stop
endif
call random_number(a(1:n,1:n))
call random_number(x(1:n))
b(1:n) = matmul(a(1:n,1:n),x(1:n)) ! compute RHS
nrhs = 1 ! number of right hand sides in b
lda = nmax ! leading dimension of a
ldb = nmax ! leading dimension of b
call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
! Note: the solution is returned in b
! and a has been changed.
! compare computed solution to original x:
print *, " x computed rel. error"
do i=1,n
err = abs(x(i)-b(i))/abs(x(i))
print '(3d16.6)', x(i),b(i),err
enddo
end program randomsys1
|
Here is the code rewritten to use dynamic memory allocation:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | ! $CLASSHG/codes/lapack/random/randomsys2.f90
program randomsys2
implicit none
real(kind=8), dimension(:),allocatable :: x,b
real(kind=8), dimension(:,:),allocatable :: a
real(kind=8) :: err
integer :: i, info, lda, ldb, nrhs, n
integer, dimension(:), allocatable :: ipiv
! initialize random number generator seed
! if you remove this, the same numbers will be generated each
! time you run this code.
call init_random_seed()
print *, "Input n ... "
read *, n
allocate(a(n,n))
allocate(b(n))
allocate(x(n))
allocate(ipiv(n))
call random_number(a)
call random_number(x)
b = matmul(a,x) ! compute RHS
nrhs = 1 ! number of right hand sides in b
lda = n ! leading dimension of a
ldb = n ! leading dimension of b
call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
! Note: the solution is returned in b
! and a has been changed.
! compare computed solution to original x:
print *, " x computed rel. error"
do i=1,n
err = abs(x(i)-b(i))/abs(x(i))
print '(3d16.6)', x(i),b(i),err
enddo
deallocate(a,b,ipiv)
end program randomsys2
|
The next version also estimates the condition number of the matrix using the LAPACK routine dgecon:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | ! $CLASSHG/codes/lapack/random/randomsys3.f90
program randomsys3
implicit none
real(kind=8), dimension(:),allocatable :: x,b,work
real(kind=8), dimension(:,:),allocatable :: a
real(kind=8) :: errnorm, xnorm, rcond, anorm, colsum
integer :: i, info, lda, ldb, nrhs, n, j
integer, dimension(:), allocatable :: ipiv
integer, allocatable, dimension(:) :: iwork
character, dimension(1) :: norm
! initialize random number generator seed
! if you remove this, the same numbers will be generated each
! time you run this code.
call init_random_seed()
print *, "Input n ... "
read *, n
allocate(a(n,n))
allocate(b(n))
allocate(x(n))
allocate(ipiv(n))
call random_number(a)
call random_number(x)
b = matmul(a,x) ! compute RHS
! compute 1-norm needed for condition number
anorm = 0.d0
do j=1,n
colsum = 0.d0
do i=1,n
colsum = colsum + abs(a(i,j))
enddo
anorm = max(anorm, colsum)
enddo
nrhs = 1 ! number of right hand sides in b
lda = n ! leading dimension of a
ldb = n ! leading dimension of b
call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
! compute 1-norm of error
errnorm = 0.d0
xnorm = 0.d0
do i=1,n
errnorm = errnorm + abs(x(i)-b(i))
xnorm = xnorm + abs(x(i))
enddo
! relative error in 1-norm:
errnorm = errnorm / xnorm
! compute condition number of matrix:
! note: uses A returned from dgesv with L,U factors:
allocate(work(4*n))
allocate(iwork(n))
norm = '1' ! use 1-norm
call dgecon(norm,n,a,lda,anorm,rcond,work,iwork,info)
if (info /= 0) then
print *, "*** Error in dgecon: info = ",info
endif
print 201, n, 1.d0/rcond, errnorm
201 format("For n = ",i4," the approx. condition number is ",e10.3,/, &
" and the relative error in 1-norm is ",e10.3)
deallocate(a,b,ipiv)
deallocate(work,iwork)
end program randomsys3
|