Repository: scipy/scipy
File: scipy/linalg/src/_linalg_inv.hh (C++)

Code before:
    }

    lwork = (4*n > lwork ? 4*n : lwork);

    /*
     * Finally, we can start allocating memory.
     *
     * The key point is that LAPACK always operates on F-ordered arrays.
     * The memory strategy thus depends on the `overwrite_a` value.
     *
     * For `overwrite_a=False` (default), we:
     *   - allocate a temp buffer (`scratch` and `data` below) once
     *   - for each slice, we
     *       - copy-and-transpose the slice into the temp buffer,
     *       - feed the buffer to LAPACK
     *       - copy-and-transpose the result back to the C-ordered result array
     *
     * For `overwrite_a=True`, we assume that
     *   - `ret_data` may point to the same memory as `ap_Am` array
     *   - the caller had ensured that the input array is Fortran-ordered
     *   - the caller wants to get the result also Fortran ordered
     *
     * It's a caller's responsibility to make sure that these pre-conditions are met,
     * none of them is checked here.
     * Therefore, if `overwrite_a = True`, we skip the copy-and-transpose steps above,
     * and fill the output `ret_data` straight with the output from LAPACK.
     *
     */
    CBLAS_INT buf_size = overwrite_a ? lwork : 2*n*n + lwork;

    T* buffer = (T *)malloc(buf_size*sizeof(T));
    if (NULL == buffer) { info = -101; return (int)info; }

    T *data=NULL, *scratch=NULL, *work=NULL;
    if (overwrite_a) {
        // work in-place 
        data = ret_data;
        work = &buffer[0];
    }
    else {
        // Chop buffer into parts, one for data and one for work
        data = &buffer[0];
        scratch = &buffer[n*n];
        work = &buffer[2*n*n];
    }

    CBLAS_INT* ipiv = (CBLAS_INT *)malloc(n*sizeof(CBLAS_INT));
    if (ipiv == NULL) {
        free(buffer);
        info = -102;
        return (int)info;

Review comment (suggestion):
```suggestion
     * Therefore, if `overwrite_a = True`, we skip the copy-and-transpose steps above,
     * and `ret_data` will simply contain the result from the LAPACK call.
```

This statement might not be 100% accurate as the buffer is just passed into LAPACK and the result is just returned blindly while this might make it seem as if it was a delibarate choice to then copy the result into an F-contiguous buffer.

Code after:
    }

    lwork = (4*n > lwork ? 4*n : lwork);

    /*
     * Finally, we can start allocating memory.
     *
     * The key point is that LAPACK always operates on F-ordered arrays.
     * The memory strategy thus depends on the `overwrite_a` value.
     *
     * For `overwrite_a=False` (default), we:
     *   - allocate a temp buffer (`scratch` and `data` below) once
     *   - for each slice, we
     *       - copy-and-transpose the slice into the temp buffer,
     *       - feed the buffer to LAPACK
     *       - copy-and-transpose the result back to the C-ordered result array
     *
     * For `overwrite_a=True`, we assume that
     *   - `ret_data` may point to the same memory as `ap_Am` array
     *   - the caller had ensured that the input array is Fortran-ordered
     *   - the caller wants to get the result also Fortran ordered
     *
     * It's a caller's responsibility to make sure that these pre-conditions are met,
     * none of them is checked here.
     * Therefore, if `overwrite_a = True`, we skip the copy-and-transpose steps above,
     * and `ret_data` will simply contain the result from the LAPACK call.
     *
     */
    CBLAS_INT buf_size = overwrite_a ? lwork : 2*n*n + lwork;

    T* buffer = (T *)malloc(buf_size*sizeof(T));
    if (NULL == buffer) { info = -101; return (int)info; }

    T *data=NULL, *scratch=NULL, *work=NULL;
    if (overwrite_a) {
        // work in-place 
        data = ret_data;
        work = &buffer[0];
    }
    else {
        // Chop buffer into parts, one for data and one for work
        data = &buffer[0];
        scratch = &buffer[n*n];
        work = &buffer[2*n*n];
    }

    CBLAS_INT* ipiv = (CBLAS_INT *)malloc(n*sizeof(CBLAS_INT));
    if (ipiv == NULL) {
        free(buffer);
        info = -102;
        return (int)info;