Les sous-programmes BLAS dgemm, dgemv et ddot ne fonctionnent pas avec les scalaires?

Question

J'ai un sous-programme Fortran qui utilise les sous-programmes de BLAS, dgemm, dgemv et ddot, qui calculent matrice * matrice, vecteur * et vecteur *. J'ai m * m matrices et m * 1 vecteurs. Dans certains cas, m = 1. Il semble que ces sous-programmes ne fonctionnent pas bien dans ces cas. Ils ne donnent pas d'erreur, mais il semble y avoir une certaine instabilité numérique dans les résultats. Je dois donc écrire quelque chose comme:

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)

Donc, ma question est la suivante: ai-je raison de dire que les sous-routines de BLAS ne fonctionnent pas correctement lorsque m = 1 ou existe-t-il simplement un problème dans mon code? Le compilateur peut-il affecter cela? J'utilise gfortran.

Answer 1

Les routines BLAS sont supposées se comporter correctement avec des objets de taille 1. Je ne pense pas que cela dépende du compilateur, mais cela pourrait dépendre de la mise en oeuvre de BLAS sur laquelle vous vous appuyez (bien que je considère cela comme une bug de la mise en oeuvre). Implémentation de référence de BLAS (à lire: non optimisée pour la cible), disponible sur Netlib , gère bien ce cas.

J'ai effectué des tests sur les tableaux de taille 1 et les tranches de taille 1 d'un tableau plus grand (comme dans votre propre code), et ils fonctionnent correctement:

 $ 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     

Choses que je considérerais pour déboguer ce problème plus loin:

• Vérifiez que la définition de votre ddot est correcte
• Substituez le BLAS de référence à votre optimisé pour vérifier s'il change quoi que ce soit (vous pouvez simplement compiler et lier dans le fichier ddot.f auquel j'ai lié plus tôt dans ma réponse)