Converting A Fortran77 Code To Matlab Code for finding eigen values/vectors

Question

I converted a written code in fortran 77 to Matlab code. This function computes the eigenvalues and eigenvectors of a matrix using QL algorithm. for some reasons I can't use the eig function's result in matlab. The obtained eigenvalues from this method is not identical to those obtained by eig function, some of them are the same but some differs. I don't know where is the problem. Thank you for any help and suggestion. I can give the input arrays, if needed for running and watching the results.

here is the fortran code:

      SUBROUTINE tqli(d,e,n,np,z)
      INTEGER n,np
      REAL d(np),e(np),z(np,np)
CU    USES pythag
      INTEGER i,iter,k,l,m
      REAL b,c,dd,f,g,p,r,s,pythag
      do 11 i=2,n
        e(i-1)=e(i)
11    continue
      e(n)=0.
      do 15 l=1,n
        iter=0
1       do 12 m=l,n-1
          dd=abs(d(m))+abs(d(m+1))
          if (abs(e(m))+dd.eq.dd) goto 2
12      continue
        m=n
2       if(m.ne.l)then
          if(iter.eq.30)pause 'too many iterations in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.*e(l))
          r=pythag(g,1.)
          g=d(m)-d(l)+e(l)/(g+sign(r,g))
          s=1.
          c=1.
          p=0.
          do 14 i=m-1,l,-1
            f=s*e(i)
            b=c*e(i)
            r=pythag(f,g)
            e(i+1)=r
            if(r.eq.0.)then
              d(i+1)=d(i+1)-p
              e(m)=0.
              goto 1
            endif
            s=f/r
            c=g/r
            g=d(i+1)-p
            r=(d(i)-g)*s+2.*c*b
            p=s*r
            d(i+1)=g+p
            g=c*r-b
C     Omit lines from here ...
            do 13 k=1,n
              f=z(k,i+1)
              z(k,i+1)=s*z(k,i)+c*f
              z(k,i)=c*z(k,i)-s*f
13          continue
C     ... to here when finding only eigenvalues.
14        continue
          d(l)=d(l)-p
          e(l)=g
          e(m)=0.
          goto 1
        endif
15    continue
      return
      END

and the following is the matlab code:

function [d,z]=tqli(d,e,n,np,z)

for i=2:n
    e(i-1)=e(i);
end
e(n)=0.;
for l=1:n
    iter=0;
    for m=l:(n-1)
        dd=abs(d(m))+abs(d(m+1));
        if ((abs(e(m))+dd)==dd)
            break
        end
    end
    m=n;
    if (m~=l)
        if (iter==30)
            disp('too many iteration in tqli')
        end
        iter=iter+1;
        g=(d(l+1)-d(l))/(2.*e(l));
        r=pythag(g,1.);
        g=d(m)-d(l)+e(l)/(g+r*sign(g));
        s=1.;
        c=1.;
        p=0.;
        for i=(m-1):-1:l
            f=s*e(i);
            b=c*e(i);
            r=pythag(f,g);
            e(i+1)=r;
            if(r==0.)
                d(i+1)=d(i+1)-p;
                e(m)=0.;
                break
            end
            s=f/r;
            c=g/r;
            g=d(i+1)-p;
            r=(d(i)-g)*s+2.*c*b;
            p=s*r;
            d(i+1)=g+p;
            g=c*r-b;
            for k=1:n
                f=z(k,i+1);
                z(k,i+1)=s*z(k,i)+c*f;
                z(k,i)=c*z(k,i)-s*f;
            end
        end
        d(l)=d(l)-p;
        e(l)=g;
        e(m)=0.;
    end
end
end

Answer 1

In your Fortran code you declare a bunch of variables as REAL. Most compilers, by default, implement these as 32-bit floating-point numbers. The corresponding variables in your Matlab version are, by default, 64-bit floating-point numbers.

Without seeing either your inputs or outputs it's difficult to tell whether this is, partly or wholly, the cause of the differences in the outputs from the two versions. But, changing the precision of floating-point numbers is often the cause of the sort of problems you report, especially in tricky numerical methods such as eigenvalue computation.

While I'm writing, I also observe that you have translated the bad practice of comparing floating-point values for equality from Fortran to Matlab. It's bad practice in Fortran, and it's bad practice in Matlab.

Be careful in Matlab when you write expressions such as

2.*c

In Fortran that multiplies each element of the scalar-or-array variable c by the REAL value 2.0. In Matlab it multiplies each element of the scalar-or-vector c by the integer 2, .* is an operator (or combination of operators if you prefer) in Matlab, meaning 'perform multiplication element-by-element'. You have changed, subtly, the semantics of your program and possibly the precise value of some of the numbers you calculate.

This leads me to my final comment: what you call your Matlab version is no more than a transliteration of the Fortran code, a line-by-line copying from one language to the other. You can get away with it some of the time, but sooner or later this naive approach will trip you up. Perhaps it already has. You really should rewrite your Matlab into Matlab as if you intended to write good Matlab.

Answer 2

I saw a couple of things that may cause trouble with your matlab translation. One was your conversion of the fortran sign. You need to use abs(r) instead of just r.

The other more serious problem I saw was your flow structure resulting from your trying to refactor the goto's. When I converted this with f2matlab (and an unreleased tool that I wrote, remgoto) it came up with the following flow structure. I hope this helps you on your way. Note this is all untested!

function [d,e,n,np,z]=tqli(d,e,n,np,z);

remg([1:2])=true;

for i = 2 : n;
 e(i-1) = e(i);
end
e(n) = 0.;
for l = 1 : n;
 while (1);
  if(remg(2))
   if(remg(1))
    iter = 0;
   end;
   remg(1)=true;
   for m = l : n - 1;
    dd = abs(d(m)) + abs(d(m+1));
    if( abs(e(m))+dd==dd )
     remg(2)=false;
     break;
    end;
   end;
   if(~(remg(2)))
    continue;
   end;
   m = fix(n);
  end;
  remg(2)=true;
  if( m~=l )
   if( iter==30 )
    disp(['too many iterations in tqli',' -- Hit Return to continue']);
    pause ;
   end;
   iter = fix(iter + 1);
   g =(d(l+1)-d(l))./(2..*e(l));
   r = pythag(g,1.);
   g = d(m) - d(l) + e(l)./(g+(abs(r).*sign(g)));
   s = 1.;
   c = 1.;
   p = 0.;
   for i = m - 1 : -1: l ;
    f = s.*e(i);
    b = c.*e(i);
    r = pythag(f,g);
    e(i+1) = r;
    if( r==0. )
     d(i+1) = d(i+1) - p;
     e(m) = 0.;
     remg(1)=false;
     break;
    end;
    s = f./r;
    c = g./r;
    g = d(i+1) - p;
    r =(d(i)-g).*s + 2..*c.*b;
    p = s.*r;
    d(i+1) = g + p;
    g = c.*r - b;
    %        Omit lines from here ...
    for k = 1 : n;
     f = z(k,i+1);
     z(k,i+1) = s.*z(k,i) + c.*f;
     z(k,i) = c.*z(k,i) - s.*f;
    end; k = fix(n+1);
    %        ... to here when finding only eigenvalues.
   end;
   if(~(remg(1)))
    continue;
   end;
   d(l) = d(l) - p;
   e(l) = g;
   e(m) = 0.;
   remg(1)=false;
   continue;
  end;
  break;
 end;
end;
end %subroutine tqli