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

Question

У меня есть подпрограмма 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) 

Так что мой фактический вопрос, я прав, что те, BLAS»Подпрограммы не работает должным образом, когда т = 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, с которым я связался ранее в своем ответе)