Le subroutine BLAS dgemm, dgemv e ddot non funzionano con gli scalari?

Question

Ho una subroutine Fortran che usa le subroutine BLAS dgemm, dgemv e ddot, che calcolano matrice * matrice, matrice * vettore e vettore * vettore. Ho m * m matrici e m * 1 vettori. In alcuni casi m = 1. Sembra che quelle subroutine non funzionino bene in quei casi. Non danno errori, ma sembra esserci una certa instabilità numerica nei risultati. Quindi devo scrivere qualcosa del tipo:

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)

Quindi la mia vera domanda è: ho ragione che quelle subroutine BLAS non funzionano correttamente quando m = 1 o c'è qualcosa di sbagliato nel mio codice? Il compilatore può influire su questo? Sto usando gfortran.

Answer 1

Si suppone che le routine BLAS si comportino correttamente con oggetti di dimensione 1. Non penso che possa dipendere dal compilatore, ma potrebbe dipendere dall'implementazione di BLAS su cui si fa affidamento (anche se lo considererei un bug dell'implementazione). L'implementazione di riferimento (leggi: non ottimizzata per il target) di BLAS, che può essere trovata su Netlib , gestisce bene il caso.

Ho eseguito alcuni test su entrambi gli array di dimensioni 1 e dimensioni 1 di array più grandi (come nel tuo codice) ed entrambi funzionano bene:

 $ 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     

Cose che prenderei in considerazione per debug ulteriormente questo problema:

• Verifica che la tua definizione ddot sia corretta
• Sostituisci il BLAS di riferimento con quello ottimizzato, per verificare se cambia qualcosa (puoi semplicemente compilare e collegare il file ddot.f che ho collegato in precedenza nella mia risposta)