C から Fortran を呼び出してみる

概要

• C言語で書いたコードから Fortran を呼び出す。
• C言語で作成した行列を Fortran で扱ってみる。
  • (都合上、使用したのは Fortran 90 です。)

情報

Fortran って？？

• 参照したページ
  • 日本NAGのFortran入門
    • NAG (Numerical Algorithms Group) が Fortran90 を世界初出した。と沿革に書いてあった。
    • ので、Fortan についてはココ見ると良さそう。

• 自分なりの要約
  • Fortran は数値計算に向いている。
  • 数値計算ライブラリが大量に揃っている。
  • 移植性に優れている。
  • 上位互換を(ほぼ)全て保っている。
    • Fortran 77 の機能は切り捨てられず、上位バージョンでも動く。

環境

• Linux、CentOS 6.4。
• Fortran は「gfortran」で、C は「gcc」でコンパイル。
• Fortran のソースコードは「Fortran 90」で記述 (.f90)。

Hello World

• そもそも Fortran 初めてだったので、Hello World から。
  • プログラムは「program <名前> ～ end program」で囲む
  • print は、始めにフォーマットを記述。その後に表示対象を書く。
  • 大文字・小文字は区別されない。

ソースコード

• hello.f90

program hello
  print *, "hello"
end program hello

コンパイル & 実行

• gfortran でコンパイル後、実行

$ gfortran hello.f90
$ ./a.out
 hello

C との連携 ( 簡単な例 )

• oracle の、CとFortranの連携について書かれたページを参考に簡単な例を作成してみる。
• int と double の値を渡して、それを処理してみる。

ソースコード

• hello.c

#include <stdio.h>

extern void fsim_( int *i, double *r );


int main() {

    int i = 100;
    double r = 23.0;

    fsim_( &i, &r );

    printf( "%d, %lf \n", i, r );

    return 0;
}

• hello.f90

subroutine fsim( i, r )
    integer i
    double precision r

    i = i + 15
    r = r * 25.8

    return
end subroutine fsim

コンパイル & 実行

• C のコードは gcc、Fortran のコードは gfortran + リンク。

$ gcc -c hello.c
$ gfortran hello.f90 hello.o
$ ./a.out

ポイント

• Fortran では、呼び出す処理を subroutine として定義する。
  • C の int は "integer"、double は "double precision"。
• C では、Fortran で定義した subroutine を extern して使用。
  • Fortran で定義した名前にアンダースコア「_」を付けて extern する。

C との連携 ( 1次元配列 )

• 1次元配列を読み込んで Fortran 上で処理してみる。
• Fortran には、配列と、配列のサイズを渡すことで処理してもらう。

ソースコード

• main.c

#include <stdio.h>

extern void add_array_( int[], int * );


int main() {

    int i;

    int length = 10;
    int array[length];
    for ( i = 0; i < length; i++ ) {
        array[i] = i + 1;
    }

    printf( "Before : \n" );
    for ( i = 0; i < length; i++ ) {
        printf( "%d ", array[i] );
    }
    printf( "\n" );

    add_array_( array, &length );

    printf( "After : \n" );
    for ( i = 0; i < length; i++ ) {
        printf( "%d ", array[i] );
    }
    printf( "\n" );

    return 0;
}

• func.f90

subroutine add_array( array, length )
    integer length
    integer array( length )

    integer i
    do i = 1, length
        array(i) = array(i) + 10
    end do

    return
end subroutine add_array

コンパイル & 実行

• 始めの例と同様

$ gcc -c main.c
$ gfortran func.f90 main.o
$ ./a.out
Before :
1 2 3 4 5 6 7 8 9 10
After :
11 12 13 14 15 16 17 18 19 20

C との連携 ( 2次元配列 ( = 行列 ) )

• 簡単な例で見ていたページでは C 上2次元配列を使用していた。
  • こちらの場合、引数渡しする際に『 a[][10] 』のように、column 数を指定する必要がある。
  • C のアドレス参照の問題による。詳しくは「C 2次元配列 引数」でググると良さそう。

• 不定サイズの行列を渡すことを考えて、次のようにしてみる。
  • C では、1次元配列として扱う。その状態で Fortran に渡す。
  • Fortran 上で、2次元配列に整形。
  • その上で、Fortran 上で処理を行う。
  • 1次元配列に戻して、C 側に戻す。

• 下記の例では、Fortran 上 2次元配列化して、列ごとに処理してみた。
  • Fortran は列優先なので、その仕組みに乗ってみた形。

ソースコード

• main.c

#include <stdio.h>
#include <stdlib.h> // Random

#define RANDOM_SEED (2)
#define THRESHOLD   (50)

extern void multiple_matrix_( double *, int *, int * );

int main() {

    int i, j;

    int row = 4, col = 5; // 4行5列
    double matrix[row * col];

    srand( RANDOM_SEED );
    for ( i = 0; i < row; i++ ) {
        for ( j = 0; j < col; j++ ) {
            matrix[j+i*col] = ( rand() % 100 < THRESHOLD ) ? 1.0 : 0.0;
        }
    }

    printf( "Before : \n" );
    for ( i = 0; i < row; i++ ) {
        for ( j = 0; j < col; j++ ) {
            printf( "%.4lf ", matrix[j+i*col] );
        }
        printf( "\n" );
    }
    printf( "\n" );

    multiple_matrix_( matrix, &row, &col );

    printf( "After : \n" );
    for ( i = 0; i < row; i++ ) {
        for ( j = 0; j < col; j++ ) {
            printf( "%.4lf ", matrix[j+i*col] );
        }
        printf( "\n" );
    }
    printf( "\n" );

    return 0;
}

• func.f90

subroutine multiple_matrix( org_matrix, row, col )
    integer row, col
    double precision org_matrix( row * col )
    double precision matrix( col, row )
    integer i, j, count

    ! Array => Matrix
    do i = 1, col
        do j = 1, row
            matrix( i, j ) = org_matrix( i + ( j - 1 ) * col )
        end do
    end do

    ! Processing
    count = 1
    do i = 1, col
        do j = 1, row
            matrix( i, j ) = matrix( i, j ) * count
            count          = count + 1
        end do
    end do

    ! Matrix => Array
    do i = 1, col
        do j = 1, row
            org_matrix( i + ( j - 1 ) * col ) = matrix( i, j )
        end do
    end do

    return
end subroutine multiple_matrix

コンパイル & 実行

• 1次元の時と同様

$ gcc -c main.c
$ gfortran func.f90 main.o
$ ./a.out
Before :
0.0000 1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 1.0000 1.0000
1.0000 0.0000 1.0000 0.0000 1.0000
1.0000 0.0000 0.0000 0.0000 1.0000

After :
0.0000 5.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 14.0000 18.0000
3.0000 0.0000 11.0000 0.0000 19.0000
4.0000 0.0000 0.0000 0.0000 20.0000

備考

• Fortan では、行列サイズが大きいとセグメンテーション違反を起こすらしい。
  • マシンの使用スタック領域サイズを増やしたりして、解消するようだ。
  • 「Fortan 巨大 配列」でググると、そのような問題が出てくる。

• ので、上記の方法は本当はあまり良くないのかも。。