Function leastSquares

src/backend/opencl/solve.cpp:99–254

template<typename T>
Array<T> leastSquares(const Array<T> &a, const Array<T> &b) {
    int M  = a.dims()[0];
    int N  = a.dims()[1];
    int K  = b.dims()[1];
    int MN = min(M, N);

    Array<T> B = createEmptyArray<T>(dim4());
    gpu_blas_trsm_func<T> gpu_blas_trsm;

    cl_event event;
    cl_command_queue queue = getQueue()();

    if (M < N) {
#define UNMQR 0  // FIXME: UNMQR == 1 should be faster but does not work

        // Least squres for this case is solved using the following
        // solve(A, B) == matmul(Q, Xpad);
        // Where:
        // Xpad == pad(Xt, N - M, 1);
        // Xt   == tri_solve(R1, B);
        // R1   == R(seq(M), seq(M));
        // transpose(A) == matmul(Q, R);

        // QR is performed on the transpose of A
        Array<T> A = transpose<T>(a, true);

#if UNMQR
        const dim4 NullShape(0, 0, 0, 0);
        dim4 endPadding(N - b.dims()[0], K - b.dims()[1], 0, 0);
        B = (endPadding == NullShape
                 ? copyArray(b)
                 : padArrayBorders(b, NullShape, endPadding, AF_PAD_ZERO));
        B.resetDims(dim4(M, K));
#else
        B = copyArray<T>(b);
#endif

        int NB       = magma_get_geqrf_nb<T>(A.dims()[1]);
        int NUM      = (2 * MN + ((M + 31) / 32) * 32) * NB;
        Array<T> tmp = createEmptyArray<T>(dim4(NUM));

        vector<T> h_tau(MN);

        int info   = 0;
        Buffer *dA = A.get();
        Buffer *dT = tmp.get();
        Buffer *dB = B.get();

        magma_geqrf3_gpu<T>(A.dims()[0], A.dims()[1], (*dA)(), A.getOffset(),
                            A.strides()[1], &h_tau[0], (*dT)(), tmp.getOffset(),
                            getQueue()(), &info);

        A.resetDims(dim4(M, M));

        magmablas_swapdblk<T>(MN - 1, NB, (*dA)(), A.getOffset(),
                              A.strides()[1], 1, (*dT)(),
                              tmp.getOffset() + MN * NB, NB, 0, queue);