svd not working for matrices with a single column

[loiseaujc]

Description

Consider the following MWE

program main
  use stdlib_kinds, only: dp
  use stdlib_linalg, only: svd
  implicit none

  integer, parameter :: m = 2, n = 1
  ! Sigular value decomposition of A.
  real(kind=dp) :: A(n, n)
  real(kind=dp) :: U(n, n), S(n), Vt(n, n)

  ! Random matrix.
  call random_number(A) 

  ! Call to stdlib.
  call svd(A, S, U, Vt)

end program main

Using gfortran 13.2.0, I get this run time error:

At line 75528 of file build/dependencies/stdlib/src/stdlib_linalg_lapack_d.F90
Fortran runtime error: Index '2' of dimension 1 of array 'a' above upper bound of 1

Error termination. Backtrace:
#0  0x5621220d2c4a in __stdlib_linalg_lapack_d_MOD_stdlib_dgesdd
	at build/dependencies/stdlib/src/stdlib_linalg_lapack_d.F90:75528
#1  0x56212190442b in __stdlib_linalg_MOD_stdlib_linalg_svd_d
	at build/dependencies/stdlib/src/stdlib_linalg_svd.f90:483
#2  0x5621218f77c1 in MAIN__
	at app/main.f90:15
#3  0x5621218f780c in main
	at app/main.f90:2
<ERROR> Execution for object " MWE_stdlib_lstsq " returned exit code  2
<ERROR> *cmd_run*:stopping due to failed executions

Expected Behaviour

Computing the SVD of a column vector is admittedly a contrived example but this piece of code is part of larger subroutine in LightKrylov iteratively computing the SVD of large-scale matrices using Lanczos Bidiagonalization. Note that if a row vector is considered instead of a column vector, the MWE runs perfectly.

Pinging @perazz, @jalvesz, @jvdp1.

Version of stdlib

0.6.1

Platform and Architecture

Linux with Ubuntu 22.04

Additional Information

No response

[jalvesz]

@loiseaujc I noticed that in your mwe you use only n to define the matrices. shouldn't m be the first dimension of every matrix? tried that here https://godbolt.org/z/vf6bxqezd and it works

[loiseaujc]

Oups, my bad. I'll have to double check in LightKrylov then. I'll close this issue for the moment.

[loiseaujc]

Actually, if m = 1, and n = 1, I still get an error.

program main
  use stdlib_kinds, only: dp
  use stdlib_linalg, only: svd
  implicit none

  integer, parameter :: m = 1, n = 1
  ! Sigular value decomposition of A.
  real(kind=dp) :: A(m, n)
  real(kind=dp) :: U(m, m), S(n), Vt(n, n)

  ! Random matrix.
  call random_number(A) 

  ! Call to stdlib.
  call svd(A, S, U, Vt)

end program main

[perazz]

Sorry all for the late reply, I'm abroad on business travel the whole week. The error apparently comes from a LAPACK limitation, as a(2,1) is always accessed: in DGESDD,

                    ! produce r in a, zeroing out other entries
                    call stdlib_dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )

but it should not take place as the requirement is LDA >= max(1,M). So, I will investigate what matrix option is being used for this edge case.

EDIT: Yes: that code segment is called part of this branch:

           if( m>=n ) then
              ! a has at least as many rows as columns. if a has sufficiently
              ! more rows than columns, first reduce using the qr
              ! decomposition (if sufficient workspace available)

so all provided matrix sizes are OK according to the inputs. stdlib_dlaset should return without doing anything because the input n-1 == 0, but of course the compiler's bound checking complains because we are referencing address a(2,1), which does not exist. So it indeed looks like you've found an issue with LAPACK: not technically a bug, but rather a F77-style edge case that modern compilers complain about. We must check if any of the other options (I.e. the ‘reduced matrices’ one) allows to continue bypassing this error

[loiseaujc]

Oh I see. I was using gesvd in LightKrylov which doesn't seem to suffer from this issue. I'll switch to stdlib_svd nonetheless and simply add a if (k > 1) call svd(A, S, U, Vt) since in practice it is very unlikely that a Krylov method converges in a single itertion. No big deal. Closing this issue then.

[perazz]

I think we need to fix this @loiseaujc, svd must work for all cases, including a 1-sized matrix