73 template <MemoryMatrix A, MemoryVector S, MemoryMatrix U, MemoryMatrix VT>
74 requires(have_same_value_type_v<A, U, VT>
and mem::have_compatible_addr_space<A, S, U, VT>
and is_blas_lapack_v<get_value_t<A>>)
75 int gesvd(A &&a, S &&s, U &&u, VT &&vt) {
76 static_assert(has_F_layout<A>
and has_F_layout<U>
and has_F_layout<VT>,
"Error in nda::lapack::gesvd: C order not supported");
78 auto dm = std::min(a.extent(0), a.extent(1));
79 if (s.size() < dm) s.resize(dm);
82 EXPECTS(a.indexmap().min_stride() == 1);
83 EXPECTS(s.indexmap().min_stride() == 1);
84 EXPECTS(u.indexmap().min_stride() == 1);
85 EXPECTS(vt.indexmap().min_stride() == 1);
88 auto gesvd_call = []<
typename... Ts>(Ts &&...args) {
89 if constexpr (mem::have_device_compatible_addr_space<A, S, U, VT>) {
90#if defined(NDA_HAVE_DEVICE)
91 lapack::device::gesvd(std::forward<Ts>(args)...);
96 lapack::f77::gesvd(std::forward<Ts>(args)...);
101 using value_type = get_value_t<A>;
102 value_type bufferSize_T{};
103 auto rwork = array<
double, 1,
C_layout, heap<mem::get_addr_space<A>>>(5 * dm);
105 gesvd_call(
'A',
'A', a.extent(0), a.extent(1), a.data(), get_ld(a), s.data(), u.data(), get_ld(u), vt.data(), get_ld(vt), &bufferSize_T, -1,
107 int bufferSize =
static_cast<
int>(std::ceil(std::real(bufferSize_T)));
110 nda::array<value_type, 1,
C_layout, heap<mem::get_addr_space<A>>> work(bufferSize);
111 gesvd_call(
'A',
'A', a.extent(0), a.extent(1), a.data(), get_ld(a), s.data(), u.data(), get_ld(u), vt.data(), get_ld(vt), work.data(), bufferSize,