High Performance Scientific Computing   AMath 483/583 Class Notes   Spring Quarter, 2011

#### Previous topic

Installing LAPACK and BLAS

#### Next topic

Parallel Computing

# LAPACK examplesΒΆ

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 ```