トップ 差分 一覧 ソース 検索 ヘルプ RSS ログイン

BugTrack-その他のメモ/4

R備忘録 /状態空間モデリング/donlp2/その他のメモ

FORTRAN - QUADPACKをCから使う?

  • 投稿者: みゅ
  • カテゴリ: FORTRAN
  • 優先度: 普通
  • 状態: 完了
  • 日時: 2009年09月23日 14時11分07秒

内容

  • QUADPACKをCから使う?
  • Cで記述した関数を積分する

サンプル

  • test02.c
#include "stdio.h"
#include "math.h"

int inc1_(int *);  // FORTRANで書かれた関数.これを呼ぶ.
void timestamp();   // QUADPACKの関数.「_」なしで作成している
int i;

int main(int  argc,  char *argv[])
{
  i = 1;
  i = inc1_(&i);
  printf("i : %d\n", i);
  timestamp();
}

// 積分される関数
// FORTRAN側でrealで定義した場合は、C側でfloatで受ける.
float c02_(float *x){
  double a = *x;
  return(cos(100.0e+00*sin(a)));
}
  • test02f.f90
    • グローバル変数定義のサンプル
module teigi
  real::pi
end module
    • Cから呼ばれる口/グローバル変数がつかえてることも確認
integer function inc1(i)
! Fortran declaring PI global
!  COMMON /pi/ pi ! Common Block and variable have the same name
!  real pi
  use teigi
  implicit none
  integer i
  pi = 3.141592653589793E+00
  inc1 = i + 1
  call timestamp()
  call test01
  return
end
    • QUADPACKについていたサンプルのtest01関数
subroutine test01

!*****************************************************************************80
!
!! TEST01 tests QAG.
!
!  Discussion:
!
!    QAG is an adaptive automatic integrator using a Gauss-Kronrod rule.
!
!    integrate cos(100*sin(x)) from 0 to pi.
!
!    The exact answer is pi * j0(100), or roughly 0.06278740.
!
!    KEY chooses the order of the integration rule, from 1 to 6.
!
  use teigi
  implicit none

  real, parameter :: a = 0.0E+00
  real abserr
  real b
  real, parameter :: epsabs = 0.0E+00
  real, parameter :: epsrel = 0.001E+00
  real, external :: f02
  real, external :: c02
  integer ier
  integer, parameter :: key = 6
  integer neval
!  real, parameter :: pi = 3.141592653589793E+00
!  COMMON /pi/ pi ! Common Block and variable have the same name
!  real pi
  real result
  real, parameter :: true = 0.06278740E+00

  b = pi

!  call qag ( f02, a, b, epsabs, epsrel, key, result, abserr, neval, ier )
  call qag ( c02, a, b, epsabs, epsrel, key, result, abserr, neval, ier )
!  ↑ここをc02としてCの関数を与えている

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Test QAG'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrand is COS(100*SIN(X))'
  write ( *, '(a,g14.6)' ) '  Integral left endpoint A =    ', a
  write ( *, '(a,g14.6)' ) '  Integral right endpoint B =   ', b
  write ( *, '(a,g14.6)' ) '  Exact integral is             ', true
  write ( *, '(a,g14.6)' ) '  Estimated integral is         ', result
  write ( *, '(a,g14.6)' ) '  Estimated integral error =    ', abserr
  write ( *, '(a,g14.6)' ) '  Exact integral error =        ', true - result
  write ( *, '(a,i8)' ) '  Number of function evaluations, NEVAL = ', neval
  write ( *, '(a,i8)' ) '  Error return code IER = ', ier

  return
end

function f02 ( x )

!*****************************************************************************80
!
!! F02 is the integrand function COS(100*SIN(X)).
!
  implicit none

  real f02
  real x

  f02 = cos ( 100.0E+00 * sin ( x ) )

  return
end
  • Makefile - これでうまくいくんじゃないかと思う
%.o: %.c
	gcc -c -o $@ $^

%.o: %.f90
	gcc -c -o $@ $^

test01 : test01.o test01f.o
	gcc $^ -o $@

test02 : test02.o test02f.o
	gfortran $^ -L. -lquadpack -o $@
    • FORTRANを含む実行ファイルを作成する場合は、FORTRAN関連のライブラリをリンクする必要があるので、gccではなく、gfortranをつかうと楽だそうで.
  • 実行
$ ./test02
September 23 2009   2:01:45.154 PM

TEST01
  Test QAG

  Integrand is COS(100*SIN(X))
  Integral left endpoint A =       0.00000
  Integral right endpoint B =      3.14159
  Exact integral is               0.627874E-01
  Estimated integral is           0.627873E-01
  Estimated integral error =      0.118180E-04
  Exact integral error =          0.141561E-06
  Number of function evaluations, NEVAL =      427
  Error return code IER =        0
i : 2
September 23 2009   2:01:45.164 PM
  • 元のFORTRANだけの場合
TEST01
  Test QAG

  Integrand is COS(100*SIN(X))
  Integral left endpoint A =       0.00000
  Integral right endpoint B =      3.14159
  Exact integral is               0.627874E-01
  Estimated integral is           0.627874E-01
  Estimated integral error =      0.117584E-04
  Exact integral error =           0.00000
  Number of function evaluations, NEVAL =      427
  Error return code IER =        0
    • 少しちがう
    • floatで計算してるからなのか?

コメント