From octave-maintainers-request at bevo dot che dot wisc dot edu  Tue Jun 22 12:15:36 1999
Subject: the latest LAPACK
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: "Ross A. Lippert" <ripper at cs dot sandia dot gov>
Cc: octave-maintainers at bevo dot che dot wisc dot edu
Date: Tue, 22 Jun 1999 12:15:18 -0500 (CDT)

On 22-Jun-1999, Ross A. Lippert <ripper at cs dot sandia dot gov> wrote:

| > (there is a shifting policy that you can confuse badly with the
| > right counterexample, and presently, octave's barfs on this).

Can you explain how this is `octave barfing'?  Isn't it more like
LAPACK barfing?

I believe there was also at least one bug in LAPACK's SVD calculation
that was reported here but for which I never saw a fix.  I did ask
someone to report it to the LAPACK maintainers as I did not have time
to pursue it, though I did manage to put together a Fortran-only
example that demonstrated the bug.  It is appended below.  Can you
verify that this works correctly now, and if not, can you please
report the problem?  If it is a bug in the way I am using the LAPACK
routines, can you tell me what the problem (and fix) is?

Thanks,

jwe

      program main
      integer info
      double precision a(3,6), s(3), u(3,3), vt(6,6), work(1000)
      data a /18*0.0/

      a(1,1) = 1.0
      a(2,2) = 1.0
      a(3,3) = 1.0

      a(1,4) = 1.0
      a(2,5) = 1.0
      a(3,6) = 1.0

      print *, 'a'
      print *, a

* this fails, returning sigma = [ -sqrt(2), -sqrt(2), -sqrt(2) ]
* instead of [ sqrt(2), sqrt(2), sqrt(2) ].

      call dgesvd ('N', 'N', 3, 6, a, 3, s, u, 3, vt, 6,
     $     work, 1000, info) 

      print *, 's'
      print *, s

      do j = 1, 6
        do i = 1, 3
          a(i,j) = 0.0
        enddo
      enddo

      a(1,1) = 1.0
      a(2,2) = 1.0
      a(3,3) = 1.0

      a(1,4) = 1.0
      a(2,5) = 1.0
      a(3,6) = 1.0

      print *, 'a'
      print *, a

* this works correctly, but shouldn't we be able to avoid computing
* *both* the right and left singular matrices if all we need is the
* vector of singular values?

      call dgesvd ('O', 'N', 3, 6, a, 3, s, u, 3, vt, 6,
     $     work, 1000, info) 

      print *, 's'
      print *, s

      end