[SOLVED] Some questions in C/Fortran code

ejspeiro · 11-30-2012, 12:00 PM

Greetings guys,

I am exploring the creation of wrappers functions for Fortran from C code. I decided to do it by implementing a simple matrix-vector multiplication code, for which I have created two source files:

1. The first one is the routine in Fortran:

Code:

      subroutine product(nr, nc, aa, bb, rr)

        integer nr
        integer nc
        double precision aa(100,100)
        double precision bb(100)
        double precision rr(100)

        do ii = 1, nr
          rr(ii) = 0
          do jj = 1, nc
            rr(ii) = rr(ii) + aa(ii,jj)*bb(jj)
          enddo
        enddo

      end subroutine

In the latter snippet of code, nr and nc stand for number of rows and number of columns, respectively, aa is the matrix, bb is the vector, and rr is the solution vector. Please excuse my Fortran code being so basic, but I am not really a Fortran programmer.

2. The second source file is the actual C driver code which invokes the latter procedure:

Code:

#include <stdio.h>

//! Execution constants:

#define II 100
#define JJ 100

//! Main module.
// This module stands for the main application we are interesting in creating
// the wrapper for, i.e., Sym.8. In this case, the main module requires a
// matrix-vector multiplication. Both the matrix and the vector, along with the
// solution vector will be feed from this main module, which is written in
// C/C++ to the function implementing such operation, which is written in
// Fortran. This example should be general enough since arrays are being passed
// as arguments, from C/C++ to Fortran.

int main() {

  double aa[II][JJ]; // Matrix aa.
  double bb[II];     // Vector bb, so we are interested in computing: aa*bb.
  double rr[II];     // Solution vector so that: rr = aa*bb.
  int nr;  // Number of rows of both the matrix aa and the vector bb.
  int nc;  // Number of columns of the matrix aa.
  int ii;  // Iterator.
  int jj;  // Iterator.


  //! Initialize input values:
  nr = 5;
  nc = 5;
  // Create matrix and vector:
  for (ii = 0; ii < nr; ii++) {
    for (jj = 0; jj < nc; jj++){
      aa[ii][jj] = (double) ii*jj;
    }
    bb[ii] = (double) ii;
  }
  
  for (ii = 0; ii < nr; ii++) {
    for (jj = 0; jj < nc; jj++) {
      fprintf(stdout, "%7.2f", aa[ii][jj]);
    }
    fprintf(stdout, "\t%7.2f", bb[ii]);
    putchar('\n');
  }

  //! Call product function written in Fortran:
  product(nr, nc, aa, bb, rr);

  fprintf(stdout, "Solution vector:\n\n\t[");
  for (ii = 0; ii < nr; ii++){
    fprintf(stdout, "%7.2f", rr[ii]);
  }
  fprintf(stdout, " ]^T\n");
  return 0;
}

//! Procedure: matrix_vector_product
// This is the wrapper function right here. We notice the following VITAL
// details:
//
// 1. The use of the classifier "extern". This increases the visibility of this
// function so that it is available from other files comprising the project.
// Please see:
// http://stackoverflow.com/questions/1...variables-in-c
//
// 2. Notice that all of the argument are memory positions, AND MEMORY POSITIONS
// ONLY! This is important because once the compiler does its work and
// translates everything into object code, the only thing different programming
// languages have in common is accessing memory. Therefore, we must provide
// arguments that are compatible to those that will be used by the OTHER
// programming languages. Since every programming language access memory, this
// already ensures compatibility, EXCEPT for that we also need to make them
// compatible in terms of how much memory do these require. For that, we use
// the native data-types within the languages, since these are nothing but
// mnemonics for this. In this case, we use int and void.
// Notice that the use of memory position is important, not only for the
// definition of the wrapper function, but also for the actual invocation of
// the EXTERNAL procedure: we also pass memory positions.
//
// 2.1. Notice that even though we declare the arguments to be memory positions,
// we do NOT use the & operator when invoking the function in the body of them
// main module. I do not know why :(
// Please see:
// http://www.linuxquestions.org/questi...31#post4840331
//
// 3. Finally, notice the use of the trailing under score character. This is
// a convention when using external function to avoid the compiler of the
// "client programming language" to think this has to be treated as a recursive
// procedure.

extern product(int *nr, int *nc, void *aa, void *bb, void *rr){

  product_(&nr, &nc, aa, bb, rr);
}

My questions are basically only 2:

Are the comments explaining the interactions between the two languages (written in the C driver file) correct?
Why am I not using the & operator when invoking the Fortran procedure form the C driver? If I use it, it does NOT work, whereas if I don't, then it does.

Notice that the codes work and yield a correct results. I would like to educate myself in terms of what have I done.

Thanks in advanced.

linosaurusroot · 12-01-2012, 07:00 AM

Fortan is always (unless they've snuck new in features in the last 20 years) "call by value" rather than "call by reference".

suicidaleggroll · 12-01-2012, 08:49 AM

I've only glanced at the code, but I do see an issue. Multi-dimensional arrays in C are row-major, while in Fortran they're column-major. You need to swap your "nc" and "nr" indexing of aa in your Fortran code. I suggest changing the dimensions so they're not equal (for testing), and filling them with arbitrary values in C (maybe make it 3x5 and fill it with ((1,2,3),(4,5,6), etc) to make sure the correct value is ending up in the index you expect in Fortran.

Also, is this F77, F90, F95...? There are major differences between them. Based on the indentation on each line I'd say you're aiming for 77, however you should be aware that ANSI F77 should be written in all caps and that's not how the DO loops should be written in F77. They should be written as below:

Code:

        DO 10 II = 1, NR
          RR(II) = 0
          DD 20 JJ = 1, NC
            RR(II) = rr(II) + AA(II,JJ)*BB(JJ)
   20     CONTINUE
   10   CONTINUE

Alternatively, you could combine the two continue statements into one (I don't like this as it makes it more difficult to read IMO, but it's valid):

Code:

        DO 10 II = 1, NR
          RR(II) = 0
          DD 10 JJ = 1, NC
            RR(II) = rr(II) + AA(II,JJ)*BB(JJ)
   10   CONTINUE

I'd also HIGHLY suggest you turn off implicit declaration and explicitly declare your ii and jj arrays. you would turn this off by placing an "implicit none" at the start of your function. Implicit declaration needs to DIAF, and if you ever plan to do more Fortran in the future, you would do well to get in the habit of ALWAYS turning it off when you write a code.

H_TeXMeX_H · 12-01-2012, 09:41 AM

I have a question: If this can all be done in C, why bother with Fortran ?

ntubski · 12-01-2012, 01:57 PM

I don't know any Fortran but there are definitely some serious problems with the C code, you should compile with warnings enabled and pay attention to them. The main problem is that you call the function product() before it is declared which means that the C compiler has to guess the types of the parameters. The first two arguments (nr and nc) are ints, so even though the definition calls them pointers, the function gets int values.

Code:

// 1. The use of the classifier "extern". This increases the visibility of this
// function so that it is available from other files comprising the project.

Functions are "extern" by default, they can be called from other files as long as they are not declared "static".

Code:

extern product(int *nr, int *nc, void *aa, void *bb, void *rr){

The return type is missing here. If you are not returning a value, you should declare a return type of void.

makyo · 12-02-2012, 10:15 PM

Hi.

This may help: http://gcc.gnu.org/onlinedocs/gcc-4....-and-Functions

Best wishes ... cheers, makyo

ejspeiro · 12-03-2012, 11:04 AM

Quote:

Originally Posted by suicidaleggroll

I've only glanced at the code, but I do see an issue. Multi-dimensional arrays in C are row-major, while in Fortran they're column-major. You need to swap your "nc" and "nr" indexing of aa in your Fortran code. I suggest changing the dimensions so they're not equal (for testing), and filling them with arbitrary values in C (maybe make it 3x5 and fill it with ((1,2,3),(4,5,6), etc) to make sure the correct value is ending up in the index you expect in Fortran.

suicidaleggroll: Thanks for this warning!
suicidaleggroll and H_TeXMeX_H: This is the actual reason of why I am so interested in developing this code, for a problem that can be done in C: To better understand these interaction aspects. My ultimate purpose is to understand the wrapping mechanisms that are available to invoke FORTRAN code from ScaLAPACK, in a C application.

Quote:

Originally Posted by suicidaleggroll

Also, is this F77, F90, F95...? There are major differences between them. Based on the indentation on each line I'd say you're aiming for 77, however you should be aware that ANSI F77 should be written in all caps and that's not how the DO loops should be written in F77. They should be written as below:

Code:

        DO 10 II = 1, NR
          RR(II) = 0
          DD 20 JJ = 1, NC
            RR(II) = rr(II) + AA(II,JJ)*BB(JJ)
   20     CONTINUE
   10   CONTINUE

Alternatively, you could combine the two continue statements into one (I don't like this as it makes it more difficult to read IMO, but it's valid):

Code:

        DO 10 II = 1, NR
          RR(II) = 0
          DD 10 JJ = 1, NC
            RR(II) = rr(II) + AA(II,JJ)*BB(JJ)
   10   CONTINUE

suicidaleggroll: Thanks for this as well, and you are right: I recycled some old code from my undergrad days and it is intended to be FORTRAN 77.

Quote:

Originally Posted by suicidaleggroll

I'd also HIGHLY suggest you turn off implicit declaration and explicitly declare your ii and jj arrays. you would turn this off by placing an "implicit none" at the start of your function. Implicit declaration needs to DIAF, and if you ever plan to do more Fortran in the future, you would do well to get in the habit of ALWAYS turning it off when you write a code.

suicidaleggroll: Thanks for this as well!

ejspeiro · 12-03-2012, 11:08 AM

Quote:

Originally Posted by ntubski

I don't know any Fortran but there are definitely some serious problems with the C code, you should compile with warnings enabled and pay attention to them. The main problem is that you call the function product() before it is declared which means that the C compiler has to guess the types of the parameters. The first two arguments (nr and nc) are ints, so even though the definition calls them pointers, the function gets int values.

Code:

// 1. The use of the classifier "extern". This increases the visibility of this
// function so that it is available from other files comprising the project.

Functions are "extern" by default, they can be called from other files as long as they are not declared "static".

Code:

extern product(int *nr, int *nc, void *aa, void *bb, void *rr){

The return type is missing here. If you are not returning a value, you should declare a return type of void.

ntubski: You are completely right! I am aware that the old "this was just a quick-n-dirty code to get some quick results" SHOULD never be an excuse, and I appreciate your comments! Specifically,

Quote:

Originally Posted by ntubski

(...) you should compile with warnings enabled and pay attention to them.

Which specific flags would you recommend? For example, I often use -Wall. Any suggestions?

Thanks.

H_TeXMeX_H · 12-03-2012, 11:33 AM

'-Wall' is usually enough to give you warnings. '-pedantic' may also be useful if you really want to be compliant to the standards ... sometimes needlessly so, for example this does not allow C++ '//' quoting.

ejspeiro · 12-03-2012, 11:34 AM

Quote:

Originally Posted by H_TeXMeX_H

'-Wall' is usually enough to give you warnings. '-pedantic' may also be useful if you really want to be compliant to the standards ... sometimes needlessly so, for example this does not allow C++ '//' quoting.

Thanks! \m/

ejspeiro · 12-03-2012, 01:42 PM

Dear all,

thank you very much for your replies. This is the improved FORTRAN 77 code:

Code:

      SUBROUTINE PRODUCT(NR, NC, AA, BB, RR)
      
        IMPLICIT NONE
        
        INTEGER II
        INTEGER JJ
        INTEGER NR
        INTEGER NC
        DOUBLE PRECISION AA(100,100)
        DOUBLE PRECISION BB(100)
        DOUBLE PRECISION RR(100)

        DO 10 JJ = 1, NC
          RR(JJ) = 0
          DO 20 II = 1, NR
            RR(JJ) = RR(JJ) + AA(JJ,II)*BB(II)
   20     CONTINUE
   10   CONTINUE

      END SUBROUTINE

And the C code:

Code:

#include <stdio.h>

//! Execution constants:

#define II 100
#define JJ 100

//! Prototypes:

void extern product(int *nr, int *nc, void *aa, void *bb, void *rr);
void product_(int *nr, int *nc, void *aa, void *bb, void *rr);

//! Main module:

// This module stands for the main application we are interesting in creating
// the wrapper for, i.e., Sym.8. In this case, the main module requires a
// matrix-vector multiplication. Both the matrix and the vector, along with the
// solution vector will be feed from this main module, which is written in
// C/C++ to the function implementing such operation, which is written in
// Fortran. This example should be general enough since arrays are being passed
// as arguments, from C/C++ to Fortran.

int main() {

  double aa[II][JJ]; // Matrix aa.
  double bb[II];     // Vector bb, so we are interested in computing: aa*bb.
  double rr[II];     // Solution vector so that: rr = aa*bb.
  int nr;  // Number of rows of both the matrix aa and the vector bb.
  int nc;  // Number of columns of the matrix aa.
  int ii;  // Iterator.
  int jj;  // Iterator.


  //! Initialize input values:
  nr = 5;
  nc = 2;
  // Create matrix and vector:
  for (ii = 0; ii < nr; ii++) {
    for (jj = 0; jj < nc; jj++){
      aa[ii][jj] = (double) ii*jj;
    }
  }
  for (ii = 0; ii < nc; ii ++) {
    bb[ii] = (double) ii;
  }
  // Print them:
  for (ii = 0; ii < nr; ii++) {
    for (jj = 0; jj < nc; jj++) {
      fprintf(stdout, "%7.2f", aa[ii][jj]);
    }
    putchar('\n');
  }
  fprintf(stdout, "times\n\n\t[");
  for (ii = 0; ii < nc; ii++){
    fprintf(stdout, "%7.2f", bb[ii]);
  }
  fprintf(stdout, " ]^T\n");
  
  //! Compute result in C:
  for (ii = 0; ii < nr; ii++) {
    rr[ii] = 0.0;
    for (jj = 0; jj < nc; jj++) {
      rr[ii] = rr[ii] + aa[ii][jj]*bb[jj];
    }
  }
  fprintf(stdout, "Solution vector (computed in C):\n\n\t[");
  for (ii = 0; ii < nr; ii++){
    fprintf(stdout, "%7.2f", rr[ii]);
  }
  fprintf(stdout, " ]^T\n");

  //! Call product function written in C (for testing purposes):
  product(&nr, &nc, aa, bb, rr);

  fprintf(stdout, "Solution vector (computed in Fortran 77):\n\n\t[");
  for (ii = 0; ii < nr; ii++){
    fprintf(stdout, "%7.2f", rr[ii]);
  }
  fprintf(stdout, " ]^T\n");
  return 0;
}

//! Procedure: product

// This is the wrapper function right here. We notice the following VITAL
// details:
//
// 1. The use of the classifier "extern". This increases the visibility of this
// function so that it is available from other files comprising the project.
// Please see:
// http://stackoverflow.com/questions/1433204/what-are-extern-variables-in-c
//
// 2. Notice that all of the argument are memory positions, AND MEMORY POSITIONS
// ONLY! This is important because once the compiler does its work and
// translates everything into object code, the only thing different programming
// languages have in common is accessing memory. Therefore, we must provide
// arguments that are compatible to those that will be used by the OTHER
// programming languages. Since every programming language access memory, this
// already ensures compatibility, EXCEPT for that we also need to make them
// compatible in terms of how much memory do these require. For that, we use
// the native data-types within the languages, since these are nothing but
// mnemonics for this. In this case, we use int and void.
// Notice that the use of memory position is important, not only for the
// definition of the wrapper function, but also for the actual invocation of
// the EXTERNAL procedure: we also pass memory positions.
//
// 2.1. Notice that even though we declare the arguments to be memory positions,
// we do NOT use the & operator when invoking the function in the body of them
// main module. I do not know why :(
// Please see:
// http://www.linuxquestions.org/questions/showthread.php?p=4840331#post4840331
//
// 3. Finally, notice the use of the trailing under score character. This is
// a convention when using external function to avoid the compiler of the
// "client programming language" to think this has to be treated as a recursive
// procedure.

void extern product(int *nr, int *nc, void *aa, void *bb, void *rr){

  product_(nr, nc, aa, bb, rr);
}

They both yield the same output:

Code:

[ejspeiro@dulcinea wrapper]$ make run
./main
   0.00   0.00
   0.00   1.00
   0.00   2.00
   0.00   3.00
   0.00   4.00
times

        [   0.00   1.00 ]^T
Solution vector (computed in C):

        [   0.00   1.00   2.00   3.00   4.00 ]^T
Solution vector (computed in Fortran 77):

        [   0.00   1.00   2.00   3.00   4.00 ]^T
[ejspeiro@dulcinea wrapper]$

Thanks! I will call this thread as solved, nonetheless, if you have further suggestions, I will totally welcome them.

ntubski · 12-03-2012, 07:18 PM

Quote:

Originally Posted by ejspeiro

Which specific flags would you recommend? For example, I often use -Wall. Any suggestions?

I usually use -Wall and -Wextra.

Quote:

if you have further suggestions, I will totally welcome them.

Don't use void* when you actually know what the type should be, don't use the do-nothing wrapping function:

Code:

#include <stdio.h>

//! Execution constants:

#define II 100
#define JJ 100

//! Prototypes:

// defined in Fortran code
extern void product_(int *nr, int *nc, double aa[II][JJ], double bb[II], double rr[II]);

//! Main module:

...

int main() {
    ...
    //! Call product_ function written in Fortran 77
    product_(&nr, &nc, aa, bb, rr);
    ...
}

What's the "^T" for?

ejspeiro · 12-03-2012, 07:59 PM

Quote:

Originally Posted by ntubski

don't use the do-nothing wrapping function:

Ummm interesting? I am curious about what the reason could be?

Quote:

Originally Posted by ntubski

What's the "^T" for?

It means transpose. I find it useful to keep track of when are vector column-size or row-wise.

Thanks.

ntubski · 12-03-2012, 09:24 PM

Quote:

Originally Posted by ejspeiro

Quote:

don't use the do-nothing wrapping function:

Ummm interesting? I am curious about what the reason could be?

Uh, why do you want a function that does nothing? There is no reason for having it, so remove it.