Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

Fortran Programming Language

Fortran to matlab infuriating problem


peace,

I have been translating some Numerical recipies programs to matlab
from C and Fortran.  Fortran is obviously closer to matlab than C.

I have successfully implemented tred2 and have been trying to
implement tqli...

I have been getting almost correct results from tqli (with treds
outputs as its input), but there is just one small error either in
syntax translation by me or I don't know what, and I cannot seem to
get it to function 100% correctly...  I have spent days at this and it
is driving me crazy now...

The fortran (77) which actually translates very nicely, i thought is

  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

my matlab is

function [d,z] = sceigen_test (d,e,z)

%input
%d is the diagonal vector of the tridiag matrix
%e is the sub and super diagonal of the tridiag matrix
%z is the transformation matrix from the tridiag program

% output,
%d is the eigenvalue vector
% z is the eigenvector matrix
% d is the eigenvale vector

% defining loop size based upon ip vector/matrix size
n = size(d,2);

%shifting e to the left, so that e(1) is not the zero anymore
for i = 2:n
    e(i-1) = e(i);
end
e(n)=0;
for l = 1:n
    iter=0;
    for m = l:n % not to n-1 for reason expained below
        if m == n
            break
        end
        dd = abs(d(m)) + abs(d(m+1));
        if ((abs(e(m))+ dd) == dd)
            break
        end
    end

   % m = n% this is implemented in the above loop, so the loop goes to
n instead of n-1

    if m ~= l
        if iter == 30
            error('bloody typical..too many iterations');
        end
        iter = iter+1;
        g = (d(l+1)-d(l))/(2.0*e(l));
        r = pythag(g, 1.0);%it works correctly
        g = d(m) - d(l) + (e(l) / (g + SIGNO(r, g)));%SIGNO = SIGN...
definitely
        s = 1.0;
        c = 1.0;
        p = 0.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.0
                d(i+1)= d(i+1)-p;
                e(m) = 0.0;
                break
            end
            s = f/r;
            c = g/r;
            g = d(i + 1) - p;
            r = (d(i) - g)*s + 2.0*c*b;
            p = s*r;
            d(i + 1) = g + p;
            g = c * r - b;
            %eigenvector part now...
            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
        if r ~= 0.0
            d(l) = d(l) - p;
            e(l) = g;
            e(m) = 0.0;
        end
    end
end

The input (example)

which is (if anyone will test for themselves)

[0.0861,0.0172,0.0466,-0.0186,-0.0384;0.0172,0.0635,-0.0146,-0.0154,0.0399; 0.0466,-0.0146,0.0555,-0.0166,-0.0563;-0.0186,-0.0154,-0.0166,0.0459,-0.000 7;-0.0384,0.0399,-0.0563,-0.0007,0.0704;]

which is a covariance matrix, so it is symmetric etc...

eigenvalues (ordered hi to low)

eig_sort =

  Columns 1 through 3

   0.171639537105982                   0
0
                   0                              0.100853395915414
0
                   0
0                              0.037197771768781
                   0
0                              0
                   0
0                              0

  Columns 4 through 5

                   0                               0
                   0                               0
                   0                               0
   0.011737362737036                    0
                   0                               0.000002142655106

I have eigenvectors also, but i want to focus on theeigenvalues for
now

Tha actual matlab generates eigenvalues (reverse ordered) are:

de =

  Columns 1 through 3

   0.000000000000000                   0
0
                   0                              0.011739410520354
0
                   0
0                              0.037197816074304
                   0
0                              0
                   0
0                              0

  Columns 4 through 5

                   0                   0
                   0                   0
                   0                   0
   0.100853446274617        0
                   0                   0.171639537313043

They are similar, but there always a constant error...  I noticed
also, that if I subtract eig_sort(5,5) from itself and add it to
eig_sort(4,4), the results become more accurate...

But this is probably just indicative of some other problem....

It would be great if someone could see what I was doing wrong and
point it out..

Peace,
Marwan

the eigenvales look untidy...

Here they are again...

mine are

0.171639537105982      0.100853395915414      0.037197771768781
0.011737362737036      0.000002142655106

matlab's are (in reverse order to my own)

 0.000000000000000       0.011739410520354        0
0.037197816074304
 0.100853446274617       0   0.171639537313043

Marwan wrote:
> peace,

> I have been translating some Numerical recipies programs to matlab
> from C and Fortran.  Fortran is obviously closer to matlab than C.

> I have successfully implemented tred2 and have been trying to
> implement tqli...

> I have been getting almost correct results from tqli (with treds
> outputs as its input), but there is just one small error either in
> syntax translation by me or I don't know what, and I cannot seem to
> get it to function 100% correctly...  I have spent days at this and it
> is driving me crazy now...

[a bunch of code snipped...]

> They are similar, but there always a constant error...  I noticed
> also, that if I subtract eig_sort(5,5) from itself and add it to
> eig_sort(4,4), the results become more accurate...

> But this is probably just indicative of some other problem....

> It would be great if someone could see what I was doing wrong and
> point it out..

Let me preface my reply by stating that I'm not trying to be a smart-ass, but....

What you are doing wrong is converting NR C and Fortran code into matlab code into the
first place.

Why not use the matlab functionality to do eigenvalue decomposition/analysis (or any other
form of linear algebra)? That's one of the things (if not the main one) matlab was built
for; and it's arguably more robust than NR routines.

Is there a reason you are converting NR codes into matlab?

cheers,

paulv

--
Paul van Delst             Ride lots.
CIMSS @ NOAA/NCEP/EMC               Eddy Merckx

On 15 May 2007 07:56:27 -0700, Marwan <marwanboust@gmail.com>
 wrote in <1179240987.251515.86@y80g2000hsf.googlegroups.com>:

> the eigenvales look untidy...
> Here they are again...
> mine are
> 0.171639537105982      0.100853395915414      0.037197771768781
> 0.011737362737036      0.000002142655106
> matlab's are (in reverse order to my own)
>  0.000000000000000       0.011739410520354        0
> 0.037197816074304
>  0.100853446274617       0   0.171639537313043

        What precision does Matlab work in?  Your fortran routine declares
everything REAL (despite your output of apparently DOUBLE PRECISION
eigenvalues).  What results do you get if you change your reals to
doubles?  (You'll probably need to change your input/output arrays to
double as well).

--
Ivan Reid, School of Engineering & Design, _____________  CMS Collaboration,
Brunel University.    Ivan.Reid@[brunel.ac.uk|cern.ch]    Room 40-1-B12, CERN
        KotPT -- "for stupidity above and beyond the call of duty".

On May 15, 10:49 am, Marwan <marwanboust@gmail.com> wrote:

> peace,

> I have been translating some Numerical recipies programs to matlab
> from C and Fortran.  Fortran is obviously closer to matlab than C.

The best newsgroup to find people to debug your Matlab code is
comp.soft-sys.matlab . Maybe it would be easier to translate Fortran
90 code using array operations to Matlab -- there is an F90 version of
Numerical Recipes.

Since Matlab's original purpose was to remove the need for users to
code their own linear algebra algorithms, I wonder why you are not
using the built-in functions. If you want Matlab code for some
algorithm for pedagogical reasons, the code may already exist on the
Matlab file exchange or some other place.

I am merely using matlab as a testing area... I actually need a fully
expanded code for implementation elsewhere...

anyway hopefully someone might be able to help.

Marwan wrote:
> I am merely using matlab as a testing area...

Hello,

For testing of what? The Fortran/C NR code? If so, compare the straight Fortran/C code
with the native Matlab results.

> I actually need a fully
> expanded code for implementation elsewhere...

Which code? The Fortran/C NR code, or its matlab translation?

And what do you mean by "expanded"? Is the final implementation going to be done using
matlab, or using the Fortran/C NR codes?

If the reason you are using matlab is to test the NR codes (which is a good idea), then
curiosity insists I ask again: why are you translating the NR Fortran/C source to matlab?
Why not just compare the NR Fortran/C output with the matlab output (using the built-in
matlab linear algebra routines) for a canned, typical dataset?

cheers,

paulv

--
Paul van Delst             Ride lots.
CIMSS @ NOAA/NCEP/EMC               Eddy Merckx

Because I will be coding for hardware and I will not have the
privilege of calling matlab functions...  so i want a code that is
fully expanded...  That is all the steps for the calculation are in
the code.

For this reason, I wanted code which i test on matlab, without any
functions.  I can then use this for testing as I write in verilog
etc...

peace...

On May 15, 12:26 pm, Marwan <marwanboust@gmail.com> wrote:

> Because I will be coding for hardware and I will not have the
> privilege of calling matlab functions...  so i want a code that is
> fully expanded...  That is all the steps for the calculation are in
> the code.

> For this reason, I wanted code which i test on matlab, without any
> functions.  I can then use this for testing as I write in verilog
> etc...

I don't follow the logic here, but...

Did you try Dr Reid's advice re: double precision?  Matlab computes
everything in double precision by default and as he notes your Fortran
routine (at least as posted) shows single precision for a floating
point variables.  Until I had both at the same precision for sure, I'd
not be looking for other differences.

On 15 May 2007 14:52:38 -0700, dpb <bozart@gmail.com>
 wrote in <1179265957.972665.11@o5g2000hsb.googlegroups.com>:

> On May 15, 12:26 pm, Marwan <marwanboust@gmail.com> wrote:
>> Because I will be coding for hardware and I will not have the
>> privilege of calling matlab functions...  so i want a code that is
>> fully expanded...  That is all the steps for the calculation are in
>> the code.
>> For this reason, I wanted code which i test on matlab, without any
>> functions.  I can then use this for testing as I write in verilog
>> etc...
> I don't follow the logic here, but...
> Did you try Dr Reid's advice re: double precision?  Matlab computes
> everything in double precision by default and as he notes your Fortran
> routine (at least as posted) shows single precision for a floating
> point variables.  Until I had both at the same precision for sure, I'd
> not be looking for other differences.

        Verilog is used (inter alia?) to program Field-Programmable Gate
Arrays.  The only time I've seen its input, it was in a C-like syntax.
So there's probably not much in the way of libraries[1], even if there were
a way to access a "subroutine block" in hardware.  Caveat, I've been
involved in projects using FPGAs, but I've not actually programmed one
myself.

http://en.wikipedia.org/wiki/FPGA

[1] Except for "code" that was always placed "inline".

--
Ivan Reid, School of Engineering & Design, _____________  CMS Collaboration,
Brunel University.    Ivan.Reid@[brunel.ac.uk|cern.ch]    Room 40-1-B12, CERN
GSX600F, RG250WD         "You Porsche. Me pass!"   DoD #484     JKLO#003, 005
WP7# 3000   LC Unit #2368 (tinlc)   UKMC#00009   BOTAFOT#16    UKRMMA#7 (Hon)
        KotPT -- "for stupidity above and beyond the call of duty".

Basically.

I will be coding in verilog.  There will not by a host of functions
for me to call, so functions like eigenvalues and eigenvectors will
need to be coded explicitly.

Yes Verilog is being used by me for FPGA programming.

Basically, my issue, is when I try to directly translate the TQLI.F OR
TQLI.C code to matlab, it just will not work for reason totally beyond
me...

TRED2.f worked perfectly...  which add to the irritation ironically...

here is hoping someone who has this done sees this post!

peace.

Marwan wrote:
> I will be coding in verilog.  There will not by a host of functions
> for me to call, so functions like eigenvalues and eigenvectors will
> need to be coded explicitly.

I have added comp.arch.fpga and comp.lang.verilog, where someone
may have worked on similar problems.

> Basically, my issue, is when I try to directly translate the TQLI.F OR
> TQLI.C code to matlab, it just will not work for reason totally beyond
> me...

Someone else can help figure out Matlab.  I will note that the hardware
implementations of algorithms tend to be somewhat different from the
usual software implementations, but you should understand the software
first.  Matlab seems to me a fine way to do that.

I recommend systolic arrays for hardware matrix processing.  The usual
reason for wanting a hardware implementation is speed, and systolic
arrays usually work pretty well for that.  Also, they scale with problem
size when needed.  About how big are your matrices and how fast do you
need the results?

-- glen

Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc