- 投稿者： みゅ
- カテゴリ： FORTRAN
- 優先度： 普通
- 状態： 完了
- 日時： 2009年09月23日 14時11分07秒
- QUADPACKをCから使う？
- 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)));
}
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
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
Ｒ備忘録 /状態空間モデリング/donlp2/その他のメモ