Le subroutine BLAS dgemm, dgemv e ddot non funzionano con gli scalari?
Domanda
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.
Soluzione
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)