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

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

1 голос
/ 17 июня 2009

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

1 Ответ

1 голос
/ 23 июля 2009

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