Closed loiseaujc closed 2 weeks ago
@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
Oups, my bad. I'll have to double check in LightKrylov
then. I'll close this issue for the moment.
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
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
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.
I think we need to fix this @loiseaujc, svd
must work for all cases, including a 1-sized matrix
Description
Consider the following MWE
Using
gfortran 13.2.0
, I get this run time error: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