Function solve

src/backend/opencl/cpu/cpu_solve.cpp:249–301

template<typename T>
Array<T> solve(const Array<T> &a, const Array<T> &b,
               const af_mat_prop options) {
    if (options & AF_MAT_UPPER || options & AF_MAT_LOWER) {
        return triangleSolve<T>(a, b, options);
    }

#ifdef AF_USE_MKL_BATCH
    if (a.dims()[2] > 1 || a.dims()[3] > 1) {
        return generalSolveBatched(a, b, options);
    }
#endif

    const dim4 NullShape(0, 0, 0, 0);

    dim4 aDims = a.dims();
    int batchz = aDims[2];
    int batchw = aDims[3];

    int M = a.dims()[0];
    int N = a.dims()[1];
    int K = b.dims()[1];

    Array<T> A = copyArray<T>(a);
    dim4 endPadding(max(M, N) - b.dims()[0], K - b.dims()[1], 0, 0);
    Array<T> B = (endPadding == NullShape
                      ? copyArray(b)
                      : padArrayBorders(b, NullShape, endPadding, AF_PAD_ZERO));

    mapped_ptr<T> aPtr = A.getMappedPtr();
    mapped_ptr<T> bPtr = B.getMappedPtr();

    for (int i = 0; i < batchw; i++) {
        for (int j = 0; j < batchz; j++) {
            auto pA = aPtr.get() + A.strides()[2] * j + A.strides()[3] * i;
            auto pB = bPtr.get() + B.strides()[2] * j + B.strides()[3] * i;

            if (M == N) {
                std::vector<int> pivot(N);
                gesv_func<T>()(AF_LAPACK_COL_MAJOR, N, K, pA, A.strides()[1],
                               &pivot.front(), pB, B.strides()[1]);
            } else {
                int sM = a.strides()[1];
                int sN = a.strides()[2] / sM;

                gels_func<T>()(AF_LAPACK_COL_MAJOR, 'N', M, N, K, pA,
                               A.strides()[1], pB, max(sM, sN));
            }
        }
    }
    if (M != N) { B.resetDims(dim4(N, K, B.dims()[2], B.dims()[3])); }

    return B;
}

#define INSTANTIATE_SOLVE(T)                                                 \
    template Array<T> solve<T>(const Array<T> &a, const Array<T> &b,         \