Подпрограммы BLAS dgemm, dgemv и ddot не работают со скалярами?

Question

У меня есть подпрограмма Fortran, которая использует подпрограммы BLAS dgemm, dgemv и ddot, которые вычисляют матрицу * matrix, матрицу * vector и вектор * vector.У меня есть 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)

Итак, мой актуальный вопрос в том, прав ли я в том, что подпрограммы этих BLAS не работают должным образом при m = 1, или просто что-то не так в моем коде?Может ли компилятор повлиять на это?Я использую gfortran.

Answer 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 файл, на который я ссылался ранее в своем ответе)