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

[Ben Abbott]

On Jan 13, 2008, at 6:05 AM, Carlo de Falco wrote:

>
> On 13/gen/08, at 00:18, John W. Eaton wrote:
>
>> | doesn't matlab rely on LAPACK too?
>>
>> I don't think it is doing precisely the same calculation for right  
>> and
>> left divide when the matrix is not square.
>
> the result of the fortran version coincides with that given by  
> matlab, though
>
>> jwe
>
> c.
>

Carlo, it looks to me like your fortran isn't doing the same  
calculation.

I'm running Octave 3.0.0+ on Mac OS X (Leopard), compiled with gfortran.

The short Octave script produced the output below

for n=730:750
   A = ones (1, n);
   B = 1;
   X = A\B;
   disp ( sprintf ('n = %d: X = %10.8f', n, X(1) ) )
endfor

octave:1> test
Parameter 10 to routine DGEBRD was incorrect
Mac OS BLAS parameter error in DGEBRD, parameter #0, (unavailable), is 0
octave(58928) malloc: *** error for object 0x66d6b0: double free
*** set a breakpoint in malloc_error_break to debug
n = 730: X = 0.00136986
n = 731: X = 0.00136799
n = 732: X = 0.00136612
n = 733: X = 0.00136426
n = 734: X = 0.00136240
n = 735: X = 0.00136054
n = 736: X = 0.00135870
n = 737: X = 0.00135685
n = 738: X = 0.00135501
n = 739: X = 0.00135318

Octave exits for n=740 (with the "Parameter 10 to routine DGEBRD was  
incorrect" message). I am surprised that the problem is with parameter  
10 of DGEBRD.

Looking over code at Netlib, http://www.netlib.org/lapack/double/dgelsd.f 
, the 10th parameter is RANK, which is an output.

In any event, I've modified Carlo's short program to do the equivalent  
(at least I think it is) in Fortran. It appears to work correctly.

$ cat test.f
       program test

       implicit none

       integer lwork, lda, ldb
       parameter (lwork = 10000000)
       parameter (lda = 1)
       parameter (ldb = 1000)
       integer m, n, nrhs, rank, iwork(lwork), info

       double  precision s, a(lda,ldb), b(ldb), rcond, work(lwork)

       integer ii, jj

       m = 1
       nrhs = 1

       rcond = 1.0e-7

       do n = 730,750
         do ii=1,m
            do jj=1,n
               a(ii,jj) = 1.0
            enddo
            b(ii) = 1.0
         enddo

         call dgelsd( m, n, nrhs, a, lda, b, ldb, s, rcond, rank,
      $     work, lwork, iwork, info )

         write(*,*) 'n = ', n, ' b=', b(1)
       enddo

       end

$ gfortran -c -g -o test.o test.f
$ gfortran -framework Accelerate -o test test.o
$ ./test
  n =          730  b=  1.369863013698631E-003
  n =          731  b=  1.367989056087544E-003
  n =          732  b=  1.366120218579236E-003
  n =          733  b=  1.364256480218276E-003
  n =          734  b=  1.362397820163487E-003
  n =          735  b=  1.360544217687082E-003
  n =          736  b=  1.358695652173912E-003
  n =          737  b=  1.356852103120761E-003
  n =          738  b=  1.355013550135502E-003
  n =          739  b=  1.353179972936401E-003
  n =          740  b=  1.351351351351351E-003
  n =          741  b=  1.349527665317138E-003
  n =          742  b=  1.347708894878702E-003
  n =          743  b=  1.345895020188427E-003
  n =          744  b=  1.344086021505375E-003
  n =          745  b=  1.342281879194634E-003
  n =          746  b=  1.340482573726542E-003
  n =          747  b=  1.338688085676040E-003
  n =          748  b=  1.336898395721924E-003
  n =          749  b=  1.335113484646190E-003
  n =          750  b=  1.333333333333332E-003

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

I've also tried g95, which gave the same result.

Can someone familiar with how Octave does its calculation comment on  
what this may indicate?

TiA
Ben
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20080113/1eca5de9/attachment.html