汇编角

当我开始调试时，我首先在 GOMP（OpenMP）中看到了崩溃。我重新编译了所有 R 和所需的软件包，并禁用了 OpenMP，幸运的是，该错误并没有消失。此示例

bdgraph.sim( n = 70, p = 5, size = 7, vis = TRUE )

导致错误

BLAS/LAPACK routine 'DPOTRS' gave error code -1

而它在较旧的 GFortran 中可以正常工作。我发现该错误是由于 LAPACK 函数 DPOTRS 中的 UPLO 已变为 0。因此，我在 UPLO 上添加了一个监视点，并转到该行

CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )

我不得不查看汇编代码来了解原因，gdb 给了我指令 (movq $0x1,0x70(%rsp))，并且我已阅读了 DPOSV 的完整且幸运的是很短的反汇编代码。这是有趣的部分

        CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
  1174d4:       48 8b 04 24             mov    (%rsp),%rax <======= rax holds LDB
  1174d8:       4c 89 7c 24 68          mov    %r15,0x68(%rsp) <=== save INFO to output param
  1174dd:       49 89 d8                mov    %rbx,%r8 <========== pass LDA as LDA
  1174e0:       4c 89 e1                mov    %r12,%rcx <========= pass A as A
  1174e3:       4c 8b 4c 24 08          mov    0x8(%rsp),%r9 <===== pass B as B
  1174e8:       4c 89 ea                mov    %r13,%rdx <========= pass NRHS as NRHS
  1174eb:       48 89 ee                mov    %rbp,%rsi <========= pass N as N
  1174ee:       4c 89 f7                mov    %r14,%rdi <========= pass UPLO as UPLO
  1174f1:       48 c7 44 24 70 01 00    movq   $0x1,0x70(%rsp) <=== pass 1 hidden arg on stack
  1174f8:       00 00 
  1174fa:       48 89 44 24 60          mov    %rax,0x60(%rsp) <=== pass LDB as LDB (stack)
      END
  1174ff:       48 83 c4 28             add    $0x28,%rsp <== remove 5 vars from stack
  117503:       5b                      pop    %rbx
  117504:       5d                      pop    %rbp
  117505:       41 5c                   pop    %r12
  117507:       41 5d                   pop    %r13
  117509:       41 5e                   pop    %r14
  11750b:       41 5f                   pop    %r15
         CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
  11750d:       e9 de 56 ef ff          jmpq   cbf0 <dpotrs_@plt> <=== tail call to dpotrs

位于 1174f1 的指令是破坏 BDgraph 软件包中的函数中 DPOSV 调用者的堆栈帧中 UPLO (uplo) 的移动。

void inverse( double A[], double A_inv[], int *p )
{
        int info, dim = *p;
        char uplo = 'U';

        // creating an identity matrix
        #pragma omp parallel for
        for( int i = 0; i < dim; i++ )
                for( int j = 0; j < dim; j++ )
                        A_inv[ j * dim + i ] = ( i == j );
  
        // LAPACK function: computes solution to A * X = B, where A is symmetric positive definite ma
        F77_NAME(dposv)( &uplo, &dim, &dim, A, &dim, A_inv, &dim, &info );
}

事实上，此移动还会破坏堆栈上的其他字节，顺便说一下，uplo 持有的字节从移动中获取了 7 个零字节中的一个。编译器将 1 写入堆栈，因为它想将其作为隐藏参数传递给 DPOTRS

SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )

DPOTRS 除了一个隐藏参数外，还接受 8 个上述参数，该隐藏参数是 UPLO 的长度，即 1。生成的代码不会调用 DPOTRS，而是修改已给定的调用参数到 DPOSV，然后跳转到 DPOTRS。这称为尾调用优化/消除（在本例中为同级调用优化），它节省了递归程序中的堆栈深度，这允许使用有限的堆栈运行深度递归代码，并且原则上可以提高性能。此 1 作为 DPOTRS 的隐藏参数应该已经传递给 DPOSV，但事实并非如此，因此堆栈上没有空间，因此堆栈损坏。

崩溃并未发生在旧版 GFortran 中，该版本不会对 DPOTRS 的调用进行尾调用优化，而是会为调用参数分配更多堆栈空间，并正常调用 DPOTRS。

Fortran 和 LAPACK 中的字符串

与 C 不同，Fortran 中的字符串不是以零结尾的。很明显，对于子例程的 CHARACTER*(*) 参数，Fortran 需要以某种方式将长度传达给被调用的函数。在某些系统中，长度是内存布局的一部分（例如 Pascal），但在 Fortran 中，它通常作为隐藏参数自动传递给函数，位于参数列表的末尾。该参数在 Fortran 程序中不可见。这些隐藏长度参数当前在 GFortran、ifort、flang 和可能的其他编译器中使用。但是，它们是特定于编译器的，例如，它们在 GFortran 7 和 GFortran8 之间会更改类型（在我的系统上，int 到 size_t，这也意味着 32 位到 64 位）。这导致 C 和 Fortran（以及不同的 Fortran 编译器）之间的互操作性问题，并且这个问题已经存在很多年了，BLAS/LAPACK/CBLAS/LAPACKE 通过不在语言之间传递此类字符串来避免此问题。

CHARACTER/CHARACTER*(1) 的子程序参数，即长度为 1 的字符串，在 BLAS 和 LAPACK 中得到了广泛使用。有时实际参数可能更长，例如 “U” 的 “Upper”，但只有第一个字符用作参数。字符串参数的实际长度永远不需要，也永远不会被访问。这就是该接口一直能够正常工作的原因（并且仍在使用其他编译器和 GFortran 的旧版本，以及 -fno-optimize-sibling-calls）：隐藏的长度参数位于参数列表的末尾，因此使用允许变参的通用调用约定时，它们在存在或不存在时都不会造成问题。现在唯一需要它们的情况是使用最新 GFortran 进行尾调用优化。

LAPACK C/Fortran 接口

BLAS/LAPACK 采用的做法是，在从 C 调用时，不会传递字符串长度参数（始终为 1）。然后，将相同的接口用于从 Fortran 调用和从 C 调用，甚至用于不同 Fortran 编译器之间的调用（可能并不总是有效，《编写 R 扩展》报告称，有时可能需要使用相同的编译器来构建 BLAS/LAPACK 和 R）。使用 BLAS/LAPACK 的应用程序通过该接口与这些库链接，而且通常是动态链接。R 也不例外：它可以使用包含的参考 BLAS/LAPACK 构建，但也可以使用外部实现，并且预计将在动态链接时做出选择：R 可以构建，并且通常会构建为，当不同的 BLAS/LAPACK 实现作为系统包安装时（例如，在 Linux 发行版中），它将由 R 使用（有关更多详细信息，请参阅 R 安装和管理）。

作为 参考 BLAS 和 LAPACK 的一部分，Netlib 还提供了 CBLAS（BLAS 的 C 接口）和 LAPACKE（LAPACK 的 C 接口）。这些是包装器，并且从底层 Fortran 代码获取略有不同的参数。但是，它们从 C 调用到 Fortran 时遵循相同的做法：它们不传递字符串长度。此外，应用程序通常直接调用 BLAS/LAPACK，而不是通过 CBLAS/LAPACKE，但是它们很难因遵循与 CBLAS/LAPACKE 相同的约定而受到指责，并且使用当前实现的这些接口并不能解决问题。

R 不使用 CBLAS/LAPACKE。BLAS 直接从 R 运行时和基本包中使用，而 LAPACK（仅）从基本包中的 stats 包中使用。R 包包括与 BLAS 和 LAPACK 直接交互的本机代码。如果最终解决此问题的方法是更改接口，R 将遵循套件，这意味着所有直接使用 BLAS/LAPACK 的 R 包都必须更新。

标准

该标准似乎并不能保证 BLAS/LAPACK 使用的做法必须得到 Fortran 编译器的支持。要获得具有该保证的 C 接口，可以使用 BIND(C)，并将长度为 1 的字符串参数声明为 CHARACTER (c_char)。Richard J. Hanson 和 Tim Hopkins 合著的“使用现代 Fortran 进行数值计算”一书有一个示例，该示例使用这些特性来告诉 Fortran 编译器生成一个可以从 C 调用且不需要将 1 作为字符串长度的函数 C_dgemm（DGEMM 来自 BLAS）。此包装函数 C_dgemm 调用由相同 Fortran 编译器编译的 Fortran 函数 dgemm_，但来自原始 Netlib 代码（示例 在此处提供）。通过此方法，可以在链中获得一个额外的函数调用，但在大多数情况下，这应该是可以忽略的开销。

原则上，Netlib 然后所有优化 BLAS/LAPACK 的实现者都可以切换到这种接口，然后所有应用程序（包括 R 和 R 包）都可以切换到这种接口。如果可以在不更改函数名称的情况下实现这一点，那就太好了（并且对于 1，情况就是这样，C_dgemm 不需要长度）。从 Fortran 到 Fortran 的调用（绕过那些 C 接口）将与现在一样，仍然处于不受标准支持的区域。

此解决方案无法在用户端（例如在 R/R 包中）实现，因为 C 接口需要由与库相同的 Fortran 编译器生成，因此无法在动态链接时在 BLAS/LAPACK 实现之间切换。

传递长度

从 C 调用 BLAS/LAPACK Fortran 代码时，也可以将 1 作为字符串长度传递，即使这是特定于编译器的，并且容易受到编译器更改以及当前根本不传递长度的做法的影响（例如，GFortran 7 和 GFortran 8 之间的类型发生了变化）。

如果在用户端采用此方法（GNU Octave 通过宏来实现），那么在 BLAS/LAPACK 实现之间切换时将再次出现问题。

最近 GFortran 9（主干）添加了一个选项 -fc-prototypes-external，用于生成 Fortran 函数的 C 原型，即包括所有隐藏参数及其正确大小（选项文档，gfortran 约定规范）。这很有用，但实际上取决于编译器，因此无法真正用于与 BLAS/LAPACK 交互的用户代码中。

检测问题

可以从 BDgraph 中提取示例，将其提取到独立示例中，以便在任何 R 代码之外演示问题

#include <stdio.h>

void dposv_(const char* uplo, const int* n, const int* nrhs,
            double* a, const int* lda,
            double* b, const int* ldb, int* info);

int main(int argc, char **argv) {
  char UPLO = 'U';  
  int N = 1;  
  int NRHS = 1;
  double A[1] = {2};
  int LDA = 1;
  double B[1] = {10};
  int LDB = 1;
  int INFO;

  dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO);
  printf("Result: %.0f Info: %d UPLO: %c\n", B[0], INFO, UPLO);
  return 0;
}

可以将其链接到最新 参考 LAPACK 的版本，该版本使用 GFortran 的默认选项（`-O2 -frecursive`）进行编译。我构建的示例如下

gfortran -O2 -g -o dposv dposv.c liblapack.a librefblas.a

输出应为

Result: 5 Info: 0

但在 Ubuntu 19.04 上的 GFortran 7、8、9 中，我得到

 ** On entry to DPOTRS parameter number  1 had an illegal value

注意：这不应用于检查已安装的编译器是否有问题。问题可能存在，但此程序可能会通过，因为堆栈的其他部分已损坏。我在 Debian 和 Fedora 上看到过这种情况。

一种更可靠的方法（但仅限于此 LAPACK 调用）是查看反汇编（objdump -dS）并搜索 dposv，是否具有上面显示的有问题的尾调用。此外，可以使用 GFortran 选项 -fdump-tree-all 从不同的编译器传递中获取输出，然后查看 dlapack.f.231t.optimized（可以使用 -fdump-tree-optimized 来获取它），搜索 dposv，然后检查是否存在（注意结尾处的 [tail call]）

dpotrs (uplo_18(D), n_22(D), nrhs_23(D), a_29(D), lda_14(D), b_31(D), ldb_15(D), info_16(D), 1); [tail call]

问题的当前范围

根据分析 R 中包含的参考 LAPACK 的 GFortran 9 编译输出，似乎有 38 个 LAPACK 函数可能会受到影响：dlarrc dpotrs dtrtrs dgetrs dlarzb dlatzm dlarfx ilaenv dtrtri dsygst dpotrf dlauum dpbtrf dlalsd dbdsdc dpbsv dpftrs dposv dpotri dppsv dpstrf dsbgv dspsv dsytri2 zdrscl zgesv zgetrf zgetrs zlaed0 zlahr2 zlalsd zlarfx zlatdf zlauum zpotrf zpotri ztrtri ztrtrs。这是采用字符参数并尾调用到另一个函数的函数列表。我通过解析 -fdump-tree-all 的输出进行了检查，但我没有检查它们是否都实际重新使用了隐藏的长度参数，以及该参数是否在堆栈中，这取决于平台。调用其中一些函数的 R 包数量超过一百个，尽管实际上只有大约 20 个被发现测试失败。

有问题的 GFortran 版本已经出现在 Linux 发行版中。

在 Ubuntu 19.04 中，所有 GFortran 7、8、9 都已经出现此问题（7.4.0-8ubuntu1、8.3.0-6ubuntu1、9-20190402-1ubuntu1）。GFortran 6 未受影响，运行良好。GFortran 8 为默认版本。在 Debian Sid 中，GFortran 7 和 GFortran 8 已经出现问题（“Debian 8.3.0-7”、“Debian 7.4.0-9”）。在 Debian Buster 中，GFortran 8（“Debian 8.3.0-6”）已经出现问题。在 Debian Stretch 中，GFortran 6（“Debian 6.3.0-18+deb9u1”）运行良好。在 Fedora 30 中，GFortran 9（“Red Hat 9.1.1-1”）已经出现问题。我只查看了 x86_64。

发行版已经包含出现此问题的 GFortran 版本，但这并不意味着这些发行版中的 LAPACK/BLACK 二进制包也会出现问题，而且我尚未发现任何出现此问题的二进制包，但我只查看了 dposv，这可能不够，如下所示，而且如果重新构建二进制文件，情况也可能发生变化。

例如，在 Ubuntu 19.04 中，liblapack-dev 3.8.0-2、libopenblas-dev 0.3.5+ds-2 和 libatlas3-base 3.10.3-8 仍然运行良好（Debian Sid 中的包也是如此）。但是，在使用默认编译器（apt-get build-dep、apt-get source、dpkg-buildpackage）重新构建时，它们肯定或很可能出现问题。liblapack-dev 仅在重新构建时才会出现错误的 dposv。OpenBLAS 和 ATLAS 使用 liblapack-pic 作为 LAPACK 代码，因此问题仅在重新构建 liblapack-pic、安装它，然后重新构建 libopenblas-dev 或 libatlas-base-dev 后才会出现。使用 OpenBLAS，然后会获得错误的 dposv。使用 ATLAS，不会出现错误，但这是一种巧合，因为 ATLAS 重新实现了选定的 LAPACK 函数，包括 dposv，然后以不同的方式编译（而且我没有详细检查其反汇编），但其他函数也取自 LAPACK（通过 liblapack-pic）。

ATLAS 在内部使用包装器从 C 正确调用 Fortran，包括长度，但仅对作者需要的几个例程使用，并且不会影响与 LAPACK 的外部接口。请注意，ATLAS 对其构建的 Fortran 源进行此操作，并且不会在链接时进行替换，因此在构建 C 包装器时，Fortran 编译器是已知的且固定的。

下一步

GFortran 开发人员正在积极解决此问题，以便 BLAS/LAPACK 仍然可以与 GFortran 配合使用，并且错误报告讨论中提供了更多更新。毕竟，一个无法与 BLAS/LAPACK 配合使用的 Fortran 编译器（抽象自其错误）几乎没有用处。

目前，每个人都应该对 GFortran 7 及更高版本使用 -fno-optimize-sibling-calls。这已经是 R-Devel 和 R-Patched 中的默认设置，但遗憾的是，在使用 GFortran 7 及更高版本构建 R 的发布版本时必须显式添加它（Debian Sid 已经有了使用此选项的 R 3.6.0-2）。不过，此设置仅修复了 R 中包含的参考 BLAS/LAPACK。

使用 -fno-optimize-sibling-calls 后，CRAN 团队发现的所有因编译器更改而导致的新包故障都已重新开始工作（这意味着，它们已使用 R、包以及最重要的是使用 -fno-optimize-sibling-calls 重新编译的 BLAS/LAPACK 进行了检查）。为了安全起见，我认为使用此选项编译所有 Fortran 代码是有意义的，并且必须对所有 BLAS/LAPACK 实现使用此选项。显然，Linux 发行版需要做很多工作才能确保在构建二进制包时使用此选项。

我认为 BLAS/LAPACK 的接口需要在参考实现中进行修复/补充，然后在经过优化的实现中进行修复/补充，然后在所有应用程序中进行修复/补充（除非作者找到一种方法，让接口保持现状，但保证符合标准）。尝试使用不同的方法在 R 包中单独修复此问题将出现我上面提到的问题，并且可能会导致混乱。

但是，使用 C 和 Fortran 混合编写的 R 包不应依赖 Fortran 编译器的非标准互操作行为（编写 R 扩展不鼓励使用字符串，即使长度为 1），除非调用 BLAS/LAPACK。

有时可以使用 LTO（链接时优化）由编译器/链接器检测到 C 和 Fortran 调用之间（以及不同源文件中的 C 声明和定义之间）的不兼容问题。在 GCC 中，它是 -Wlto-type-mismatch 警告。由于缺少隐藏的字符串长度参数，在调用 BLAS/LAPACK 时必须预料到这些警告，但如果可能，应修复所有其他情况，并鼓励包作者尝试（请参阅 R 的 --enable-lto 配置选项）。请注意，在使用 LTO 构建参考 LAPACK 分发版时，会针对 CBLAS 和 LAPACKE 获得这些警告（即使在 GFortran 更改之后，这些不兼容性不再无害，情况也是如此）。我们还在查看针对 R 本身发出的这些警告。

对于好奇的人来说，可以在 GFortran 开发人员的讨论 GCC PR#90329 中找到有趣的指针和更多最新信息，以及他们针对参考 lapack 打开的 问题。