@@ -2169,35 +2169,6 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in
21692169 end
21702170 return (S, U, Vᴴ), info[]
21712171 end
2172- function gesdvd! ( # SafeSVD implementation
2173- A:: AbstractMatrix{$elty} , Ac:: AbstractMatrix{$elty} = copy (A), # only modifies Ac!
2174- S:: AbstractVector{$relty} = similar (A, $ relty, min (size (A)... )),
2175- U:: AbstractMatrix{$elty} = similar (A, $ elty, size (A, 1 ), min (size (A)... )),
2176- Vᴴ:: AbstractMatrix{$elty} = similar (A, $ elty, min (size (A)... ), size (A, 2 ))
2177- )
2178- require_one_based_indexing (A, U, Vᴴ, S)
2179- chkstride1 (A, U, Vᴴ, S)
2180- m, n = size (A)
2181- minmn = min (m, n)
2182- work = Vector {$elty} (undef, 1 )
2183- if eltype (A) <: Complex
2184- if length (U) == 0 && length (Vᴴ) == 0
2185- lrwork = (LAPACK. version () <= v " 3.6" ) ? 7 * minmn : 5 * minmn
2186- else
2187- lrwork = minmn * max (5 * minmn + 5 , 2 * max (m, n) + 2 * minmn + 1 )
2188- end
2189- rwork = Vector {$relty} (undef, lrwork)
2190- else
2191- rwork = nothing
2192- end
2193- (S, U, Vᴴ), info = _gesdd_body! (Ac, S, U, Vᴴ, work, rwork)
2194- if info > 0
2195- copy! (Ac, A)
2196- (S, U, Vᴴ), info = _gesvd_body! (Ac, S, U, Vᴴ, work, rwork)
2197- end
2198- chklapackerror (info)
2199- return S, U, Vᴴ
2200- end
22012172 # ! format: off
22022173 function gesvdx! (
22032174 A:: AbstractMatrix{$elty} ,
@@ -2428,4 +2399,38 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in
24282399 end
24292400end
24302401
2402+ # SafeSVD implementation:
2403+ # attempts gesdd and falls back to gesvd, requiring two independent copies of A
2404+ # Here A is never modified, and Ac is the matrix that will be passed to LAPACK
2405+ function gesdvd! (
2406+ A:: AbstractMatrix , Ac:: AbstractMatrix{T} , # only modifies Ac!
2407+ S:: AbstractVector{Tr} = similar (A, real (T), min (size (A)... )),
2408+ U:: AbstractMatrix{T} = similar (A, T, size (A, 1 ), min (size (A)... )),
2409+ Vᴴ:: AbstractMatrix{T} = similar (A, T, min (size (A)... ), size (A, 2 ))
2410+ ) where {T <: BlasFloat , Tr <: BlasReal }
2411+ @assert Tr == real (Tr)
2412+ require_one_based_indexing (A, U, Vᴴ, S)
2413+ chkstride1 (A, U, Vᴴ, S)
2414+ m, n = size (A)
2415+ minmn = min (m, n)
2416+ work = Vector {T} (undef, 1 )
2417+ if eltype (A) <: Complex
2418+ if length (U) == length (Vᴴ) == 0
2419+ lrwork = (LAPACK. version () <= v " 3.6" ) ? 7 * minmn : 5 * minmn
2420+ else
2421+ lrwork = minmn * max (5 * minmn + 5 , 2 * max (m, n) + 2 * minmn + 1 )
2422+ end
2423+ rwork = Vector {T} (undef, lrwork)
2424+ else
2425+ rwork = nothing
2426+ end
2427+ (S, U, Vᴴ), info = _gesdd_body! (Ac, S, U, Vᴴ, work, rwork)
2428+ if info > 0
2429+ copy! (Ac, A)
2430+ (S, U, Vᴴ), info = _gesvd_body! (Ac, S, U, Vᴴ, work, rwork)
2431+ end
2432+ chklapackerror (info)
2433+ return S, U, Vᴴ
2434+ end
2435+
24312436end
0 commit comments