0

I have a piece of code written in Fortran and in Matlab. They do exactly the same calculation, namely

  1. Construct a tanh -field and find its Laplacian
  2. Multiply some terms together

The result of this multiplication yields a matrix, whose (4,4)th and (6,6)th I subtract.

  • In Fortran their difference is ~1e-20
  • In Matlab their difference is identically zero.

This issue is very critical, as I test if this number is less than zero. Question: Is there a way to perform the calculations such that I get the same precision in Matlab as in Fortran?

I list the codes below:

MatLAB

clear all

weights = [4./9, 1./9,1./9,1./9,1./9, 1./36,1./36,1./36,1./36];
dir_x   = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
dir_y   = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTANTS
length_y = 11; length_x = length_y;
y_center = 5; x_center  = y_center;


densityHigh = 1.0;
densityLow  = 0.1;
radius  = 3.0;
c_width = 1.0;

average_density = 0.5*(densityHigh+densityLow);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





for x=1:length_x
    for y=1:length_y
        for i=1:9
            fIn(i, x, y) = weights(i)*densityHigh;
            test_radius = sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center));
            if(test_radius <= (radius+c_width))
                fIn(i, x, y) = weights(i)*( average_density - 0.5*(densityHigh-densityLow)*tanh(2.0*(radius-sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center))/c_width)) );
            end
        end
    end
end 

ref_density_2d = ones(length_x)*average_density;
for i=1:length_x
    ref_density(:,:,i) = abs(ref_density_2d(:, i)');
end


          rho = sum(fIn);
laplacian_rho = (+1.0*(circshift(rho(1,:,:), [0, -1, -1]) + circshift(rho(1,:,:), [0, +1, -1]) + circshift(rho(1,:,:), [0, -1, +1]) + circshift(rho(1,:,:), [0, +1, +1])) + ...
                 +4.0*(circshift(rho(1,:,:), [0, -1, +0]) + circshift(rho(1,:,:), [0, +1, +0]) + circshift(rho(1,:,:), [0, +0, -1]) + circshift(rho(1,:,:), [0, +0, +1])) + ...
                -20.0*rho(1,:,:));

psi   = 4.0*0.001828989483310*(rho-densityLow).*(rho-densityHigh).*(rho-ref_density) - laplacian_rho*(1.851851851851852e-04)/6.0;

psi(1,4,4)-psi(1,6,6)

Fortran

 PROGRAM main

 IMPLICIT NONE

 INTEGER, PARAMETER :: DBL = KIND(1.D0)
 REAL(KIND = DBL), DIMENSION(1:11,1:11) :: psi, rho
 INTEGER :: i, j, m, ie, iw, jn, js
 REAL(KIND = DBL) :: R, rhon, lapRho

 INTEGER, DIMENSION(1:11,1:11,1:4) :: ni

 REAL(KIND = DBL) :: kappa, kappa_6, kappa_12, kappaEf, beta, beta4



 beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))
 kappa    = 1.5D0*0.0001*1.D0/( (1.0 - 0.1)**2 ) 


!-------- Define near neighbours and initialize the density rho ----------------
 DO j = 1, 11
   DO i = 1, 11

! Initialize density
      rho(i,j) = 1.D0
        R =  DSQRT( ( DBLE(i)-5.0 )**2 + ( DBLE(j)-5.0 )**2 )
        IF (R <= (DBLE(3.0) + 1.D0)) THEN
          rho(i,j) = 0.55D0 - 0.5*0.9*TANH(2.D0*(DBLE(3.0) - R)/1.D0)
        END IF

 !Generate neighbors array
      ni(i,j,1) = i + 1
      ni(i,j,2) = j + 1
      ni(i,j,3) = i - 1
      ni(i,j,4) = j - 1
   END DO
 END DO


! Fix neighbours at edges
 ni(1,:,3) = 11
 ni(11,:,1) = 1
 ni(:,1,4) = 11
 ni(:,11,2) = 1



!--------- Differential terms for the stress form of the interfacial force -----
 DO j = 1, 11
   DO i = 1, 11

! Identify neighbors
     ie = ni(i,j,1)
     jn = ni(i,j,2)
     iw = ni(i,j,3)
     js = ni(i,j,4)

! Laplacian of the density rho
     lapRho = 4.D0*( rho(ie,j ) + rho(iw,j ) + rho(i ,jn) + rho(i ,js) )       &
             + rho(ie,jn) + rho(ie,js) + rho(iw,jn) + rho(iw,js) - 20.D0*rho(i,j)

! Define the chemical potential Psi
     psi(i,j) = 4.D0*beta*( rho(i,j) - 0.55 )*( rho(i,j) - 0.1 )*( rho(i,j) - 1.0 ) &
              - kappa*lapRho/6.D0
   END DO
 END DO


write(*,*) psi(6,6)-psi(4,4)


 END PROGRAM
Alexander Vogt
  • 17,879
  • 13
  • 52
  • 68
BillyJean
  • 1,537
  • 1
  • 22
  • 39
  • 6
    Despite having been advised otherwise in the answer to another version of this question (http://stackoverflow.com/questions/29867870/precision-of-calculations) you are still using single-precision numbers in your Fortran code. That is `0.1` is, by default, a single-precision real number in Fortran (unless you have used compiler options to adjust the usual behaviour and I bet you haven't, have you you little rascal ?). That might account for the differences you observe. By this stage in your quest you should have, at least, investigated this possibility. – High Performance Mark Apr 27 '15 at 15:54
  • I changed it so it uses double precision, but the conclusion is the same -- the results still differ. I spent a long time investigating, trust me... – BillyJean Apr 27 '15 at 15:58
  • 3
    I'm sorry, but no one is likely to go through that much code (in two languages), run it, and debug it. Narrow down your problem. Find out the first place where your results diverge. In addition to what @HighPerformanceMark said above, see my comment to your previous question as well – you need check that functions like `abs`, `sqrt`, `tanh`, etc. (and even division) return the same result in all cases. And by "check" I don't mean simply print out to a few decimal places. – horchler Apr 27 '15 at 15:58
  • @horchler I know, I can't expect anyone to do that... but seriously, it is so frustrating that I can't make it work. I also tried following up on the `tanh`-thing, and they differ on the 19th decimal or so. But there's no way for me to remedy this... ugh – BillyJean Apr 27 '15 at 16:01
  • This time I at least provided the Fortran-code for interested readers. – BillyJean Apr 27 '15 at 16:02
  • But then that basically means that I have no way to compare the `tanh`-functions etc. – BillyJean Apr 27 '15 at 16:05
  • @HighPerformanceMark: actually it's [supposedly 15–17 places](http://en.wikipedia.org/wiki/Double-precision_floating-point_format). I don't know how FORTRAN's `write` works, but perhaps the OP is not using the best means of displaying values? Something like `sprintf('%.17f',eps)` instead of `sprintf('%.17g',eps)`? – horchler Apr 27 '15 at 16:08
  • Evaluation of the two constants at the beginning of the Fortran code, i.e., 12.0*0.0001/(1.0*( (1.0 - 0.1)**4 )) and 1.5*0.0001*1.0/( (1.0 - 0.1)**2 ), using Python2.6.1 gives 0.0018289894833104711 and 0.0001851851851851852. I guess these numbers are used manually at the end of Matlab code. But the first figure seems to be written only 13 significant digits while the latter with 15 digits in Matlab code. Although I don't think a difference in 19th dp is reproducible with double precision, the above literal figures might slightly increase the difference in the output? – roygvib Apr 27 '15 at 16:27
  • @HighPerformanceMark has already pointed it out, sorry... – roygvib Apr 27 '15 at 16:30
  • I would rather like to see the Fortran code in the other question... – Alexander Vogt Apr 27 '15 at 16:35
  • 2
    I retracted my closure of this question, and voted to close the other one. @AlexanderVogt is right that his answer, below, is the right answer but belongs with this question not the other one which lacks Fortran code. – High Performance Mark Apr 27 '15 at 17:47

1 Answers1

6

You are still not consequently using double precision throughout your code, e.g.:

beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))

and many more. If I force the compiler to use double precision as default for floats (for gfortran the compile option is -fdefault-real-8), the result from your code is:

0.00000000000000000000000000000000000

So you need to fix your code. The cited line, for instance, should read: 

beta     = 12.D0*0.0001D0/(1.D0*( (1.0D0 - 0.1D0)**4 ))

[Although I despise the notation D0, but that's a different story]

Alexander Vogt
  • 17,879
  • 13
  • 52
  • 68