## dpstf2 returns incorrect value of rank

Post here if you want to report a bug to the LAPACK team

### dpstf2 returns incorrect value of rank

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
Victor_K

Posts: 5
Joined: Sat Nov 27, 2010 9:48 am

### Re: dpstf2 returns incorrect value of rank

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.
Julien Langou

Posts: 828
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA