在Fortran90中BLAS函数返回零

Question

我正在学习如何在Fortran90中使用BLAS，并使用子程序SAXPY和函数SNRM2编写了一个简单的程序。程序计算两点之间的距离，从另一个向量减去一个向量，然后取结果的欧几里德范数。

我将SNRM2的返回值指定为external，根据对类似问题“调用BLAS函数”的回答。我的完整计划：

program test
implicit none

real :: dist
real, dimension(3) :: a, b
real, external :: SNRM2

a = (/ 3.0, 0.0, 0.0 /)
b = (/ 0.0, 4.0, 0.0 /)

call SAXPY(3, -1.0, a,1, b,1)
print *, 'difference vector: ', b

dist = 6.66  !to show that SNRM2 is doing something
dist = SNRM2(3, b, 1) 
print *, 'length of diff vector: ', dist

end program test

该方案的结果是：

difference vector:   -3.00000000       4.00000000       0.00000000    
length of diff vector:    0.00000000

差分向量是正确的，但是长度应该是5。那么为什么SNRM2要返回一个零的值呢？

我知道变量dist是由SNRM2修改的，所以我不怀疑我的openBLAS安装是坏的。我正在运行macos10.13，并安装了所有的自制软件。

我正在使用gfortran进行编译，启用了许多标志，并且没有收到警告：

gfortran test.f90 -lblas -g -fimplicit-none -fcheck=all -fwhole-file -fcheck=all -fbacktrace -Wall -Wextra -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wconversion -Wno-unused-parameter -pedantic -o test

我试着看snrm2.f代码，但我没有看到任何潜在的问题。

我还尝试用real(4)或real(selected_real_kind(6))声明变量，但行为没有变化。

谢谢!

Answer 1

根据这个页面，苹果加速框架附带的BLAS中的单个精确例程似乎存在一些问题。在我的Mac ( On 10.11)上，gfortran-8.1 (通过Homebrew安装)+默认BLAS (在系统中)给出了错误的结果：

$ gfortran-8 test.f90 -lblas
or
$ gfortran-8 test.f90 -L/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Versions/Current/ -lBLAS
$ ./a.out
 difference vector:   -3.00000000       4.00000000       0.00000000    
 length of diff vector:    0.00000000  

当显式地链接到OpenBLAS时(通过Homebrew安装)提供了正确的结果：

$ gfortran-8 test.f90 -L/usr/local/Cellar/openblas/0.2.20_2/lib -lblas
$ ./a.out
 difference vector:   -3.00000000       4.00000000       0.00000000    
 length of diff vector:    5.00000000 

上面的页面表明，当以一种与旧的g77样式不兼容的方式链接到系统BLAS时，就会出现问题。实际上，附加-ff2c选项提供了正确的结果：

$ gfortran-8 -ff2c test.f90 -lblas
$ ./a.out
 difference vector:   -3.00000000       4.00000000       0.00000000    
 length of diff vector:    5.00000000  

但是我想最好使用最新的OpenBLAS (比使用-ff2c选项).

下面是C中的单独测试(以检查问题不特定于gfortran)。

// test.c
#include <stdio.h>
float snrm2_( int*, float*, int* );

int main()
{
    float b[3] = { -3.0f, 4.0f, 0.0f };
    int n = 3, inc = 1;

    float dist = snrm2_( &n, b, &inc );

    printf( "b = %10.7f %10.7f %10.7f\n", b[0], b[1], b[2] );
    printf( "dist = %10.7f\n", dist );
    return 0;
}

$ gcc-8 test.c -lblas
$ ./a.out
b = -3.0000000  4.0000000  0.0000000
dist =  0.0000000

$ gcc-8 test.c -lblas -L/usr/local/Cellar/openblas/0.2.20_2/lib
$ ./a.out
b = -3.0000000  4.0000000  0.0000000
dist =  5.0000000

据我所尝试，双精度版本(DNRM2)甚至适用于系统BLAS，因此问题似乎只适用于单精度版本(如上面页面中所建议的)。