BLAS sub-rotinas dgemm, dgemv e ddot não funciona com escalares?
Pergunta
Eu tenho uma sub-rotina Fortran que usa sub-rotinas 'BLAS dgemm, dgemv e ddot, que matriz de calcular * matriz, matriz * vector e vector * vetor. I ter de m * m * m as matrizes e vectores 1. Em alguns casos m = 1. Parece que esses sub-rotinas não funciona bem nesses casos. Eles não dá erros, mas parece haver alguma instabilidade numérica em resultados. Então eu tenho que escrever 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)
Então, minha pergunta real é, estou certo de que sub-rotinas desses BLAS não funciona adequadamente quando m = 1 ou há simplesmente errado alguma coisa no meu código? Pode o compilador afetar isso? Estou usando gfortran.
Solução
BLAS rotinas devem se comportar corretamente com objetos de tamanho 1. Eu não acho que isso pode depender de compilador, mas poderia possível dependem da implementação de Blas Você está confiando em (embora eu considerá-lo um bug da implementação). A referência (leia-se: não-alvo-optimizado) aplicação de BLAS, que pode ser encontrado no Netlib , alças que caso bem.
Eu fiz alguns testes em ambas as matrizes de tamanho 1, eo tamanho-1 fatias de matriz maior (como em seu próprio código), e ambos funcionam bem:
$ 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
As coisas que eu considero para depurar este problema ainda mais:
- Verifique se a sua definição
ddot
está correto - Substitua os BLAS referência ao seu otimizado um, para verificar se ela muda nada (você pode apenas compilar e link no arquivo
ddot.f
I ligada ao anteriormente na minha resposta)