Documentation issue for *lasd0

[mtowara]

Description
The documentation for the dlasd0 routine fails to mention that both U and VT are inputs and need to be initialized to Identity. Else one gets garbled or completely meaningless data.

This also applies to MKL and potentially to other LAPACK implementations.

Minimal example:

! compile with gfortran -llapack
! SUBROUTINE dlasd0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO )

program test
  Implicit None
  ! Parameters
  Integer, Parameter :: dp = kind(0.0D0)
  Integer, Parameter :: n = 6, sqre = 0, smlsiz=3
  ! Local Scalars
  Integer :: info, ldu, ldvt, lwork, i
  ! Local Arrays
  Real (Kind=dp), Allocatable :: d(:), e(:), u(:, :), vt(:, :), work(:)
  Integer, Allocatable :: iwork(:)

  ldu = n
  ldvt = n
  Allocate(d(n), e(n-1), u(ldu,n), vt(ldu,n), iwork(8*n), work(3*n**2+2*n))

  do i = 1, n
    D(i) = i+1
  end do

  ! initialize U and VT to identity
  CALL dlaset( 'Full', n, n, 0.0D0, 1.0D0, u, ldu )
  CALL dlaset( 'Full', n, n, 0.0D0, 1.0D0, vt, ldvt )

  print *, "D input: ", d(:)
  Call dlasd0(n,sqre,d,e,u,ldu,vt,ldvt,smlsiz,iwork,work,info)
  print *, "D output: ", d(:)

  print *, "U:"
  do i = 1, n
    print *, u(i,:)
  end do

  print *, "VT:"
  do i = 1, n
    print *, vt(i,:)
  end do
end program test

Without initialization (dlaset):

 D input:      2.0000000000000000        3.0000000000000000        4.0000000000000000        5.0000000000000000        6.0000000000000000        7.0000000000000000     
 D output:    6.2172489379008766E-015   7.0000000000000000        6.0000000000000000        4.0000000000000000        3.0000000000000000        2.0000000000000000     
 U:
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
 VT:
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000 

Correct behavior with dlaset:

 D input:    2.0000000000000000        3.0000000000000000        4.0000000000000000        5.0000000000000000        6.0000000000000000        7.0000000000000000     
 D output:    5.0000000000000000        7.0000000000000000        6.0000000000000000        4.0000000000000000        3.0000000000000000        2.0000000000000000     
 U:
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000     
   1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
 VT:
   0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000     
   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
   1.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue

Best,
Markus with NAG Ltd.

[langou]

Hi Markus, (1) do you want to propose a PR? (2) do you suggest the same edit for BDSDC? Julien.

[mtowara]

Sorry, forgot about this.
I can propose a change, are we okay with documenting just the usecase where U=VT=I?
(I suspect the routine might be also used with pre-factored matrices, in which case something different than I comes in)

[langou]

Let us write something like: For the description for U, "If U=I in input, then U, in output, is the left singular vector matrix of B." And, for the description for VT, "If VT=I in input, then VT, in output, is the right singular vector matrix of B."

I suspect the routine might be also used with pre-factored matrices, in which case something different than I comes in.

Maybe. I do not know. I do not think it is used this way in LAPACK, and I do not know if this is a possible usage of the routine.

[mtowara]

Sorry for dragging my feet on this.
I could not reproduce the issue for BDSDC with a quick test.
PR: #855