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

[John W. Eaton]

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