Function svd

src/backend/opencl/svd.cpp:65–195

template<typename T, typename Tr>
void svd(Array<T> &arrU, Array<Tr> &arrS, Array<T> &arrVT, Array<T> &arrA,
         bool want_vectors = true) {
    dim4 idims    = arrA.dims();
    dim4 istrides = arrA.strides();

    const int m      = static_cast<int>(idims[0]);
    const int n      = static_cast<int>(idims[1]);
    const int ldda   = static_cast<int>(istrides[1]);
    const int lda    = m;
    const int min_mn = std::min(m, n);
    const int ldu    = m;
    const int ldvt   = n;

    const int nb    = magma_get_gebrd_nb<T>(n);
    const int lwork = (m + n) * nb;

    cpu_lapack_lacpy_func<T> cpu_lapack_lacpy;
    cpu_lapack_bdsqr_work_func<T> cpu_lapack_bdsqr_work;
    cpu_lapack_ungbr_work_func<T> cpu_lapack_ungbr_work;
    cpu_lapack_lamch_func<Tr> cpu_lapack_lamch;

    // Get machine constants
    static const double eps    = cpu_lapack_lamch('P');
    static const double smlnum = std::sqrt(cpu_lapack_lamch('S')) / eps;
    static const double bignum = 1. / smlnum;

    Tr anrm = abs(getScalar<T>(reduce_all<af_max_t, T, T>(arrA)));

    T scale                = scalar<T>(1);
    static const int ione  = 1;
    static const int izero = 0;

    bool iscl = false;
    if (anrm > 0. && anrm < smlnum) {
        iscl  = true;
        scale = scalar<T>(calc_scale<Tr>(anrm, smlnum));
    } else if (anrm > bignum) {
        iscl  = true;
        scale = scalar<T>(calc_scale<Tr>(anrm, bignum));
    }

    if (iscl == 1) { multiply_inplace(arrA, abs(scale)); }

    int nru  = 0;
    int ncvt = 0;

    // Instead of copying U, S, VT, and A to the host and copying the results
    // back to the device, create a pointer that's mapped to device memory where
    // the computation can directly happen
    T *mappedA = static_cast<T *>(getQueue().enqueueMapBuffer(
        *arrA.get(), CL_FALSE, CL_MAP_READ, sizeof(T) * arrA.getOffset(),
        sizeof(T) * arrA.elements()));
    std::vector<T> tauq(min_mn), taup(min_mn);
    std::vector<T> work(lwork);
    Tr *mappedS0 = (Tr *)getQueue().enqueueMapBuffer(
        *arrS.get(), CL_TRUE, CL_MAP_WRITE, sizeof(Tr) * arrS.getOffset(),
        sizeof(Tr) * arrS.elements());
    std::vector<Tr> s1(min_mn - 1);