dpstf2 returns incorrect value of rank

[Victor_K]

Hi lapackers,

Pivoted Cholesky subroutines seems to have a bug. In a case when diagonal elements of the input matrix are negative the return value of rank is expected to be 0 but it returns 1. The following test proves this
program test
implicit none
integer n,i,j,lda,rank,info
parameter (n=4,lda=n)
integer piv(n)
real*8 A(n,n),B(n,n),C(n,n),work(2*n)
data A /-1d0, 0d0, 0d0, 0d0,
& 0d0, -1d0, 0d0, 0d0,
& 0d0, 0d0, -1d0, 0d0,
& 0d0, 0d0, 0d0, -1d0/
data B / 0d0, 0d0, 0d0, 0d0,
& 0d0, -1d0, 0d0, 0d0,
& 0d0, 0d0, -1d0, 0d0,
& 0d0, 0d0, 0d0, -1d0/
data C / 1d0, 0d0, 0d0, 0d0,
& 0d0, -1d0, 0d0, 0d0,
& 0d0, 0d0, -1d0, 0d0,
& 0d0, 0d0, 0d0, -1d0/

print *,"test 1"
do i=1,n
print 4, (B(i,j), j=1,n)
end do
call dpstf2('L', n, B, lda, piv, rank,-1d0, work, info)
print *,"info=",info, "rank=", rank
print 44, (piv(i), i=1,n)
print *,"factored matrix:"
do i=1,n
print 4, (B(i,j), j=1,n)
end do

print *,""
print *,"test 2"
do i=1,n
print 4, (C(i,j), j=1,n)
end do
call dpstf2('L', n, C, lda, piv, rank,-1d0, work, info)
print *,"info=",info, "rank=", rank
print 44, (piv(i), i=1,n)
print *,"factored matrix:"
do i=1,n
print 4, (C(i,j), j=1,n)
end do

print *,""
print *,"test 3"
do i=1,n
print 4, (A(i,j), j=1,n)
end do
call dpstf2('L', n, A, lda, piv, rank,-1d0, work, info)
print *,"info=",info, "rank=", rank
print 44, (piv(i), i=1,n)
print *,"factored matrix:"
do i=1,n
print 4, (A(i,j), j=1,n)
end do

stop
4 format(1x, 4(D11.4,1x))
44 format(1x, 4(I2,1x))
end

The output of the test:

test 1
0.0000D+00 0.0000D+00 0.0000D+00 0.0000D+00
0.0000D+00 -0.1000D+01 0.0000D+00 0.0000D+00
0.0000D+00 0.0000D+00 -0.1000D+01 0.0000D+00
0.0000D+00 0.0000D+00 0.0000D+00 -0.1000D+01
info= 1 rank= 0
1 2 3 4
factored matrix:
0.0000D+00 0.0000D+00 0.0000D+00 0.0000D+00
0.0000D+00 -0.1000D+01 0.0000D+00 0.0000D+00
0.0000D+00 0.0000D+00 -0.1000D+01 0.0000D+00
0.0000D+00 0.0000D+00 0.0000D+00 -0.1000D+01

test 2
0.1000D+01 0.0000D+00 0.0000D+00 0.0000D+00
0.0000D+00 -0.1000D+01 0.0000D+00 0.0000D+00
0.0000D+00 0.0000D+00 -0.1000D+01 0.0000D+00
0.0000D+00 0.0000D+00 0.0000D+00 -0.1000D+01
info= 1 rank= 1
1 2 3 4
factored matrix:
0.1000D+01 0.0000D+00 0.0000D+00 0.0000D+00
0.0000D+00 -0.1000D+01 0.0000D+00 0.0000D+00
0.0000D+00 0.0000D+00 -0.1000D+01 0.0000D+00
0.0000D+00 0.0000D+00 0.0000D+00 -0.1000D+01

test 3
-0.1000D+01 0.0000D+00 0.0000D+00 0.0000D+00
0.0000D+00 -0.1000D+01 0.0000D+00 0.0000D+00
0.0000D+00 0.0000D+00 -0.1000D+01 0.0000D+00
0.0000D+00 0.0000D+00 0.0000D+00 -0.1000D+01
info= 1 rank= 1
1 2 3 4
factored matrix:
NaN 0.0000D+00 0.0000D+00 0.0000D+00
NaN NaN 0.0000D+00 0.0000D+00
NaN 0.0000D+00 -0.1000D+01 0.0000D+00
NaN 0.0000D+00 0.0000D+00 -0.1000D+01

Please look at the test case 3. Test cases 1 and 2 work well.

The same problem relates to dpstrf.

Thanks
Victor

[Julien Langou]

Hi Victor, thanks for your email and pointing this problem to us. This is fixed in the LAPACK svn repository as revision 1536. The xPSTF2/xPSTRF subroutines now perform a quick return with

Code: Select all

info= 1 rank= 0

for all of your three matrices. No NaNs are created in case #3. Cheers, Julien.