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

Code before:
     *       - 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;
    }

    // {ge,po,tr}con need rwork or iwork
    void *irwork;
    if constexpr(type_traits<T>::is_complex) {
        irwork = malloc(3*n*sizeof(real_type));   // {po,tr}con need at least 3*n
    } else {
        irwork = malloc(n*sizeof(CBLAS_INT));
    }
    if (irwork == NULL) {
        free(buffer);
        free(ipiv);
        info = -102;

Review comment (suggestion):
Something I missed on my first pass; I completely agree with the logic here. However, this would only work for 2D inputs right? Makes the logic a bit tangled in the iteration over the slices later on. Perhaps the `data = ret_data` statement can be moved into that loop and by doing something along the lines of:
```c++
if (overwrite_a) {
    data = compute_slice_ptr(idx, ret_data, n, n, strides[ndim-2], strides[ndim-1]);
} else {
    --- copy-and-transpose ---
}
```

This makes it somewhat easier to gauge: now you need the context of the memory allocations to now that there was some branching there that already set the `data` variable correctly + you'd also need the context of it only working for 2D inputs as otherwise this logic would break down.

Code after:
     *       - 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;
    }

    // {ge,po,tr}con need rwork or iwork
    void *irwork;
    if constexpr(type_traits<T>::is_complex) {
        irwork = malloc(3*n*sizeof(real_type));   // {po,tr}con need at least 3*n
    } else {
        irwork = malloc(n*sizeof(CBLAS_INT));
    }
    if (irwork == NULL) {
        free(buffer);
        free(ipiv);
        info = -102;