Подпрограммы BLAS dgemm, dgemv и ddot не работают со скалярами?
Вопрос
У меня есть подпрограмма 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.
Решение
Предполагается, что процедуры 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
файл, на который я ссылался ранее в своем ответе)