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

[Ben Abbott]

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