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.

Foi útil?

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)
Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top