Quote:

> OK:

> linear : 0.5407 +/- 3.8e-4

> linear_lapack : 0.5364 +/- 3.8e-4

> Note that the standard deviations are rather larger than in my previous

> posts -- I was making a mistake calculating them. But the difference

> between the two is still large at 10 sigma.

> How about you write a test code where the timing relates to the matrix

> inversion only?

Here goes. I split up Davids benchmark, and added another one, as well

as converting it to f95 so that I could compile it with Lahey. Sorry

for the large number of modules, I basically cut-and-glued together

stuff that was lying around on my disk without regard for making a

small self-contained module.

Well, since what separates the wheat from the chaff in linear solvers

really is the accuracy (as someone said, what good is computing the

wrong answer really fast), I added a test for a somewhat more

ill-conditioned matrix than Davids matrix, namely the well-known

Hilbert matrix. For a matrix sufficiently big so that the timings can

actually be measured:

% ./testlinsolve

Hilbert matrix of size 500 is ill-conditioned.

rcond: 4.740663050630338E-21

Condition number: 2.109409568492821E+20

la_gesv solved system in 0.1499999999999999 seconds.

Average cumulative error is: 1577.575979219314

David Frank solved system in 2.570000000000000 seconds.

Average cumulative error is: 4.358666707056873E+17

We see that la_gesv (which uses LU) is more than an order of magnitude

faster. Well, even LU can't get reasonable result with a Hilbert

matrix this size, but for size 11 (size 12 happens to be about the

maximum la_gesv can handle in this case) we see:

% ./testlinsolve

Hilbert matrix of size 11 is well-conditioned.

rcond: 1.913805532331239E-15

Condition number: 522519129089298.4

la_gesv solved system in 0.000000000000000E+00 seconds.

Average cumulative error is: 3.605202376467460E-03

David Frank solved system in 0.000000000000000E+00 seconds.

Average cumulative error is: 18642906443420.82

While the runtimes are neglible in both cases, Franks code has some

serious issues with accuracy.

Here is the code (in addition you need blas, lapack, lapack95):

module conf

implicit none

! At least IEEE 754 double precision.

integer, parameter :: wp = selected_real_kind(15, 307)

end module conf

! ----------------------

module dfrank_linear

use conf

implicit none

contains

! dat(0 = dependent variable

! dat(1:NC = independent variables

subroutine run_frank_test ()

integer,parameter :: NS=100000 ! #data sets

integer,parameter :: NC=50 ! #coef for independent var.

real(wp) :: dat(0:NC,NS), coef(0:NC)

integer :: i

real(wp) :: time1,time2, test(0:NC) = (/ (i, i=nc+1,1,-1) /)

integer :: n

real(wp) :: a(nc+1, nc+1), b(nc+1)

test(1:nc) = test(1:nc)/100

call random_seed()

call random_number(dat)

do n = 1,NS ! embed test coef within random data sets

dat(1:NC,n) = dat(1:NC,n)-0.5d0 + test(1:NC)

dat(0,n) = test(0) + sum(dat(1:NC,n)*test(1:NC))

end do

call make_frank_matrix (dat, a, b)

call cpu_time(time1) ! start benchmark time

coef = linear(a, b, NC)

call cpu_time(time2) ! stop benchmark time

write (*,91) 'e.g. pentium 2.8ghz, CVF, time = ',time2-time1

write (*,91) 'coef recovered = ',coef(0:4), coef(nc)

write (*,91) 'coef avg error = ',SUM(abs(coef-test))/(NC+1)

91 format (a,99(f0.4,1X))

end subroutine run_frank_test

subroutine make_frank_matrix (dat, a, b)

real(wp), intent(in) :: dat(0:,:)

real(wp), intent(out) :: a(:, :), b(:)

real(wp) :: x

integer :: i, j, k, n, nc, ns

nc = size (dat, 1) - 1

print *, nc

ns = size (dat, 2)

print *, ns

a = 0. ; b = 0.

a(1,1) = ns

do n = 1,ns

do k = 1,nc

x = dat(k,n)

j = k+1

a(j,1) = a(j,1) + x

a(1,j) = a(j,1)

do i = 2,nc+1

a(i,j) = a(i,j) + x * dat(i-1,n)

end do

end do

b(1) = b(1) + dat(0,n)

do k = 1,nc

b(k+1) = b(k+1) + dat(k,n) * dat(0,n)

end do

end do

end subroutine make_frank_matrix

! -----------------------

function linear(a, b, nc) result (coef)

implicit none

integer :: nc

real(wp) :: coef(0:nc)

real(wp) :: a(nc+1,nc+1), b(nc+1), temp(nc+1), ai(nc+1,nc+1)

real(wp) :: c, d

integer :: i, j, k, m, ipvt(nc+1)

! - - - invert 'a' matrix - - -

ipvt = (/ (i, i=1,nc+1) /)

DO k = 1,nc+1

m = k-1+MAXLOC(ABS(a(k:nc+1,k)),dim=1)

IF (m /= k) THEN

ipvt( (/m,k/) ) = ipvt( (/k,m/) )

a( (/m,k/),:) = a( (/k,m/),:)

END IF

d = 1/a(k,k)

temp = a(:,k)

DO j = 1,nc+1

c = a(k,j)*d

a(:,j) = a(:,j)-temp*c

a(k,j) = c

END DO

a(:,k) = temp*(-d)

a(k,k) = d

END DO

ai = a(:,ipvt)

! - - - use inverse to evaluate coef - - -

do k = 0,nc

coef(k) = ai(k+1,1) * b(1)

do i = 1,nc

coef(k) = coef(k) +ai(k+1,i+1) * b(i+1)

end do

end do

end function linear

end module dfrank_linear

!****h* /testlinsolve

! PURPOSE

! Test linear solvers.

!****

program testlinsolve

use conf

use f95_lapack

use dfrank_linear

implicit none

! call run_frank_test ()

call run_hilbert_test ()

contains

!****t* testlinsolve/inv_cond

! PURPOSE

! Try to calculate the inverse condition number of a matrix.

!****

subroutine inv_cond (a, rcond, illcond)

real(wp), intent(inout) :: a(:,:)

real(wp), intent(out) :: rcond

logical, optional :: illcond

real(wp) :: s(size(a, 1))

!Do SVD via LAPACK95

call la_gesdd (a, s)

! Condition number of the matrix is the ratio between

! the largest to the smallest singular value. If this

! approaches infinity, the matrix is ill-conditioned,

! i.e. the input vectors were not linearly independent.

! For numerical reasons calculate the inverse of the

! condition number. If this number approaches the machine precision,

! we're screwed.

rcond = minval (s) / maxval (s)

if (rcond < epsilon (rcond)) then

if (present (illcond)) then

illcond = .true.

end if

else

if (present (illcond)) then

illcond = .false.

end if

end if

end subroutine inv_cond

!****t* testlinsolve/hilbert

! PURPOSE

! Make a Hilbert matrix, useful for testing how well linear solvers

! deal with ill-conditioned matrices.

!****

subroutine make_hilbert (a)

real(wp), intent(out) :: a(:,:)

integer :: i, j

forall (i = 1:size (a, 1), j = 1:size (a, 2))

a(i, j) = 1.0_wp / real (i + j - 1, wp)

end forall

end subroutine make_hilbert

!****t* testlinsolve/run_hilbert_test

! PURPOSE

! Test linear solvers using the Hilber matrix.

!****

subroutine run_hilbert_test ()

integer, parameter :: n = 12

real(wp) :: a(n, n), acpy(n,n), rcond, b(n), time1, &

time2, x1(n), b2(n,1), x(n)

logical :: illcond

call make_hilbert (a)

call random_seed ()

call random_number (x) !x is the correct solution

b = matmul (a, x)

b2(:,1) = b

acpy = a

! print *, 'Correct solution is: ', x

call inv_cond (a, rcond, illcond)

if (illcond) then

print *, 'Hilbert matrix of size ', n, ' is ill-conditioned.'

else

print *, 'Hilbert matrix of size ', n, ' is well-conditioned.'

end if

print *, 'rcond: ', rcond

print *, 'Condition number: ', 1.0_wp/rcond

a = acpy

call cpu_time (time1)

call la_gesv (a, b2)

call cpu_time (time2)

print *, 'la_gesv solved system in ', time2 - time1, ' seconds.'

print *, 'Average cumulative error is: ', sum (abs (b2(:,1)-x)) / n

! print *, 'Solution is: ', b2(:,1)

a = acpy

call cpu_time (time1)

x1 = linear (a, b, n-1)

call cpu_time (time2)

print *, 'David Frank solved system in ', time2 - time1, ' seconds.'

print *, 'Average cumulative error is: ', sum (abs (x1-x)) / n

! print *, 'Solution is: ', x1

end subroutine run_hilbert_test

end program testlinsolve