in LAPACK 3.1, and I think I've discovered a bug. The bug doesn't cause
incorrect results, but it can cause the iterative refinement procedure to
fail, thus meaning that DSGESV falls back to using double precision
(and therefore spoiling the good performance of the routine).
The bug is demonstrated by the program below, which calls
DSGESV to solve a trivial 2 by 2 system of equations with
two identical right-hand sides:
- Code: Select all
A = (3 6) B = (3 3)
(6 3) (3 3)
- Code: Select all
program test
implicit none
double precision a(3,2), b(3,2), x(3,2), work(3,2)
double precision xx(2,2)
integer ipiv(2)
real swork(8)
integer iter, info
external dsgesv
a(1,1) = 3
a(1,2) = 6
a(2,1) = 6
a(2,2) = 3
b(1,1) = 3
b(1,2) = 3
b(2,1) = 3
b(2,2) = 3
write (*,*) 'Calling DSGESV with N = NRHS = 2'
write (*,*) 'Call with LDX = N+1:'
call dsgesv(2,2,a,3,ipiv,b,3,x,3,work,swork,iter,info)
C
C Expected results:
C x(1,1) = 1.0d0/3
C x(1,2) = 1.0d0/3
C x(2,1) = 1.0d0/3
C x(2,2) = 1.0d0/3
C
write (*,*) 'x(1,1) = ', x(1,1)
write (*,*) 'x(1,2) = ', x(1,2)
write (*,*) 'x(2,1) = ', x(2,1)
write (*,*) 'x(2,2) = ', x(2,2)
write (*,*) 'iter = ', iter
write (*,*) 'info = ', info
if (iter.eq.-31) then
write (*,*) 'Iterative refinement failed, so double was used.'
else
write (*,*) 'Iterative refinement took ', iter, ' iterations'
end if
a(1,1) = 3
a(1,2) = 6
a(2,1) = 6
a(2,2) = 3
b(1,1) = 3
b(1,2) = 3
b(2,1) = 3
b(2,2) = 3
write (*,*)
write (*,*) 'Call with LDX = N:'
call dsgesv(2,2,a,3,ipiv,b,3,xx,2,work,swork,iter,info)
write (*,*) 'xx(1,1) = ', xx(1,1)
write (*,*) 'xx(1,2) = ', xx(1,2)
write (*,*) 'xx(2,1) = ', xx(2,1)
write (*,*) 'xx(2,2) = ', xx(2,2)
write (*,*) 'iter = ', iter
write (*,*) 'info = ', info
if (iter.eq.-31) then
write (*,*) 'Iterative refinement failed'
else
write (*,*) 'Iterative refinement took ', iter, ' iterations'
end if
end
When I run this program after linking against a build of LAPACK 3.1.1
compiled from the regular netlib source, I get this output:
- Code: Select all
Calling DSGESV with N = NRHS = 2
Call with LDX = N+1:
x(1,1) = 0.333333333333333
x(1,2) = 0.333333333333333
x(2,1) = 0.333333333333333
x(2,2) = 0.333333333333333
iter = -31
info = 0
Iterative refinement failed, so double was used.
Call with LDX = N:
xx(1,1) = 0.333333333333333
xx(1,2) = 0.333333333333333
xx(2,1) = 0.333333333333333
xx(2,2) = 0.333333333333333
iter = 2
info = 0
Iterative refinement took 2 iterations
The bug is at this call of DAXPY in DSGESV:
- Code: Select all
CALL DAXPY(N*NRHS,ONE,WORK,1,X,1)
This call assumes that the 2D array X is stored in contiguous
memory elements, which is not the case if the leading dimension
of array X is larger than N (which it is in the first example of
my demo program).
The fix is simple; just replace the call of DAXPY by these lines:
- Code: Select all
DO I = 1, NRHS
CALL DAXPY(N,ONE,WORK(1,I),1,X(1,I),1)
END DO
so that each column of X is updated independently. With this
fix, the demo program gets the same results for both calls
of DSGESV.
This bug also applies to ZCGESV, where the call of ZAXPY can
be fixed in a similar way.
