Fortran to Matlab confused to if statement

Question

I have to convert a fortran program to matlab and I'm facing a problem. Although the fortran results are correct when I run the matlab script stucks. I believe the problem is in the first IF statement. Am I missing something to my matlab conversion ? Thanks in advance.

This is the full program with results both in Fortran and Matlab. The latest version is @Rotem 's solution. I, also, tried to add the fix to f1 and f2 as @BenBarrowes mentioned but the program stucks again. Thanks in advance :)

Fortran version

program asxm_2
implicit none
real(8) a,b, par1, sum,sum1,sum2, x,x1,x2, h, f,fa,fb,f1,f2,par2,e
integer n,j,y,k
real(8), allocatable, dimension (:) :: partitionS, valueS, errorS
a=0.
b=5.+85
par1=100*(atan(b)-atan(a)) 
fa=100/(1+a**2) 
fb=100/(1+b**2) 

print*, 'METHOD SIMPSON'
do n=1,1000000
    h=(b-a)/n
    sum1=0. 
    sum2=0. 
    x1=a 
    x2=a 
    do j=1,n-1 
        k=j 
        if(k/2/=j/2.) then 
            if(j==1) x1=x1+h 
            if(j>1) x1=x1+2*h 

            f1=100/(1+x1**2)
            sum1=sum1+f1
        else
            x2=x2+2*h 
            f2=100/(1+x2**2) 
            sum2=sum2+f2
        endif
    enddo
    par2= (h/3)*(fa+4*sum1+2*sum2+fb) 
    e=par1-par2
    if(abs(e)<=0.001) exit
enddo

y=n 
allocate(partitionS(y),valueS(y), errorS(y))
do n=1,y
     h=(b-a)/n
     sum1=0.
     sum2=0.
     x1=a
     x2=a
     do j=1,n-1
         k=j
         if(k/2==j/2.) then
             x2=x2+2*h
             f2=100/(1+x2**2)
             sum2=sum2+f2
         else
             if(j==1) x1=x1+h
             if(j>1) x1=x1+2*h
             f1=100/(1+x1**2)
             sum1=sum1+f1
         endif
     enddo
     partitionS(n)=n 
     valueS(n)= (h/3)*(fa+4*sum1+2*sum2+fb) 
     errorS(n)=par1-valueS(n)
 enddo
 print*, 'Below are the results'
 print*, partitionS(y), valueS(y), errorS(y)
 deallocate(partitionS, valueS, errorS)
 end

Fortran Results

Below are the results
332.00000000000000        155.96759681601489        9.7047371403391480E-004

Matlab version

a = 0;
b = 5.+85;
par1 = 100*(atan(b)-atan(a));
fa = 100/(1+a.^2);
fb = 100/(1+b.^2);

fprintf('METHOD SIMPSON\n');
for n = 1:1000000
    h=(b-a)/n;
    sum1=0;
    sum2=0;
    x1 = a;
    x2 = a;
    for j = 1:n-1
        k = j;
        if fix(k/2) ~= j/2 
            if j == 1
                x1 = x1+h;
            end
            if j > 1
                x1 = x1+2*h;
            end   
            f1 = 100/(1+x1.^2);
            sum1 = sum1 + f1;
        else
            x2 = x2+2*h;
            f2 = 100/(1+x2.^2);
            sum2 = sum2 + f2;
        end
    end
    par2 = (h/3)*(fa+4*sum1+2*sum2+fb);
    e = par1 - par2;
    if abs(e)<0.001
        break;
    end
end



y=n;
partitionS = zeros (n);
valueS= zeros (n);
errorS = zeros (n);


for n = 1:y
    h=(b-a)/n;
    sum1=0;
    sum2=0;
    x1=a;
    x2=a;
    for j = 1:n-1
        k = j;
        if fix(k/2) == j/2
            x2 = x2 + 2*h;
            f2 = 100/(1+x2.^2);
            sum2 = sum2 + f2;
        else
            if j == 1
                x1 = x1 + h;
            end
            if j > 1
                x1 = x1 + 2*h;
            end
            f1 = 100/(1+x2.^2);
            sum1 = sum1 + f1;
        end
    end
    partitionS(n) = n;
    valueS(n)= (h/3)*(fa+4*sum1+2*sum2+fb);
    errorS(n)=par1-valueS(n);
end

fprintf('Below are the results\n');
fprintf('%.25f\n',partitionS(n));
fprintf('%.25f\n',valueS(n));
fprintf('%.25f\n',errorS(n));

MATLAB Results

Below are the results
332.0000000000000000000000000
174.0415303853845900000000000
-18.0729630956556660000000000

Is Fortran's `k` [an integer](http://stackoverflow.com/q/24488518/)? — francescalus, Mar 25 '17 at 20:17
Are you sure it's an infinite loop? There's no `while` statement. — sco1, Mar 25 '17 at 21:17
If k/2 /= j/2.0 tests whether j is odd, can we also translate it to "if mod( j, 2 ) == 1" in Matlab? — roygvib, Mar 26 '17 at 03:27
I updated my answer before declaring I found it: You wrote `f1 = 100/(1+x2.^2);` instead of `f1 = 100/(1+x1.^2);`. Please read my answer to the end. (I hope that's it). — Rotem, Mar 26 '17 at 19:23
That was the bug for getting different results @Rotem . The problem have solved :) Many thanks to all for the help. — ktsigkounis, Mar 26 '17 at 20:19

Rotem · Accepted Answer · 2017-03-26T16:42:42.030

Like francescalus commented, it looks like problem is related to integer arithmetic in Fortran.

You may modify the first if statement in Matlab implementation as follows:

if fix(k/2) ~= j/2

In your second part, there is a typo error in the Matlab code.
You wrote x2 instead of x1.

Correct code:

f1 = 100/(1+x1.^2); %Instead of f1 = 100/(1+x2.^2);

Minor flaw:

if abs(e)<=0.001 %Instead of if abs(e)<0.001

I know very basic Fortran, so I executed both Matlab and Fortran code versions side by side.
I executed the code step by step using the debugger.
I used some arbitrary input values.

The problem is related to the first Fortran if statement: (k/2/=j/2.)
When k is an integer k/2 evaluates to floor(k/2), and j/2. evaluates to floating point (assume k is positive).
(I used fix Matlab function, in case k can also be negative).

Example:

integer j, k
j=3
k=3

print *, k/2
print *, j/2.
print *, k/2/=j/2.

Result:

           1
   1.500000
 T

In Matlab, the default type is double.

j=3;
k=3;

disp(k/2)
disp(j/2)
disp(k/2 ~= j/2)

Result:

As you can see, in Fortran condition evaluates to true, and in Matlab to false.

Complete Matlab code:

a = 0;
b = 5.+85;
par1 = 100*(atan(b)-atan(a));
fa = 100/(1+a.^2);
fb = 100/(1+b.^2);

fprintf('METHOD SIMPSON\n');
for n = 1:1000000
    h=(b-a)/n;
    sum1=0;
    sum2=0;
    x1 = a;
    x2 = a;
    for j = 1:n-1
        k = j;
        if fix(k/2) ~= j/2 
            if j == 1
                x1 = x1+h;
            end
            if j > 1
                x1 = x1+2*h;
            end   
            f1 = 100/(1+x1.^2);
            sum1 = sum1 + f1;
        else
            x2 = x2+2*h;
            f2 = 100/(1+x2.^2);
            sum2 = sum2 + f2;
        end
    end
    par2 = (h/3)*(fa+4*sum1+2*sum2+fb);
    e = par1 - par2;
    if abs(e)<=0.001
        break;
    end
end



y=n;
partitionS = zeros (n);
valueS= zeros (n);
errorS = zeros (n);


for n = 1:y
    h=(b-a)/n;
    sum1=0;
    sum2=0;
    x1=a;
    x2=a;
    for j = 1:n-1
        k = j;
        if fix(k/2) == j/2
            x2 = x2 + 2*h;
            f2 = 100/(1+x2.^2);
            sum2 = sum2 + f2;
        else
            if j == 1
                x1 = x1 + h;
            end
            if j > 1
                x1 = x1 + 2*h;
            end
            f1 = 100/(1+x1.^2);%f1 = 100/(1+x2.^2);
            sum1 = sum1 + f1;
        end
    end
    partitionS(n) = n;
    valueS(n)= (h/3)*(fa+4*sum1+2*sum2+fb);
    errorS(n)=par1-valueS(n);
end

fprintf('Below are the results\n');
fprintf('%.25f\n',partitionS(n));
fprintf('%.25f\n',valueS(n));
fprintf('%.25f\n',errorS(n));

Matlab output:

METHOD SIMPSON
Below are the results
332.0000000000000000000000000
155.9675968160148900000000000
0.0009704737140339148000000

Ι tried both ways (yours: `fix(k/2)` and @roygvib 's : `mod(j,2) == 1`) and the program doesn't stuck. Although I don't get the same results as Fortran program. — ktsigkounis, Mar 26 '17 at 10:18
Can you please post a minimal executable version of your code. The Fortran code must include declarations of all variables (e.g integer, real*8 ...). All variables must be initialized to relevant values before the loop. The loop should be executed in the main function. Also post Matlab version, with the same values of variables. For reference, post your Matlab and Fortran execution result (final value of `sum1` and `sum2` at least). — Rotem, Mar 26 '17 at 10:47

score 1 · Answer 2 · answered Mar 26 '17 at 11:31

I made a small fortran program based on your posts. Then put it through my f2matlab fortran source to matlab source converter (matlab file exchange). Here is the fortran:

program kt_f
implicit none
integer j,n,k,f1,f2
real x1,x2,h,sum1,sum2

n=100
k=50

do j=1,n-1
 k=j
 if(k/2/=j/2.) then
  if(j==1) x1=x1+h
  if(j>1) x1=x1+2*h
  f1=100/(1+x1**2)
  sum1=sum1+f1
 else
  x2=x2+2*h
  f2=100/(1+x2**2)
  sum2=sum2+f2
 endif
enddo

print *,'sum1=',sum1
print *,'sum2=',sum2

end program kt_f

When I compile and run this, the output is:

sum1=   5000.000    
sum2=   4900.000

Here is the matlab source produced. Note that in addition to the fix in the if statement, you need another fix in the line with the 100/ because this is an integer division as well. Here is the matlab code:

function kt_f(varargin)
 clear global; clear functions;
 global GlobInArgs nargs
 GlobInArgs={mfilename,varargin{:}}; nargs=nargin+1;
 persistent f1 f2 h_fv j k n sum1 sum2 x1 x2 ; 

 if isempty(f1), f1=0; end;
 if isempty(f2), f2=0; end;
 if isempty(h_fv), h_fv=0; end;
 if isempty(j), j=0; end;
 if isempty(k), k=0; end;
 if isempty(n), n=0; end;
 if isempty(sum1), sum1=0; end;
 if isempty(sum2), sum2=0; end;
 if isempty(x1), x1=0; end;
 if isempty(x2), x2=0; end;

 n = 100;
 k = 50;

 for j = 1: n - 1;
  k = fix(j);
  if(fix(k./2) ~= (j./2.));
   if(j == 1);
    x1 = x1 + h_fv;
   end;
   if(j > 1);
    x1 = x1 + 2.*h_fv;
   end;
   f1 = fix(100./(1+x1.^2));
   sum1 = sum1 + f1;
  else;
   x2 = x2 + 2.*h_fv;
   f2 = fix(100./(1+x2.^2));
   sum2 = sum2 + f2;
  end;
 end;

'sum1=',sum1
'sum2=',sum2
end %program kt_f

This gives the same output as the fortran. Please check and see whether this solves you issue.

Fortran to Matlab confused to if statement

2 Answers2