Function lstsq

numpy/linalg/umath_linalg.cpp:3994–4066

template<typename typ>
static void
lstsq(char **args, npy_intp const *dimensions, npy_intp const *steps,
             void *NPY_UNUSED(func))
{
using ftyp = fortran_type_t<typ>;
using basetyp = basetype_t<typ>;
    GELSD_PARAMS_t<ftyp> params;
    int error_occurred = get_fp_invalid_and_clear();
    fortran_int n, m, nrhs;
    fortran_int excess;

    INIT_OUTER_LOOP_7

    m = (fortran_int)dimensions[0];
    n = (fortran_int)dimensions[1];
    nrhs = (fortran_int)dimensions[2];
    excess = m - n;

    if (init_gelsd(&params, m, n, nrhs, dispatch_scalar<ftyp>{})) {
        LINEARIZE_DATA_t a_in, b_in, x_out, s_out, r_out;

        init_linearize_data(&a_in, n, m, steps[1], steps[0]);
        init_linearize_data_ex(&b_in, nrhs, m, steps[3], steps[2], fortran_int_max(n, m));
        init_linearize_data_ex(&x_out, nrhs, n, steps[5], steps[4], fortran_int_max(n, m));
        init_linearize_data(&r_out, 1, nrhs, 1, steps[6]);
        init_linearize_data(&s_out, 1, fortran_int_min(n, m), 1, steps[7]);

        BEGIN_OUTER_LOOP_7
            int not_ok;
            linearize_matrix((typ*)params.A, (typ*)args[0], &a_in);
            linearize_matrix((typ*)params.B, (typ*)args[1], &b_in);
            params.RCOND = (basetyp*)args[2];
            not_ok = call_gelsd(&params);
            if (!not_ok) {
                delinearize_matrix((typ*)args[3], (typ*)params.B, &x_out);
                *(npy_int*) args[5] = params.RANK;
                delinearize_matrix((basetyp*)args[6], (basetyp*)params.S, &s_out);

                /* Note that linalg.lstsq discards this when excess == 0 */
                if (excess >= 0 && params.RANK == n) {
                    /* Compute the residuals as the square sum of each column */
                    int i;
                    char *resid = args[4];
                    ftyp *components = (ftyp *)params.B + n;
                    for (i = 0; i < nrhs; i++) {
                        ftyp *vector = components + i*m;
                        /* Numpy and fortran floating types are the same size,
                         * so this cast is safe */
                        basetyp abs = abs2((typ *)vector, excess,
dispatch_scalar<typ>{});
                        memcpy(
                            resid + i*r_out.column_strides,
                            &abs, sizeof(abs));
                    }
                }
                else {
                    /* Note that this is always discarded by linalg.lstsq */
                    nan_matrix((basetyp*)args[4], &r_out);