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

У меня есть подпрограмма 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 не работают должным образом, когда m = 1, или в моем коде что-то не так? Может ли компилятор повлиять на это? Я использую гфортран.


gfortran matrix fortran scalar
person Community    schedule 17.06.2009    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Подпрограммы 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, на который я ссылался ранее в своем ответе)
person F'x    schedule 23.07.2009