Fortran's MKL dgemm gibt null Ergebnis

Im folgenden Fortran-Programm verwende ich Intels MKL-Bibliothek, um Matrixmultiplikationen unter Verwendung von dgemm durchzuführen. Anfangs habe ich die Subroutine matmul verwendet und korrekte Ergebnisse erhalten. Als ich matmul zu dgemm in der folgenden Schleife übersetzte, erhielt ich alle Nullvektoren anstelle der richtigen Ausgabe. Ich schätze Ihre Hilfe.Fortran's MKL dgemm gibt null Ergebnis

program spectral_norm  
implicit none 
! 
integer, parameter :: n = 5500, dp = kind(0.0d0) 
real(dp), allocatable :: A(:, :), u(:), v(:), Au(:), Av(:) 
integer :: i, j 

allocate(u(n), v(n), A(n, n), Au(n), Av(n)) 

do j = 1, n 
    do i = 1, n 
     A(i, j) = Ac(i, j) 
    end do 
end do 

u = 1 
do i = 1, 10 
    call dgemm('N','N', n, 1, n, 1.0, A, n, u, n, 0.0, Au, n) 
    call dgemm('N','N', n, 1, n, 1.0, Au, n, A, n, 0.0, v, n) 
    call dgemm('N','N', n, 1, n, 1.0, A, n, v, n, 0.0, Av, n) 
    call dgemm('N','N', n, 1, n, 1.0, Av, n, A, n, 0.0, u, n) 
    !v = matmul(matmul(A, u), A) 
    !u = matmul(matmul(A, v), A) 
end do 

write(*, "(f0.9)") sqrt(dot_product(u, v)/dot_product(v, v)) 

contains 

pure real(dp) function Ac(i, j) result(r) 
integer, intent(in) :: i, j 
r = 1._dp/((i+j-2) * (i+j-1)/2 + i) 
end function 

end program spectral_norm

Dies gibt NaN, während die korrekte Ausgabe von matmul1.274224153 ist.

Quelle

2017-06-27 AboAmmar

Verwenden MKL-Module. Tjey kann Ihnen helfen, viele Fehler zu identifizieren und auch einfache Fortran 90-Schnittstellen zu diesen Subroutinen zu enthalten. –

@VladimirF - Ich glaube nicht, dass es ein Problem mit fehlenden Modulen ist. Ich habe die 'include-Verzeichnisse' in VStudio richtig eingestellt und die Option 'MKL verwenden' gesetzt. Wenn ich versuche, zwei andere Matrizen zu multiplizieren, funktioniert es. Zum Beispiel gibt 'ruf dgemm (' t ',' n ', n, n, n, 1.0, A, n, A, n, 0.0, A, n)' die richtige Matrix. – AboAmmar

Ich denke, dass Sie 1.d0 und 0.d0 verwenden sollten, da Sie dgemm verwenden, die Argumente mit doppelter Genauigkeit erfordern. 'dgem ('N', 'N', n, 1, n, 1. d0, A, n, u, n, 0.d0, Au, n)' –

Nun, vielen Dank für Ihre Anregungen. Ich denke, ich habe die Ursache des Fehlers herausgefunden. Die Reihenfolge der Multiplikation wurde in zwei Fällen umgekehrt, es sollte stattdessen A * Au und A * Av gewesen sein. Dies liegt daran, A hat Bestellung n x n und beide Au ans Av haben Bestellung n x 1. Daher können wir Au * A oder Av * A aufgrund von nicht übereinstimmenden Dimensionen nicht multiplizieren. Ich habe die korrigierte Version unten gepostet.

program spectral_norm  
implicit none 
! 
integer, parameter :: n = 5500, dp = kind(0.d0) 
real(dp), allocatable :: A(:,:), u(:), v(:), Au(:), Av(:) 
integer :: i, j 

allocate(u(n), v(n), A(n, n), Au(n), Av(n)) 

do j = 1, n 
    do i = 1, n 
     A(i, j) = Ac(i, j) 
    end do 
end do 

u = 1 
do i = 1, 10 
    call dgemm('N', 'N', n, 1, n, 1._dp, A, n, u, n, 0._dp, Au, n) 
    call dgemm('T', 'N', n, 1, n, 1._dp, A, n, Au, n, 0._dp, v, n) 
    call dgemm('N', 'N', n, 1, n, 1._dp, A, n, v, n, 0._dp, Av, n) 
    call dgemm('T', 'N', n, 1, n, 1._dp, A, n, Av, n, 0._dp, u, n) 
end do 

write(*, "(f0.9)") sqrt(dot_product(u, v)/dot_product(v, v)) 

contains 

pure real(dp) function Ac(i, j) result(r) 
    integer, intent(in) :: i, j 
    r = 1._dp/((i+j-2) * (i+j-1)/2 + i) 
end function 

end program spectral_norm

Dies gibt die korrekten Ergebnisse:

1.274224153 
Elapsed time 0.5150000  seconds

Quelle

2017-06-27 08:57:40 AboAmmar

Fortran's MKL dgemm gibt null Ergebnis

Antwort

Verwandte Themen