Las subrutinas BLAS dgemm, dgemv y ddot no funcionan con escalares?

Question

Tengo una subrutina Fortran que usa las subrutinas BLAS dgemm, dgemv y ddot, que calculan matriz * matriz, matriz * vector y vector * vector. Tengo m * m matrices y m * 1 vectores. En algunos casos m = 1. Parece que esas subrutinas no funcionan bien en esos casos. No dan errores, pero parece haber cierta inestabilidad numérica en los resultados. Entonces tengo que escribir algo como:

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)

Entonces, mi pregunta real es: ¿tengo razón en que esas subrutinas de BLAS no funcionan correctamente cuando m = 1 o hay algo que está mal en mi código? ¿Puede el compilador afectar esto? Estoy usando gfortran.

Answer 1

Se supone que las rutinas BLAS se comportan correctamente con objetos de tamaño 1. No creo que pueda depender del compilador, pero podría depender de la implementación de BLAS en la que confía (aunque lo consideraría un error de la implementación). La implementación de referencia (léase: no optimizada para el objetivo) de BLAS, que se puede encontrar en Netlib , maneja ese caso bien.

He realizado algunas pruebas tanto en matrices de tamaño 1 como en porciones de tamaño 1 de matriz más grande (como en su propio código), y ambas funcionan bien:

 $ 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     

Cosas que consideraría para depurar este problema aún más:

• Compruebe que su definición de ddot sea correcta
• Sustituya la referencia BLAS por la optimizada, para verificar si cambia algo (puede compilar y vincular en el archivo ddot.f al que he vinculado anteriormente en mi respuesta)