MATLAB and Fortran gives different results (precision)
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have 2 small programs in MATLAB and Fortran, respectively, which should be identical
MATLAB version
f=zeros(size(r1,1),3);
for ii=1:size(r1,1)
for iii=1:size(r2,1)
if all(r1(ii,:)==r2(iii,:))
continue
end
f(ii,:)=f(ii,:)+ScaFun(r1(ii,:),e1(ii),r2(iii,:),e2(iii));
end
end
function f=ScaFun(r1,e1,r2,e2)
r21 = r1-r2;
R = sqrt(r21(1)^2+r21(2)^2+r21(3)^2);
f = e1.*e2.*r21./(R^3);
end
Fortran:
function InvSq(r1,e1,r2,e2) result(f)
real*8, allocatable, intent(in) :: r1(:,:), e1(:), r2(:,:), e2(:)
integer*8 :: n1, n2, ii, jj
real*8, allocatable :: f(:,:)
n1 = size(r1,1)
n2 = size(r2,1)
allocate(f(n1,3),source = 0.0_8)
do ii = 1,n1
jj_loop: do jj = 1,n2
if (all(r1(ii,:)==r2(jj,:))) cycle jj_loop !skip self-to-self interaction (which InvSq_scalar() will return NaN)
f(ii,:) = f(ii,:)+InvSq_scalar(r1(ii,:),e1(ii),r2(jj,:),e2(jj))
end do jj_loop
end do
end function InvSq
function InvSq_scalar(r1,e1,r2,e2) result(f21)
real*8, intent(in) :: r1(3), e1, r2(3), e2
real*8 :: r21(3), R
real*8 :: f21(3)
r21 = r1-r2
R = sqrt(r21(1)**2+r21(2)**2+r21(3)**2);
f21 = e1*e2*r21/(R**3)
end function InvSq_scalar
However, for randomly generated inputs (e.g., sizein=1e3)
r1=2*rand(sizein,3)-1;
e1=2*rand(sizein,1)-1;
r2=2*rand(sizein,3)-1;
e2=2*rand(sizein,1)-1;
The (absolute) difference between MATLAB and Fortran results are on the order of 1e-13 for most, with a few as high as 1e-12.
Now, I understand that some of the inputs may be badly scaled (i.e., when r1 very close to r2) and this could be problematic for floating point arithmetic. However, they should still use the same ieee doubles and excute the same operations in the same order and should give identical results. Is there something else under the hood makes this happening?
0 commentaires
Réponses (1)
Walter Roberson
le 26 Déc 2022
FORTRAN and MATLAB do not necessarily use the same floating point rounding mode for calculations.
Control over rounding mode is not standardized in FORTRAN. See for example IBM https://www.ibm.com/docs/en/xcafbg/9.0.0?topic=SS3KZ4_9.0.0/com.ibm.xlf111.bg.doc/xlfopg/fpround.html versus Intel https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/language-reference/a-to-z-reference/h-to-i/ieee-set-rounding-mode.html
MATLAB has an undocumented function to set rounding behaviour (for code within MATLAB; might not apply to the high performance libraries that are automatically called.) See https://stackoverflow.com/questions/55682034/how-to-change-the-rounding-mode-for-floating-point-operations-in-matlab
6 commentaires
Walter Roberson
le 27 Déc 2022
In MATLAB, sum() over a larger array might not give the same result as looping adding one value at a time. sum() for sufficiently large arrays is passed to the high performance math libraries. The values to be totalled are partitioned, and one core is allocated per partion, and then the partial results are added. This will typically give different rounding results.
Voir également
Catégories
En savoir plus sur Fortran with MATLAB dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!