复数场景下Lapacke实现QR分解时Q矩阵结果异常问题求助

我使用Lapacke库在C语言中实现复数数据的QR分解，参考社区公开代码编写了如下函数：

// Q - 输入：待分解矩阵 / 输出：Q矩阵
// R - 输出：R矩阵
// rows - 输入：Q矩阵的行数
// columns - 输入：Q矩阵的列数
// 始终满足 rows >= columns 条件
void QR(lapack_complex_double * Q, lapack_complex_double * R, size_t rows, size_t columns){
    size_t i;
    lapack_complex_double* tau = malloc(columns*sizeof(lapack_complex_double));
    LAPACKE_zgeqrf(LAPACK_ROW_MAJOR, (int) rows, (int) columns, Q, (int) columns, tau); // 返回压缩格式的Q、R
    // 拷贝上三角矩阵R（维度为columns x columns）
    for(i = 0; i < columns; ++i)
        memcpy(R+i*columns+i, Q+i*columns+i, (columns-i)*sizeof(lapack_complex_double));
    LAPACKE_zungqr(LAPACK_ROW_MAJOR, (int) rows, (int) columns, (int) columns, Q, (int) columns, tau); // 生成Q矩阵
    free(tau);
}

有其他开发者此前使用zungqr函数时遇到过同类问题，切换为zunmqr后问题解决。我本次测试得到的R矩阵和Mathematica的QRDecomposition函数、Python的numpy.linalg.qr函数计算结果完全一致，说明LAPACKE_zgeqrf运行正常，问题应该出在LAPACKE_zungqr的调用逻辑上。

我使用如下5×5复数矩阵作为测试输入：

1 + 3j  6 + 8j    11 + 13j  16 + 18j  21 + 23j  
2 + 4j  7 + 9j    12 + 14j  17 + 19j  22 + 24j  
3 + 5j  8 + 10j   13 + 15j  18 + 20j  23 + 25j  
4 + 6j  9 + 11j   14 + 16j  19 + 21j  24 + 26j  
5 + 7j  10 + 12j  15 + 17j  20 + 22j  25 + 27j

我的C代码输出的Q矩阵前3列和Python、Mathematica结果的误差在10^-11量级，吻合度极高，但后2列误差分别达到10^-3、10^-2量级，差异肉眼可见。随着矩阵尺寸增大，Q矩阵的差异会进一步扩大，通常只有前2-3列吻合度较好。希望有相关经验的开发者帮忙定位问题。

内容的提问来源于stack exchange，提问作者Konstantin