|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
|
|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
Some examples using the LAPACK routine dgesv for solving a linear system are in $UWHPSC/codes/lapack.
See Linear Algebra software for more about LAPACK and other software.
Here is the first example with static array allocation as in Fortran 77: Note that the lda parameter is used to indicate the size of the array that was declared, even though n may be smaller.
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 | ! $UWHPSC/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 | ! $UWHPSC/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 | ! $UWHPSC/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
|