BLASサブルーチンdgemm、dgemv、ddotはスカラーで動作しませんか？

Question

BLASのサブルーチンdgemm、dgemvおよびddotを使用するFortranサブルーチンがあり、マトリックス*マトリックス、マトリックス*ベクトル、およびベクトル*ベクトルを計算します。 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を使用しています。

Answer 1

BLASルーチンは、サイズ1のオブジェクトで正しく動作するはずです。コンパイラに依存することはできないと思いますが、依存している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を最適化されたBLASに置き換えて、それが何かを変更するかどうかを確認します（先ほど答えたリンク先の ddot.f ファイルをコンパイルしてリンクできます）