blas error and crash when doing ldivide/rdivide of a scalar by a vector

[Ben Abbott]

On Jan 13, 2008, at 3:19 PM, Ben Abbott wrote:

>
> On Jan 13, 2008, at 12:52 PM, John W. Eaton wrote:
>
>> On 13-Jan-2008, Ben Abbott wrote:
>>
>> | Looking over code at Netlib, http://www.netlib.org/lapack/double/dgelsd.f
>> | , the 10th parameter is RANK, which is an output.
>>
>> The error is in parameter 10 of DGEBRD, not DGELSD.
>>
>> | If I shrink lwork to 1000, then I get a similar (unrelated) error,
>> | referring to "lwork". This appears to confirm that lapack is
>> | complaining about "rank".
>> |
>> | Parameter 12 to routine DGELSD was incorrect
>> | Mac OS BLAS parameter error in DGELSD, parameter #0,
>> (unavailable), is 0
>>
>> It complains about that because it thinks LWORK is too small.
>>
>> | I've also tried g95, which gave the same result.
>>
>> It's a bug in Lapack, not the compiler.
>>
>> | Can someone familiar with how Octave does its calculation comment  
>> on
>> | what this may indicate?
>>
>> I took your example program and simplified it to be just:
>>
>>     program test
>>
>>     integer lwork, lda, ldb
>>     parameter (lwork = 10000, liwork = 11, m = 1, n = 750)
>>     integer m, n, nrhs, rank, iwork(liwork), info
>>
>>     double  precision s, a(m,n), b(n), rcond, work(lwork)
>>
>>     integer j
>>
>>     nrhs = 1
>>
>>     do j = 1, n
>>       a(1,j) = 1.0
>>     enddo
>>     b(1) = 1.0
>>
>>     call dgelsd (m, n, nrhs, a, m, b, n, s, rcond, rank,
>>    $     work, -1, iwork, info)
>>     print *, 'requested work size: ', work(1)
>>
>>     call dgelsd (m, n, nrhs, a, m, b, n, s, rcond, rank,
>>    $     work, lwork, iwork, info)
>>
>>      end
>>
>> In Octave, we use a call with lwork = -1 to get DGELSD to tell us  
>> what
>> work size to use, then we allocate space for WORK, then call DGELSD  
>> to
>> do the real work.  When I run the program above liked to ATLAS on my
>> system, it says
>>
>> $ g77 foo.f -llapack
>> $ ./a.out
>>   117.
>>  ** On entry to DGEBRD parameter number 10 had an illegal value
>>
>> so apparently it thinks WORK should be 117 elements.  But using that
>> value (by changing lwork to be 117 in the PARAMETER statement) causes
>> the
>>
>>  ** On entry to DGEBRD parameter number 10 had an illegal value
>>
>> error.
>>
>> Linking this program with the reference LAPACK code that is
>> distributed with Octave gives a different result for the expected  
>> work
>> size:
>>
>> $ g77 *.f -lblas
>> $ ./a.out
>>   741.
>>
>> However, using that value in the PARAMETER statement results in this
>> error:
>>
>>  ** On entry to DGELSD parameter number 12 had an illegal value
>>
>> It seems odd to me that DGELSD won't accept the value of LWORK that  
>> it
>> requests.  I have not tried debugging this any further to see
>> precisely what is going on.
>>
>> We can't just use a "big" value for LWORK in Octave, since the
>> required size is problem dependent.  We need to compute the necessary
>> size.  There is a formula in the comments in dgelsd.f, but it is
>> somewhat complex, so we use the query method.
>>
>> Since Octave seems to be calling DGELSD correctly to get the work
>> array size, but the returned value seems wrong, then I think the bug
>> is in LAPACK and should be fixed there.
>>
>> It might also be good to check for the same bug in ZGELSD.
>>
>> jwe
>
> I spent some time on Google ... the bug has been reported and
> apparently fixed. However, the post is from 2005, so I'm a bit
> skeptical.
>
> 	https://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=80&start=0&postdays=0&postorder=asc&highlight=dgelsd
>
> I modified jwe's short program to use ZGELSD, and it failed with the
> same error.
>
> 	Parameter 10 to routine ZGEBRD was incorrect
> 	Mac OS BLAS parameter error in ZGEBRD, parameter #0, (unavailable),
> is 0
>
> Ben

I ran several examples and found that lwork must be greater than LDB+1.

	lwork = max(nint(work(1)), ldb+2)

Where work(1) is obtained by interrogating DGELSD.

Would it be acceptable to add a check for that in Octave?

Ben