Fix LAPACKE_?lacpy_work row-major triangular copy corrupting the untouched triangle (#729)

[meng004]

Problem

LAPACKE_?lacpy_work (s/d/c/z) with matrix_layout == LAPACK_ROW_MAJOR
and a triangular uplo ('U'/'L') corrupts the part of the destination B
that the operation is supposed to leave untouched, and reads uninitialized heap
memory in the process. As reported in #729, only the requested triangle of B
should change; instead the complementary region of B is overwritten with
garbage. The column-major path is correct.

Root cause

LAPACKE/src/lapacke_slacpy_work.c (and the d/c/z siblings) implement the
row-major branch by transposing through temporaries:

} else if( matrix_layout == LAPACK_ROW_MAJOR ) {
    ...
    a_t = LAPACKE_malloc(...);   /* full m*n temporary */
    b_t = LAPACKE_malloc(...);   /* full m*n temporary, NOT initialized */
    LAPACKE_sge_trans( matrix_layout, m, n, a, lda, a_t, lda_t );
    LAPACK_slacpy( &uplo, &m, &n, a_t, &lda_t, b_t, &ldb_t ); /* writes only the triangle of b_t */
    LAPACKE_sge_trans( LAPACK_COL_MAJOR, m, n, b_t, ldb_t, b, ldb ); /* transposes the WHOLE b_t back */
    ...
}

The kernel LAPACK_?lacpy with a triangular uplo writes only the requested
triangle of b_t, so the complementary region of b_t stays uninitialized. The
final LAPACKE_?ge_trans then transposes the entire m-by-n b_t back
into the caller's B — including that never-written region — overwriting the
bytes of B that uplo requires to be preserved with uninitialized heap
contents. (In the file as shipped this is
LAPACKE/src/lapacke_slacpy_work.c:76, the "Transpose output matrices" call.)

Fix

A row-major m-by-n matrix with leading dimension lda is bit-identical to a
column-major n-by-m matrix — its mathematical transpose. Copying triangle
uplo of the row-major matrix is therefore exactly copying the opposite
triangle (U↔L) of that transpose, with the m/n extents swapped. So the
row-major branch can call the Fortran kernel directly on a/b with uplo
flipped and m/n swapped:

} else if( matrix_layout == LAPACK_ROW_MAJOR ) {
    char uplo_t = uplo;
    /* leading-dimension checks unchanged */
    if( uplo == 'U' || uplo == 'u' ) uplo_t = 'L';
    else if( uplo == 'L' || uplo == 'l' ) uplo_t = 'U';
    LAPACK_slacpy( &uplo_t, &n, &m, a, &lda, b, &ldb );
    info = 0;
}

Like the already-correct column-major path, the kernel now writes exactly the
requested triangle of B and never touches the complementary region. This is
minimal and, in fact, net-negative in size: it removes the two LAPACKE_mallocs,
the two LAPACKE_?ge_trans calls, and the error-exit ladder rather than adding
special cases. The existing row-major leading-dimension checks (lda < n →
-6, ldb < n → -8) are preserved verbatim, the column-major and
invalid-layout branches are untouched, the public signatures/headers are
unchanged, and uplo values other than 'U'/'L' (e.g. 'A', a full copy)
are layout-agnostic and pass straight through.

Test evidence

Validated with a two-arm build in which both arms share the identical
liblapack/libblas; only the four LAPACKE *lacpy_work objects differ
(pristine vs. patched). The oracle copies triangle uplo of the same logical
matrix A into a sentinel-filled B once in LAPACK_COL_MAJOR (ground truth)
and once in LAPACK_ROW_MAJOR; the two logical results must be identical.

check	before (pristine)	after (patched)
row-vs-col logical result, uplo='U'	mismatch (maxdiff 9.0)	identical (0.0)
row-vs-col logical result, uplo='L'	mismatch (maxdiff 29.0)	identical (0.0)
copied triangle vs. A	correct (0.0)	correct (0.0)
complementary region vs. sentinel	corrupted (9.0 / 29.0)	preserved (0.0)
full-copy control uplo='A' (row vs col)	identical	identical
rectangular M=3,N=5, uplo='U'/'L'	mismatch (7.0 / 22.0)	correct

The copied-triangle content was already correct before the fix; the sole defect
is the complementary region, which the fix restores to being preserved. The
uplo='A' full-copy control and the column-major reference are identical on both
arms, confirming the change does not perturb the paths it should not.

How to reproduce

Minimal C reproducer (double precision; the same holds for s/c/z):

#include <lapacke.h>
#include <stdio.h>
int main(void) {
    const int n = 5;
    double A[25], Bcol[25], Brow[25];
    for (int i = 0; i < 25; i++) A[i] = i + 1;      /* same logical A */
    for (int i = 0; i < 25; i++) Bcol[i] = Brow[i] = -9.0; /* sentinel */

    /* col-major: copy upper triangle of A into B (ground truth) */
    LAPACKE_dlacpy_work(LAPACK_COL_MAJOR, 'U', n, n, A, n, Bcol, n);
    /* row-major: same logical operation on the same logical A */
    LAPACKE_dlacpy_work(LAPACK_ROW_MAJOR, 'U', n, n, A, n, Brow, n);

    /* Brow, read row-major, must equal Bcol read col-major, element for element.
       Before the fix the strictly-lower part of Brow is garbage (not -9). */
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            if (Brow[i*n + j] != Bcol[j*n + i])
                printf("mismatch at (%d,%d): row=%g col=%g\n",
                       i, j, Brow[i*n + j], Bcol[j*n + i]);
    return 0;
}

Build with -llapacke -llapack -lblas. Before the fix it prints mismatches in
the complementary (strictly-lower) triangle; after the fix it prints nothing.

---

[mgates3]

Other than perhaps shortening the comments, this looks good. Nice that it greatly simplified the code and eliminated malloc and copy! (Though I didn't personally validate the code, and LAPACKE unfortunately lacks testers. LAPACK++ tests only column-major.)

[mgates3]

I would put comments like this in the PR description, not in the code itself. It clutters up the code -- years from now, the old behavior is irrelevant.

Same comment in other precisions.

[langou]

I agree with @mgates3 that this is a great patch. (1) This is an elegant way to handle the row major triangular matrices and, (2) in addition, there was a bug so that fixes it, great.

Also, my preference (similar to @mgates') would be to leave out the lengthy comment.