!!!FORTRAN② - QUADPACKをCから使う？
*投稿者： みゅ
*カテゴリ： FORTRAN
*優先度： 普通
*状態： 完了
*日時： 2009年09月23日 14時11分07秒
//{{bugstate}}
!!内容
*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関連のライブラリをリンクする必要があるので、ｇｃｃではなく、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で計算してるからなのか？
!!コメント
//{{comment}}