Thread: Converting ancient FORTRAN code to C code

  1. #1
    Registered User Micko's Avatar
    Join Date
    Nov 2003
    Posts
    715

    Converting ancient FORTRAN code to C code

    Hello C experts,

    I have found Fortran code that I need for my research (it is about mathematical triangular matrix factorizations). I have never seen Fortran code in my life and would greatly appreciated if someone more familiar can rewrite this in C. I also haven't done any programming in C for yoears, but I still can read and understand C code.
    The Fortran code is in attachment.

    Thank you very much in advanced.
    Attached Images Attached Images Converting ancient FORTRAN code to C code-fortran-code-png 
    Gotta love the "please fix this for me, but I'm not going to tell you which functions we're allowed to use" posts.
    It's like teaching people to walk by first breaking their legs - muppet teachers! - Salem

  2. #2
    Registered User
    Join Date
    May 2012
    Location
    Arizona, USA
    Posts
    714
    I'm not fluent in FORTRAN, but here's what I do know.

    The numbers on the left side are line numbers. Those can be used to jump to a different place in the code.

    FORTRAN's "DO" statement is like C's "for" statement. The first number (such as 90 in the first "DO" statement) is the line number of the "terminal statement", or the last statement in the loop. "CONTINUE" (on line 90) just makes the loop continue to the next iteration, but it could be some other type of statement (such as the statement in line 60, which is also a terminal statement for the second "DO" statement).

    Here's what "I = 1,M" means in the "DO" statement:

    • "I" is the iterator variable and is initialized with the value 1.
    • After each iteration, "I" is incremented by the step, which is 1 by default. If "I" is greater than "M", the loop ends (jump to the statement after the terminal statement); otherwise, continue the next iteration.


    I don't know if this is the case in the program you posted, but FORTRAN treats all variables named "I" through "N" as integers and all others as "real" (floating point) numbers.

    I think "W(I,L)" is an array access (probably equivalent to "W[I][L]" in C).

    Besides those "quirks" of the FORTRAN language, the conversion should be fairly straightforward.
    Last edited by christop; 03-16-2020 at 02:06 PM.

  3. #3
    Registered User Micko's Avatar
    Join Date
    Nov 2003
    Posts
    715
    Thank you, I managed to figure out the most of it, but I'm not sure whether there are any nested loop here.
    For example:
    Code:
    m = n+k;
    for j=n:-1:1 
        D = 0;
        for i=1:m
            f(i)=W(j,i)*V(i);
            D = D+F(i)*W(j,i);
            U(j,j)=D;
            if(j==1) break;
    The code is in Scilab, which is Sw package I used most of the time for simple scripts.
    Gotta love the "please fix this for me, but I'm not going to tell you which functions we're allowed to use" posts.
    It's like teaching people to walk by first breaking their legs - muppet teachers! - Salem

  4. #4
    Registered User
    Join Date
    Dec 2017
    Posts
    928
    With Fortran's implicit typing, variables i,j,k,l,m,n are integers; others are real.
    Variables do not need to be declared (will equal 0 on first use ... I assume).
    The loops include their ending value.
    Arrays are 1-based.

    Here's a start. I left the loops working on 1-based arrays, so the arrays need to be declared 1 bigger than their nominal size.
    Code:
    void f( int n, int k )
    {
    // Array sizes (at least):
    // w[n+1][m+1], v[m+1], f[m+1], u[n+1][n+1];
    // these are probably also supposed to be inputs/outputs
     
        int m = n + k;
     
        for ( int j = n; j > 0; --j )
        {
            double d = 0.0;
     
            for ( int i = 1; i <= m; ++i )
            {
                f[i] = w[j][i] * v[i];
                d += f[i] * w[j][i];
            }
     
            u[j][j] = d;
     
            if ( j == 1 )
                 break; // end of loop
     
            if ( d == 0.0 )
            {
                for (int i = 1; i < j; ++i)
                    u[i][j] = 0;
                continue;
            }
     
            double lambda = 1.0 / d;
     
            for ( int i = 1; i < j; ++i )
            {
                double beta = 0.0;
                for ( int L = 1; L <= m; ++L )
                    beta += w[i][L] * f[L];
                beta *= lambda;
                u[i][j] = beta;
                for ( int L = 1; L <= m; ++L )
                    w[i][L] -= beta * w[j][L];
            }
        }
    }
    This is how I interpreted the picture of the Fortran code:
    Code:
    .   M = N + K
    
        DO 90 J = N, 1, -1
    
            D = 0.
            DO 60  I = 1, M
                F(I) = W(J, I) * V(I)
    60          D = D + F(I) * W(J, I)
            U(J, J) = D
    
            IF ( J = 1 ) GO TO 90
            IF ( D = 0. ) Go TO 80
    
            λ = 1. / D
    
            DO 75  I = 1, J - 1
              β = 0.
    
              DO 65  L = 1, M
    65            β = β + W(I, L) * F(L)
    
              β = λ * β
              U(I, J) = β
    
              DO 75  L = 1, M
    75            W(I, L) = W(I, L) - β * W(J, L)
    
            GO TO 90
    
    80      DO 85  I = 1, J-1
    85          U(I, J) = 0.
    
    90      CONTINUE
    EDIT: Here it is with 0-based arrays and a driver.
    Code:
    #include <stdio.h>
     
    #define N   5         // max n value
    #define K   2         // max k value
    #define M  (N + K)
     
    void func( int n, int k,
               double w[][M], double v[], double f[], double u[][N] )
    {
     
        int m = n + k;
     
        for ( int j = n; j-- > 0; )
        {
            double d = 0.0;
            for ( int i = 0; i < m; ++i )
            {
                f[i] = w[j][i] * v[i];
                d += f[i] * w[j][i];
            }
            u[j][j] = d;
     
            if ( j == 0 )
                break;
     
            if ( d == 0.0 )
                for (int i = 0; i < j - 1; ++i)
                    u[i][j] = 0;
            else
            {
                double lambda = 1.0 / d;
     
                for ( int i = 0; i < j - 1; ++i )
                {
                    double beta = 0.0;
                    for ( int x = 0; x < m; ++x )
                        beta += w[i][x] * f[x];
                    beta *= lambda;
                    u[i][j] = beta;
     
                    for ( int x = 0; x < m; ++x )
                        w[i][x] -= beta * w[j][x];
                }
            }
        }
    }
     
    void prn( double *a, int size )
    {
        for ( int m = 0; m < size; ++m )
            printf("%6.2f ", a[m]);
        putchar('\n');
    }
     
    void prn2d( double *a, int rows, int cols )
    {
        for ( int r = 0; r < rows; ++r )
        {
            prn(a, cols);
            a += cols;
        }
    }
     
    int main()
    {
        double w[N][M] =
        {
            { 1, 2, 3, 4, 5, 6, 7 },
            { 7, 6, 5, 4, 3, 2, 1 },
            { 1, 1, 1, 1, 1, 1, 1 },
            { 2, 2, 2, 2, 2, 2, 2 },
            { 3, 3, 3, 3, 3, 3, 3 }
        };
        double v[M] =
            { 1, 2, 3, 4, 5, 6, 7 };
        double f[M] =
            { 7, 6, 5, 4, 3, 2, 1 };
        double u[N][N] =
        {
            { 1, 2, 3, 4, 5 },
            { 1, 2, 3, 4, 5 },
            { 1, 2, 3, 4, 5 },
            { 1, 2, 3, 4, 5 },
            { 1, 2, 3, 4, 5 }
        };
     
        func(N, K, w, v, f, u);
     
        printf("w:\n");
        prn2d((double*)w, N, M);
        printf("u:\n");
        prn2d((double*)u, N, N);
        printf("v:\n");
        prn(v, M);
        printf("f:\n");
        prn(f, M);
     
        return 0;
    }
    Last edited by john.c; 03-16-2020 at 03:49 PM.
    If you want the truth to stand clear before you, never be for or against. - Sent-ts'an

  5. #5
    Registered User Micko's Avatar
    Join Date
    Nov 2003
    Posts
    715
    Thank you very much john.c! It is very helpful.
    Gotta love the "please fix this for me, but I'm not going to tell you which functions we're allowed to use" posts.
    It's like teaching people to walk by first breaking their legs - muppet teachers! - Salem

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Converting code
    By bananen in forum C Programming
    Replies: 4
    Last Post: 04-19-2017, 08:25 AM
  2. Converting, pumping Raw C code, Source code
    By secretcode in forum C Programming
    Replies: 5
    Last Post: 06-03-2016, 09:01 AM
  3. Converting C code to DLX assembly code
    By xyz3 in forum C Programming
    Replies: 2
    Last Post: 05-17-2010, 02:01 PM
  4. C Code for converting generic to binary code
    By vanilla in forum C Programming
    Replies: 5
    Last Post: 11-05-2009, 03:34 AM
  5. Replies: 1
    Last Post: 10-18-2005, 10:20 AM

Tags for this Thread