svd - check the size of storage space

[jvdp1]

The procedure gesdd might return values for storage space larger than huge(lwork), resulting in a negative value fro lwork. Here is a proposition to solve this issue.

[perazz]

LGTM @jvdp1. Have you experienced cases where this actually happens?

Internally LAPACK is converting the integer(ilp) size request to a real value (can't permalink due to large file):

              work( 1 ) = stdlib_droundup_lwork( maxwrk )
              if( lwork<minwrk .and. .not.lquery ) then
                 info = -12
              end if

So maybe we should give a warning or return an error message when this happens?

[jvdp1]

Here is a small program to repeat the issue:

program test_svd
  use, intrinsic:: iso_fortran_env, only: wp => real64, real64
  use stdlib_stats_distribution_normal, only: rvs_normal
  use stdlib_linalg, only: svd
  implicit none
  integer, parameter :: m = 50008, n = 25153
  integer, parameter :: neig = 10
  real(wp), allocatable :: A(:,:)
  real(wp), allocatable :: tu(:,:), ts(:), tvt(:,:)

  allocate(A(m,n), source = 10._wp)

  A = rvs_normal(A, 15._wp)

  print*,'size ',m, n

 !Compute SVD(A)
 associate(k => min(m, n))
  allocate(ts(k))
  allocate(tu(m,k))
  allocate(tvt(k,n))
 end associate

 call svd(A, ts, tu, tvt, full_matrices=.false.)

 print*,'ts ',ts(1:neig)

end program

[perazz]

OK I see, thank you @jvdp1. It seems like the test program requires ~10G elements, which is more than the max number of elements that an int32 integer can store. So I think that the current fix is definitely OK to avoid overflow, but we should think about providing the ilp64 option with support for int64 integer sizes.

That could be always provided by templating the ilp integer pointer. However, I'm worried that doubling the whole implementation would slow down build times even further. On the other hand, having to decide between ilp or ilp64 at compile time is not very user friendly IMHO.

[jalvesz]

@perazz maybe one idea to balance out these cases would be to provide a flag with_linalg_ilp64 and with_linalg_ilp32_64. The former would allow replacing the ilp by ilp64, the latter, enabling both. Have yet to think this through but it might enable the choice of a compact library or a bigger one that covers everything, depending on the needs. Don't know if it should be cpp only, fypp only or a combination of both that could enable it.

[jvdp1]

Not sure what would be the best solution. Note that I got often (what I think) similar issues with the LAPACK95 interface of Intel MKL. Currently with Intel MKL, the user can only link to the 32bit or 64bit-integer libraries (AFAIK), and this is IMO quite a limitation. Having a solution as proposed by @jalvesz could be nice for the user, but could be a nightmare for the developpers/maintainers.

[perazz]

@jvdp1 @jalvesz I suggest we limit this PR to your bugfix and merge it soon.
Concerning the larger ilp64 effort, I think it deserves a separate PR.

[jvdp1]

Thank you. I agree with you Federico. Let wait a couple of days before merging it. Le mar. 22 oct. 2024 à 08:53, Federico Perini ***@***.***> a écrit :

…

@jvdp1 <https://github.com/jvdp1> I suggest we limit this PR to your bugfix and merge it soon. Concerning the larger ilp64 effort, I think it deserves a separate PR. — Reply to this email directly, view it on GitHub <#881 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD5RO7EBYCCTUGXLPYD3A7TZ4XY6VAVCNFSM6AAAAABQFLKIECVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRYGQYDEMJYGI> . You are receiving this because you were mentioned.Message ID: ***@***.***>