1

我有一个 Fortran 子程序，它使用 BLAS 的子程序 dgemm、dgemv 和 ddot，它们计算矩阵 * 矩阵、矩阵 * 向量和向量 * 向量。我有 m * m 矩阵和 m * 1 向量。在某些情况下，m=1。在这些情况下，这些子程序似乎不能很好地工作。他们没有给出错误，但结果似乎存在一些数值上的不稳定性。所以我必须写一些类似的东西：

if(m>1) then 
  vtuni(i,t) = yt(i,t) - ct(i,t) - ddot(m, zt(i,1:m,(t-1)*tvar(3)+1), 1, arec, 1)
else 
   vtuni(i,t) = yt(i,t) - ct(i,t) - zt(i,1,(t-1)*tvar(3)+1)*arec(1)

所以我的实际问题是，当 m=1 时，那些 BLAS 的子例程不能正常工作，或者我的代码中是否有问题，我是对的吗？编译器会影响这个吗？我正在使用 gfortran。

4

1 回答 1

1

BLAS 例程应该与大小为 1 的对象正确运行。我认为它不能依赖于编译器，但它可能依赖于你所依赖的 BLAS 的实现（尽管我认为它是执行）。BLAS 的参考（阅读：非目标优化）实现可以在Netlib上找到，可以很好地处理这种情况。

我已经对大小为 1 的数组和大小为 1 的更大数组的切片（如您自己的代码）进行了一些测试，它们都工作正常：

 $ cat a.f90 
 implicit none
 double precision :: u(1), v(1)
 double precision, external :: ddot
 u(:) = 2
 v(:) = 3
 print *, ddot(1, u, 1, v, 1)
 end
 $ gfortran a.f90 -lblas && ./a.out
  6.0000000000000000     

 $ cat b.f90                       
 implicit none
 double precision, allocatable :: u(:,:,:), v(:)
 double precision, external :: ddot
 integer :: i, j
 allocate(u(3,1,3),v(1))
 u(:,:,:) = 2
 v(:) = 3
 i = 2
 j = 2
 print *, ddot(1, u(i,1:1,j), 1, v, 1)
 end
 $ gfortran b.f90 -lblas && ./a.out
  6.0000000000000000

我会考虑进一步调试这个问题的事情：

  • 检查您的ddot定义是否正确
  • 将参考 BLAS 替换为您的优化版本，以检查它是否更改了任何内容（您可以编译并链接到ddot.f我之前在我的答案中链接到的文件）
于 2009-07-23T08:17:34.047 回答