1

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を使用しています。

4

1 に答える 1

1

BLASルーチンは、サイズ1のオブジェクトで正しく動作するはずです。コンパイラに依存することはないと思いますが、依存しているBLASの実装に依存する可能性があります（ただし、これはバグだと思います）実装）。NetlibにあるBLASのリファレンス（読み取り：ターゲットに最適化されていない）実装は、そのケースをうまく処理します。

サイズ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 に答える