0 Replies - 3447 Views - Last Post: 29 July 2012 - 07:52 AM Rate Topic: -----

#1 maryam yousefi  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 1
  • Joined: 29-July 12

question about converting a fortran77 code to Matlab code

Posted 29 July 2012 - 07:52 AM

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 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


Is This A Good Question/Topic? 0
  • +

Page 1 of 1