Thread: 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.

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.

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.

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

Thank you very much john.c! It is very helpful.