以下例子均来自网络，只是稍作了编辑，方便今后查阅。
子目录
（一） Fortran调用C语言
（二） C语言调用Fortran
（三） C++ 调用Fortran
（四） Fortran 调用 C++
需要说明的是，（一）和（二）对GCC编译器的版本要求并不高；而（三）和（四）对GCC编译器的要求比较高，需要GCC在4.7及以上才能编译通过，这是由于自Fortran 2003一代语言起，增加了一个名为“iso_c_binding”的模块，可以很方便地在Fortran和C++之间传递参数，简化了两者的混合编程。
————————————————————我是分割线———————————————————————————————
（一） Fortran调用C语言
源程序来自网络，增加了一个二维数组类型的调用示例。
Fortran文件（main.f90）的内容如下：
program main
implicit none
interface
subroutine sub_c(n1,n2,n3)
integer :: n1
real :: n2
real() :: n3
end subroutine real() function func_d(var2d)
real() :: var2d(,)
end function real() function func_c(n3)
real() :: n3
end function end interface integer :: n1
real :: n2
real() :: n3,n4
real() :: var2d(,) n1=
n2=5.0
call sub_c(n1,n2,n3)
n4=func_c(n3)
write(*,*) "n1=",n1
write(*,*) "n2=",n2
write(*,*) "n3=",n3
write(*,*) "n4=",n4 var2d=
var2d(,)=dble(n1)
var2d(,)=dble(n2)
write(*,*) "var2d(1,:): ",var2d(,:)
write(*,*) "var2d(2,:): ",var2d(,:)
n2=func_d(var2d)
write(*,*) "var2d(1,:): ",var2d(,:)
write(*,*) "var2d(2,:): ",var2d(,:) end program main
C文件（sub.c）的内容如下：
#include <math.h>
#include <stdio.h> void sub_c_(int *,float *,double *);
double func_c_(double *);
double func_d_(double [][]); void sub_c_(int *n1,float *n2,double *n3)
{
*n3=pow(*n2,*n1);
} double func_c_(double *n3)
{
double n4;
n4=sqrt(*n3);
return n4;
} double func_d_(double var2d[][])
{
double NN;
NN=77.0;
printf("%5.2f %5.2f %5.2f \n",var2d[][],var2d[][],var2d[][]);
printf("%5.2f %5.2f %5.2f \n",var2d[][],var2d[][],var2d[][]);
var2d[][]=NN;
printf("%5.2f %5.2f %5.2f \n",var2d[][],var2d[][],var2d[][]);
printf("%5.2f %5.2f %5.2f \n",var2d[][],var2d[][],var2d[][]);
return NN;
}
Makefile 文件的内容如下：
all:
gcc -c sub.c -o sub.o
gfortran -c main.f90 -o main.o
gfortran -o main.exe main.o sub.o clean:
rm *.o *.exe
编译并运行的结果如下：
[She@she-centos7 TestABC]$ make
gcc -c sub.c -o sub.o
gfortran -c main.f90 -o main.o
gfortran -o main.exe main.o sub.o
[She@she-centos7 TestABC]$ ./main.exe
n1=
n2= 5.00000000
n3= 125.00000000000000
n4= 11.180339887498949
var2d(,:): 3.0000000000000000 5.0000000000000000 0.0000000000000000
var2d(,:): 0.0000000000000000 0.0000000000000000 0.0000000000000000
3.00 5.00 0.00
0.00 0.00 0.00
3.00 77.00 0.00
0.00 0.00 0.00
var2d(,:): 3.0000000000000000 77.000000000000000 0.0000000000000000
var2d(,:): 0.0000000000000000 0.0000000000000000 0.0000000000000000
注意：二维数组在Fortran和C中存储的方式不同，Fortran中是“列优先”，而C中是“行优先”。
在参数传递时，容易导致两者的引用错误。使用时要务必注意这个区别！！！
（二） C语言调用Fortran
示例1. 用结构体在C和Fortran之间传递参数
这种方法不需要对函数体进行预先声明，直接引用即可。由于结构体本身的最小存储单元的限制，在结构体内部的数据类型最好先作优化，以节省内存开销，原理可参见《C语言结构体占用空间内存大小解析》。
c2fortran.c
#include <string.h>
struct test{
int n;
float m;
};
/* c2fortran.c */
int main(int argc,char *argv[])
{
struct test t;
t.n=;
t.m=0.003;
int i;
float e = 2.71828;
char hello[];
int length = sizeof(hello);
strcpy(hello,"Hello Fortran from C");
for(i=strlen(hello); i<length; i++)
hello[i] = ' ';
showhie_(hello,&length,&e,&t);
return();
}
showhie.f
C showhie.f
C
SUBROUTINE SHOWHIE(HELLO,LENGTH,E,T)
CHARACTER*(*) HELLO
INTEGER LENGTH
REAL E
TYPE TEST
INTEGER N
REAL M
END TYPE
TYPE (TEST) :: T
C
WRITE(*,) HELLO(:LENGTH),LENGTH,E,T%N,T%M
FORMAT(3X,A,2X,I3,4X,F7.,2X,I3,2X,F7.)
RETURN
END SUBROUTINE SHOWHIE
Makefile_test 的内容：
fortran2c : c2fortran.o showhie.o
gfortran c2fortran.o showhie.o -o c2fortran showhie.o : showhie.f
gfortran -c showhie.f -o showhie.o c2fortran.o : c2fortran.c
gcc -c c2fortran.c -o c2fortran.o clean :
rm *.o
编译并运行的结果如下：
[She@she-centos7 TestABC]$ make -f Makefile_test
gcc -c c2fortran.c -o c2fortran.o
gfortran -c showhie.f -o showhie.o
gfortran c2fortran.o showhie.o -o c2fortran
[She@she-centos7 TestABC]$ ./c2fortran
Hello Fortran from C 2.71828 0.00300
main.c
#include <stdio.h> void sub_fortran_(int *,float *,double *);
double function_fortran_(double *); int main()
{
int num_int;
float num_float;
double num_double;
double num;
num_int=;
num_float=5.0;
sub_fortran_(&num_int,&num_float,&num_double);
num=function_fortran_(&num_double);
printf("num_int=%d\nnum_float=%f\nnum_double=%f\nnum=%f",num_int,num_float,num_double,num);
return ;
}
sub.f90
subroutine Sub_Fortran(NumInt,NumFloat,NumDouble)
implicit none
integer :: NumInt
real :: NumFloat
real() :: NumDouble
NumDouble=NumFloat**NumInt
end subroutine real() function Function_Fortran(NumDouble)
implicit none
real() :: NumDouble
! Function_Fortran=sqrt(NumDouble) ! 用gfortran编译报错，pgf90也同样无法编译，不知何故...
Function_Fortran=NumDouble**0.5
end function
Makefile_CcallFortran
all:
gcc –o main.o –c main.c
gfortran –o sub.o –c sub.f90
gcc –o main2.exe main.o sub.o clean:
rm *.o main2.exe
编译及运行结果：
num=15625.000000[She@she-centos7 TestABC]$ make -f Makefile_CcallFortran
gcc -o main.o -c main.c
gfortran -o sub.o -c sub.f90
gcc -o main2.exe main.o sub.o
sub.o：在函数‘function_fortran_’中：
sub.f90:(.text+0x75)：对‘sqrt’未定义的引用
collect2: 错误：ld 返回
make: *** [all] 错误 [She@she-centos7 TestABC]$ make -f Makefile_CcallFortran
gcc -o main.o -c main.c
gfortran -o sub.o -c sub.f90
gcc -o main2.exe main.o sub.o
[She@she-centos7 TestABC]$ ./main2.exe
num_int=
num_float=5.000000
num_double=125.000000
num=15625.000000
（三） C++ 调用Fortran
本例要求：gcc 和 gfortran 的编译器版本为4.7及以上。
fortran_matrix_multiply.f90
subroutine matrix_multiply(A, rows_A, cols_A, B, rows_B, cols_B, C, rows_C, cols_C) bind(c)
use iso_c_binding
integer(c_int) :: rows_A, cols_A, rows_B, cols_B, rows_C, cols_C
real(c_double) :: A(rows_A, cols_A), B(rows_B, cols_B), C(rows_C, cols_C) C = matmul(A, B)
end subroutine matrix_multiply
cpp_main_1.cpp
#include <iostream>
#include <vector>
#include <random>
#include <ctime> using namespace std; //Fortran subroutine definition "translated" to C++
extern "C" {
void matrix_multiply(double *A, int *rows_A, int *cols_A, double *B, int *rows_B, int *cols_B, double *C, int *rows_C, int *cols_C);
} //Fill a vector with random numbers in the range [lower, upper]
void rnd_fill(vector<double> &V, double lower, double upper) { //use system clock to create a random seed
size_t seed (clock()); //use the default random engine and an uniform distribution
default_random_engine eng(seed);
uniform_real_distribution<double> distr(lower, upper); for( auto &elem : V){
elem = distr(eng);
}
} //Print matrix V(rows, cols) storage in column-major format
void print_matrix(vector<double> const &V, int rows, int cols) { for(int i = ; i < rows; ++i){
for(int j = ; j < cols; ++j){
std::cout << V[j * rows + i] << " ";
}
std::cout << std::endl;
}
std::cout << std::endl;
} int main() {
size_t N = ;
vector<double> A(N * N), B(N * N), C(N * N); //Fill A and B with random numbers in the range [0,1]:
rnd_fill(A, 0.0, 1.0);
rnd_fill(B, 0.0, 1.0); //Multiply matrices A and B, save the result in C
int sz = N;
matrix_multiply(&A[], &sz, &sz, &B[], &sz, &sz, &C[], &sz, &sz); //print A, B, C on the standard output
print_matrix(A, sz, sz);
print_matrix(B, sz, sz);
print_matrix(C, sz, sz); return ;
}
Makefile
all:
gfortran -c fortran_matrix_multiply.f90
g++ -c -std=c++ cpp_main_1.cpp
gfortran fortran_matrix_multiply.o cpp_main_1.o -lstdc++ -o mix_example.out clean:
rm *.o *.out
编译及运行结果如下：
[She@she-centos7 eg2]$ make
gfortran -c fortran_matrix_multiply.f90
g++ -c -std=c++ cpp_main_1.cpp
gfortran fortran_matrix_multiply.o cpp_main_1.o -lstdc++ -o mix_example.out
[She@she-centos7 eg2]$ ./mix_example.out
0.131538 0.678865 0.0345721
0.45865 0.934693 0.5297
0.218959 0.519416 0.00769819 0.131538 0.678865 0.0345721
0.45865 0.934693 0.5297
0.218959 0.519416 0.00769819 0.336233 0.741784 0.364408
0.60501 1.46015 0.515041
0.268717 0.638137 0.282764
本机的 gcc 和 gfortran 的编译器版本为4.8：
[She@she-centos7 eg2]$ gfortran -v
使用内建 specs。
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/usr/libexec/gcc/x86_64-redhat-linux/4.8./lto-wrapper
目标：x86_64-redhat-linux
配置为：../configure --prefix=/usr --mandir=/usr/share/man --infodir=/usr/share/info --with-bugurl=http://bugzilla.redhat.com/bugzilla --enable-bootstrap --enable-shared --enable-threads=posix --enable-checking=release --with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions --enable-gnu-unique-object --enable-linker-build-id --with-linker-hash-style=gnu --enable-languages=c,c++,objc,obj-c++,java,fortran,ada,go,lto --enable-plugin --enable-initfini-array --disable-libgcj --with-isl=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/isl-install --with-cloog=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/cloog-install --enable-gnu-indirect-function --with-tune=generic --with-arch_32=x86-64 --build=x86_64-redhat-linux
线程模型：posix
gcc 版本 4.8. (Red Hat 4.8.-) (GCC)
[She@she-centos7 eg2]$ gcc -v
使用内建 specs。
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/libexec/gcc/x86_64-redhat-linux/4.8./lto-wrapper
目标：x86_64-redhat-linux
配置为：../configure --prefix=/usr --mandir=/usr/share/man --infodir=/usr/share/info --with-bugurl=http://bugzilla.redhat.com/bugzilla --enable-bootstrap --enable-shared --enable-threads=posix --enable-checking=release --with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions --enable-gnu-unique-object --enable-linker-build-id --with-linker-hash-style=gnu --enable-languages=c,c++,objc,obj-c++,java,fortran,ada,go,lto --enable-plugin --enable-initfini-array --disable-libgcj --with-isl=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/isl-install --with-cloog=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/cloog-install --enable-gnu-indirect-function --with-tune=generic --with-arch_32=x86-64 --build=x86_64-redhat-linux
线程模型：posix
gcc 版本 4.8. (Red Hat 4.8.-) (GCC)
需要注意的是，这个程序对 gcc 和 gfortran 的版本比较敏感。它在CentOS 6.5 x64 & gcc version 4.4.7 20120313 (Red Hat 4.4.7-17) (GCC) 上测试时，编译报错，它并不支持 g++ 的 "-std=c++11" 选项，它自带的 C++ 编译选项为 “-std=c++98”，无法方便地实现 C++ 和 Fortran 之间的调用。
（四） Fortran 调用 C++
示例1. 来自《FORTRAN90 Program Calls C++ Function》
主程序 kronrod_prb.f90
program main !*****************************************************************************
!
!! MAIN is the main program for KRONROD_PRB.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! September
!
! Author:
!
! John Burkardt
!
use, intrinsic :: iso_c_binding! Fortran 2003 语言新加的模块。 implicit none
!
! TIMESTAMP is provided by the C++ library, and so the following
! INTERFACE block must set up the interface.
!
interface
subroutine timestamp ( ) bind ( c )
use iso_c_binding
end subroutine timestamp
end interface call timestamp ( ) write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'KRONROD_PRB:'
write ( *, '(a)' ) ' FORTRAN90 version'
write ( *, '(a)' ) ' FORTRAN90 test program calls C++ functions.' call test01 ( )
call test02 ( )
!
! Terminate.
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'KRONROD_PRB:'
write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' '
call timestamp ( ) stop
end
subroutine test01 ( ) !*****************************************************************************
!
!! TEST01 tests the code for the odd case N = .
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
use, intrinsic :: iso_c_binding implicit none
!
! KRONROD is provided by the C++ library, and so the following
! INTERFACE block must be set up to describe how data is to
! be passed.
!
interface
subroutine kronrod ( n, eps, x, w1, w2 ) bind ( c )
use iso_c_binding
integer ( c_int ), VALUE :: n
real ( c_double ), VALUE :: eps
real ( c_double ) :: x(*)
real ( c_double ) :: w1(*)
real ( c_double ) :: w2(*)
end subroutine kronrod
end interface integer ( c_int ), parameter :: n = real ( c_double ) eps
integer ( c_int ) i
integer ( c_int ) i2
real ( c_double ) s
real ( c_double ) w1(n+)
real ( c_double ) w2(n+)
real ( c_double ) :: wg(n) = (/ &
0.555555555555555555556D+, &
0.888888888888888888889D+, &
0.555555555555555555556D+ /)
real ( c_double ) :: wk(*n+) = (/ &
0.104656226026467265194D+, &
0.268488089868333440729D+, &
0.401397414775962222905D+, &
0.450916538658474142345D+, &
0.401397414775962222905D+, &
0.268488089868333440729D+, &
0.104656226026467265194D+ /)
real ( c_double ) x(n+)
real ( c_double ) :: xg(n) = (/ &
-0.77459666924148337704D+, &
0.0D+, &
0.77459666924148337704D+ /)
real ( c_double ) :: xk(*n+) = (/ &
-0.96049126870802028342D+, &
-0.77459666924148337704D+, &
-0.43424374934680255800D+, &
0.0D+, &
0.43424374934680255800D+, &
0.77459666924148337704D+, &
0.96049126870802028342D+ /) write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST01'
write ( *, '(a)' ) ' Request KRONROD to compute the Gauss rule'
write ( *, '(a)' ) ' of order 3, and the Kronrod extension of'
write ( *, '(a)' ) ' order 3+4=7.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Compare to exact data.' eps = 0.000001D+ call kronrod ( n, eps, x, w1, w2 ) write ( *, '(a)' ) ' '
write ( *, '(a,i2)' ) ' KRONROD returns 3 vectors of length ', n +
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' I X WK WG'
write ( *, '(a)' ) ' '
do i = , n +
write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) i, x(i), w1(i), w2(i)
end do write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Gauss Abscissas'
write ( *, '(a)' ) ' Exact Computed'
write ( *, '(a)' ) ' '
do i = , n
if ( * i <= n + ) then
i2 = * i
s = -1.0D+
else
i2 = * ( n + ) - * i
s = +1.0D+
end if
write ( *, '(2x,i4,2x,g14.6,2x,g14.6)' ) i, xg(i), s * x(i2)
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Gauss Weights'
write ( *, '(a)' ) ' Exact Computed'
write ( *, '(a)' ) ' '
do i = , n
if ( * i <= n + ) then
i2 = * i
else
i2 = * ( n + ) - * i
end if
write ( *, '(2x,i4,2x,g14.6,2x,g14.6)' ) i, wg(i), w2(i2)
end do write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Gauss Kronrod Abscissas'
write ( *, '(a)' ) ' Exact Computed'
write ( *, '(a)' ) ' '
do i = , * n +
if ( i <= n + ) then
i2 = i
s = -1.0D+
else
i2 = * ( n + ) - i
s = +1.0D+
end if
write ( *, '(2x,i4,2x,g14.6,2x,g14.6)' ) i, xk(i), s * x(i2)
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Gauss Kronrod Weights'
write ( *, '(a)' ) ' Exact Computed'
write ( *, '(a)' ) ' '
do i = , * n +
if ( i <= n + ) then
i2 = i
else
i2 = * ( n + ) - i
end if
write ( *, '(2x,i4,2x,g14.6,2x,g14.6)' ) i, wk(i), w1(i2)
end do return
end
subroutine test02 ( ) !*****************************************************************************
!
!! TEST02 tests the code for the even case N = .
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! November
!
! Author:
!
! John Burkardt
!
use, intrinsic :: iso_c_binding implicit none
!
! KRONROD is provided by the C++ library, and so the following
! INTERFACE block must be set up to describe how data is to
! be passed.
!
interface
subroutine kronrod ( n, eps, x, w1, w2 ) bind ( c )
use iso_c_binding
integer ( c_int ), VALUE :: n
real ( c_double ), VALUE :: eps
real ( c_double ) :: x(*)
real ( c_double ) :: w1(*)
real ( c_double ) :: w2(*)
end subroutine kronrod
end interface integer ( c_int ), parameter :: n = real ( c_double ) eps
integer ( c_int ) i
integer ( c_int ) i2
real ( c_double ) s
real ( c_double ) w1(n+)
real ( c_double ) w2(n+)
real ( c_double ) x(n+) write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST02'
write ( *, '(a)' ) ' Request KRONROD to compute the Gauss rule'
write ( *, '(a)' ) ' of order 4, and the Kronrod extension of'
write ( *, '(a)' ) ' order 4+5=9.' eps = 0.000001D+ call kronrod ( n, eps, x, w1, w2 ) write ( *, '(a)' ) ' '
write ( *, '(a,i2)' ) ' KRONROD returns 3 vectors of length ', n +
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' I X WK WG'
write ( *, '(a)' ) ' '
do i = , n +
write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) i, x(i), w1(i), w2(i)
end do return
end
C++函数 kronrod.cpp
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <ctime>
#include "kronrod.hpp" using namespace std; //****************************************************************************80 void abwe1 ( int n, int m, double eps, double coef2, bool even, double b[],
double *x, double *w ) //****************************************************************************80
//
// Purpose:
//
// ABWE1 calculates a Kronrod abscissa and weight.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 03 August 2010
//
// Author:
//
// Original FORTRAN77 version by Robert Piessens, Maria Branders.
// C++ version by John Burkardt.
//
// Reference:
//
// Robert Piessens, Maria Branders,
// A Note on the Optimal Addition of Abscissas to Quadrature Formulas
// of Gauss and Lobatto,
// Mathematics of Computation,
// Volume 28, Number 125, January 1974, pages 135-139.
//
// Parameters:
//
// Input, int N, the order of the Gauss rule.
//
// Input, int M, the value of ( N + 1 ) / 2.
//
// Input, double EPS, the requested absolute accuracy of the
// abscissas.
//
// Input, double COEF2, a value needed to compute weights.
//
// Input, bool EVEN, is TRUE if N is even.
//
// Input, double B[M+1], the Chebyshev coefficients.
//
// Input/output, double *X; on input, an estimate for
// the abscissa, and on output, the computed abscissa.
//
// Output, double *W, the weight.
//
{
double ai;
double b0;
double b1;
double b2;
double d0;
double d1;
double d2;
double delta;
double dif;
double f;
double fd;
int i;
int iter;
int k;
int ka;
double yy; if ( *x == 0.0 )
{
ka = ;
}
else
{
ka = ;
}
//
// Iterative process for the computation of a Kronrod abscissa.
//
for ( iter = ; iter <= ; iter++ )
{
b1 = 0.0;
b2 = b[m];
yy = 4.0 * (*x) * (*x) - 2.0;
d1 = 0.0; if ( even )
{
ai = m + m + ;
d2 = ai * b[m];
dif = 2.0;
}
else
{
ai = m + ;
d2 = 0.0;
dif = 1.0;
} for ( k = ; k <= m; k++ )
{
ai = ai - dif;
i = m - k + ;
b0 = b1;
b1 = b2;
d0 = d1;
d1 = d2;
b2 = yy * b1 - b0 + b[i-];
if ( !even )
{
i = i + ;
}
d2 = yy * d1 - d0 + ai * b[i-];
} if ( even )
{
f = ( *x ) * ( b2 - b1 );
fd = d2 + d1;
}
else
{
f = 0.5 * ( b2 - b0 );
fd = 4.0 * ( *x ) * d2;
}
//
// Newton correction.
//
delta = f / fd;
*x = *x - delta; if ( ka == )
{
break;
} if ( r8_abs ( delta ) <= eps )
{
ka = ;
}
}
//
// Catch non-convergence.
//
if ( ka != )
{
cout << "\n";
cout << "ABWE1 - Fatal error!\n";
cout << " Iteration limit reached.\n";
cout << " Last DELTA was " << delta << "\n";
exit ( );
}
//
// Computation of the weight.
//
d0 = 1.0;
d1 = *x;
ai = 0.0;
for ( k = ; k <= n; k++ )
{
ai = ai + 1.0;
d2 = ( ( ai + ai + 1.0 ) * ( *x ) * d1 - ai * d0 ) / ( ai + 1.0 );
d0 = d1;
d1 = d2;
} *w = coef2 / ( fd * d2 ); return;
}
//****************************************************************************80 void abwe2 ( int n, int m, double eps, double coef2, bool even, double b[],
double *x, double *w1, double *w2 ) //****************************************************************************80
//
// Purpose:
//
// ABWE2 calculates a Gaussian abscissa and two weights.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 03 August 2010
//
// Author:
//
// Original FORTRAN77 version by Robert Piessens, Maria Branders.
// C++ version by John Burkardt.
//
// Reference:
//
// Robert Piessens, Maria Branders,
// A Note on the Optimal Addition of Abscissas to Quadrature Formulas
// of Gauss and Lobatto,
// Mathematics of Computation,
// Volume 28, Number 125, January 1974, pages 135-139.
//
// Parameters:
//
// Input, int N, the order of the Gauss rule.
//
// Input, int M, the value of ( N + 1 ) / 2.
//
// Input, double EPS, the requested absolute accuracy of the
// abscissas.
//
// Input, double COEF2, a value needed to compute weights.
//
// Input, bool EVEN, is TRUE if N is even.
//
// Input, double B[M+1], the Chebyshev coefficients.
//
// Input/output, double *X; on input, an estimate for
// the abscissa, and on output, the computed abscissa.
//
// Output, double *W1, the Gauss-Kronrod weight.
//
// Output, double *W2, the Gauss weight.
//
{
double ai;
double an;
double delta;
int i;
int iter;
int k;
int ka;
double p0;
double p1;
double p2;
double pd0;
double pd1;
double pd2;
double yy; if ( *x == 0.0 )
{
ka = ;
}
else
{
ka = ;
}
//
// Iterative process for the computation of a Gaussian abscissa.
//
for ( iter = ; iter <= ; iter++ )
{
p0 = 1.0;
p1 = *x;
pd0 = 0.0;
pd1 = 1.0;
ai = 0.0;
for ( k = ; k <= n; k++ )
{
ai = ai + 1.0;
p2 = ( ( ai + ai + 1.0 ) * (*x) * p1 - ai * p0 ) / ( ai + 1.0 );
pd2 = ( ( ai + ai + 1.0 ) * ( p1 + (*x) * pd1 ) - ai * pd0 )
/ ( ai + 1.0 );
p0 = p1;
p1 = p2;
pd0 = pd1;
pd1 = pd2;
}
//
// Newton correction.
//
delta = p2 / pd2;
*x = *x - delta; if ( ka == )
{
break;
} if ( r8_abs ( delta ) <= eps )
{
ka = ;
}
}
//
// Catch non-convergence.
//
if ( ka != )
{
cout << "\n";
cout << "ABWE2 - Fatal error!\n";
cout << " Iteration limit reached.\n";
cout << " Last DELTA was " << delta << "\n";
exit ( );
}
//
// Computation of the weight.
//
an = n; *w2 = 2.0 / ( an * pd2 * p0 ); p1 = 0.0;
p2 = b[m];
yy = 4.0 * (*x) * (*x) - 2.0;
for ( k = ; k <= m; k++ )
{
i = m - k + ;
p0 = p1;
p1 = p2;
p2 = yy * p1 - p0 + b[i-];
} if ( even )
{
*w1 = *w2 + coef2 / ( pd2 * (*x) * ( p2 - p1 ) );
}
else
{
*w1 = *w2 + 2.0 * coef2 / ( pd2 * ( p2 - p0 ) );
} return;
}
//****************************************************************************80 void kronrod ( int n, double eps, double x[], double w1[], double w2[] ) //****************************************************************************80
//
// Purpose:
//
// KRONROD adds N+1 points to an N-point Gaussian rule.
//
// Discussion:
//
// This subroutine calculates the abscissas and weights of the 2N+1
// point Gauss Kronrod quadrature formula which is obtained from the
// N point Gauss quadrature formula by the optimal addition of N+1 points.
//
// The optimally added points are called Kronrod abscissas. The
// abscissas and weights for both the Gauss and Gauss Kronrod rules
// are calculated for integration over the interval [-1,+1].
//
// Since the quadrature formula is symmetric with respect to the origin,
// only the nonnegative abscissas are calculated.
//
// Note that the code published in Mathematics of Computation
// omitted the definition of the variable which is here called COEF2.
//
// Storage:
//
// Given N, let M = ( N + 1 ) / 2.
//
// The Gauss-Kronrod rule will include 2*N+1 points. However, by symmetry,
// only N + 1 of them need to be listed.
//
// The arrays X, W1 and W2 contain the nonnegative abscissas in decreasing
// order, and the weights of each abscissa in the Gauss-Kronrod and
// Gauss rules respectively. This means that about half the entries
// in W2 are zero.
//
// For instance, if N = 3, the output is:
//
// I X W1 W2
//
// 1 0.960491 0.104656 0.000000
// 2 0.774597 0.268488 0.555556
// 3 0.434244 0.401397 0.000000
// 4 0.000000 0.450917 0.888889
//
// and if N = 4, (notice that 0 is now a Kronrod abscissa)
// the output is
//
// I X W1 W2
//
// 1 0.976560 0.062977 0.000000
// 2 0.861136 0.170054 0.347855
// 3 0.640286 0.266798 0.000000
// 4 0.339981 0.326949 0.652145
// 5 0.000000 0.346443 0.000000
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 03 August 2010
//
// Author:
//
// Original FORTRAN77 version by Robert Piessens, Maria Branders.
// C++ version by John Burkardt.
//
// Reference:
//
// Robert Piessens, Maria Branders,
// A Note on the Optimal Addition of Abscissas to Quadrature Formulas
// of Gauss and Lobatto,
// Mathematics of Computation,
// Volume 28, Number 125, January 1974, pages 135-139.
//
// Parameters:
//
// Input, int N, the order of the Gauss rule.
//
// Input, double EPS, the requested absolute accuracy of the
// abscissas.
//
// Output, double X[N+1], the abscissas.
//
// Output, double W1[N+1], the weights for
// the Gauss-Kronrod rule.
//
// Output, double W2[N+1], the weights for
// the Gauss rule.
//
{
double ak;
double an;
double *b;
double bb;
double c;
double coef;
double coef2;
double d;
bool even;
int i;
int k;
int l;
int ll;
int m;
double s;
double *tau;
double x1;
double xx;
double y; b = new double[((n+)/)+];
tau = new double[(n+)/]; m = ( n + ) / ;
even = ( * m == n ); d = 2.0;
an = 0.0;
for ( k = ; k <= n; k++ )
{
an = an + 1.0;
d = d * an / ( an + 0.5 );
}
//
// Calculation of the Chebyshev coefficients of the orthogonal polynomial.
//
tau[] = ( an + 2.0 ) / ( an + an + 3.0 );
b[m-] = tau[] - 1.0;
ak = an; for ( l = ; l < m; l++ )
{
ak = ak + 2.0;
tau[l] = ( ( ak - 1.0 ) * ak
- an * ( an + 1.0 ) ) * ( ak + 2.0 ) * tau[l-]
/ ( ak * ( ( ak + 3.0 ) * ( ak + 2.0 )
- an * ( an + 1.0 ) ) );
b[m-l-] = tau[l]; for ( ll = ; ll <= l; ll++ )
{
b[m-l-] = b[m-l-] + tau[ll-] * b[m-l+ll-];
}
} b[m] = 1.0;
//
// Calculation of approximate values for the abscissas.
//
bb = sin ( 1.570796 / ( an + an + 1.0 ) );
x1 = sqrt ( 1.0 - bb * bb );
s = 2.0 * bb * x1;
c = sqrt ( 1.0 - s * s );
coef = 1.0 - ( 1.0 - 1.0 / an ) / ( 8.0 * an * an );
xx = coef * x1;
//
// Coefficient needed for weights.
//
// COEF2 = 2^(2*n+1) * n// * n// / (2n+1)//
//
coef2 = 2.0 / ( double ) ( * n + );
for ( i = ; i <= n; i++ )
{
coef2 = coef2 * 4.0 * ( double ) ( i ) / ( double ) ( n + i );
}
//
// Calculation of the K-th abscissa (a Kronrod abscissa) and the
// corresponding weight.
//
for ( k = ; k <= n; k = k + )
{
abwe1 ( n, m, eps, coef2, even, b, &xx, w1+k- );
w2[k-] = 0.0; x[k-] = xx;
y = x1;
x1 = y * c - bb * s;
bb = y * s + bb * c; if ( k == n )
{
xx = 0.0;
}
else
{
xx = coef * x1;
}
//
// Calculation of the K+1 abscissa (a Gaussian abscissa) and the
// corresponding weights.
//
abwe2 ( n, m, eps, coef2, even, b, &xx, w1+k, w2+k ); x[k] = xx;
y = x1;
x1 = y * c - bb * s;
bb = y * s + bb * c;
xx = coef * x1;
}
//
// If N is even, we have one more Kronrod abscissa to compute,
// namely the origin.
//
if ( even )
{
xx = 0.0;
abwe1 ( n, m, eps, coef2, even, b, &xx, w1+n );
w2[n] = 0.0;
x[n] = xx;
} delete [] b;
delete [] tau; return;
} //****************************************************************************80 double r8_abs ( double x ) //****************************************************************************80
//
// Purpose:
//
// R8_ABS returns the absolute value of an R8.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 14 November 2006
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, double X, the quantity whose absolute value is desired.
//
// Output, double R8_ABS, the absolute value of X.
//
{
double value; if ( 0.0 <= x )
{
value = x;
}
else
{
value = - x;
}
return value;
}
//****************************************************************************80 void timestamp ( ) //****************************************************************************80
//
// Purpose:
//
// TIMESTAMP prints the current YMDHMS date as a time stamp.
//
// Example:
//
// 31 May 2001 09:45:54 AM
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 July 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// None
//
{
# define TIME_SIZE static char time_buffer[TIME_SIZE];
const struct std::tm *tm_ptr;
size_t len;
std::time_t now; now = std::time ( NULL );
tm_ptr = std::localtime ( &now ); len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ); std::cout << time_buffer << "\n"; return;
# undef TIME_SIZE
}
C++头文件依赖 kronrod.hpp
extern "C"
{
void abwe1 ( int n, int m, double eps, double coef2, bool even, double b[],
double *x, double *w );
void abwe2 ( int n, int m, double eps, double coef2, bool even, double b[],
double *x, double *w1, double *w2 );
void kronrod ( int n, double eps, double x[], double w1[], double w2[] );
double r8_abs ( double x );
void timestamp ( );
}
两个版本的Makefile全文
注意：两个版本的Makefile，都可以正常运行并获得结果。区别在于版本一的Makefile 无法对这个.cpp文件进行正常调试，比如无法识别 pgdbg 中的 print 命令等。
Makefile 版本一（cpp文件无法调试）
all:
g++ -g -c -std=c++11 kronrod.cpp -o kronrod.o
pgf90 -g -c kronrod_prb.f90 -o kronrod_prb.o
pgf90 -g -lstdc++ -o main.exe kronrod.o kronrod_prb.o clean:
rm *.o *.exe
编译及运行过程如下：
[She@she-centos7 F90_call_CPP_eg1]$ make clean && make
rm *.o *.exe
g++ -g -c -std=c++ kronrod.cpp -o kronrod.o
pgf90 -g -c kronrod_prb.f90 -o kronrod_prb.o
pgf90 -g -lstdc++ -o main.exe kronrod.o kronrod_prb.o
[She@she-centos7 F90_call_CPP_eg1]$ pgdbg ./main.exe
pgdbg 调试这个 .cpp 文件，断点设在第 101 行：
pgdbg> Loaded: /home/She/Programs/TestABC/F90_call_CPP_eg1/./main.exe
MAIN_
pgdbg> () breakpoint set at: _IO_marker::abwe1 line: "kronrod.cpp"@ address: 0x40150F pgdbg> Breakpoint at 0x40150F, function _IO_marker::abwe1, file kronrod.cpp, line
#: ai = m + m + ; pgdbg> print ai
ERROR: Cannot evaluate [0x709a9d8] DWARF_CALL_FRAME_CFA: pgdbg> print m
ERROR: Cannot evaluate [0x7088e60] DWARF_CALL_FRAME_CFA: pgdbg> q
Makefile 版本二（cpp文件可以正常调试）
all:
pgc++ -g -c -std=c++11 kronrod.cpp -o kronrod.o
pgf90 -g -c kronrod_prb.f90 -o kronrod_prb.o
pgf90 -g -lstdc++ -o main.exe kronrod.o kronrod_prb.o clean:
rm *.o *.exe
编译及调试运行过程如下：
[She@she-centos7 F90_call_CPP_eg1]$ make clean && make
rm *.o *.exe
pgc++ -g -c -std=c++ kronrod.cpp -o kronrod.o
"kronrod.cpp", line : warning: variable "len" was set but never used
size_t len;
^ pgf90 -g -c kronrod_prb.f90 -o kronrod_prb.o
pgf90 -g -lstdc++ -o main.exe kronrod.o kronrod_prb.o
[She@she-centos7 F90_call_CPP_eg1]$ pgdbg ./main.exe
pgdbg 调试这个 .cpp 文件，断点同样设在第 101 行：
PGI Debugger 17.4-0 x86-64 (Workstation, 16 Process)
PGI Compilers and Tools
Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved. pgdbg> Loaded: /home/She/Programs/TestABC/F90_call_CPP_eg1/./main.exe
MAIN_
pgdbg> (1) breakpoint set at: abwe1 line: "kronrod.cpp"@101 address: 0x401765 pgdbg> Breakpoint at 0x401765, function abwe1, file kronrod.cpp, line 101
#101: ai = m + m + 1; pgdbg> print m
2
pgdbg>
不带调试参数"-g"的直接编译及运行结果：
[She@she-centos7 F90_call_CPP_eg1]$ make
g++ -c -std=c++ kronrod.cpp -o kronrod.o
pgf90 -c kronrod_prb.f90 -o kronrod_prb.o
pgf90 -lstdc++ -o main.exe kronrod.o kronrod_prb.o
[She@she-centos7 F90_call_CPP_eg1]$ ./main.exe
June :: PM KRONROD_PRB:
FORTRAN90 version
FORTRAN90 test program calls C++ functions. TEST01
Request KRONROD to compute the Gauss rule
of order , and the Kronrod extension of
order +=. Compare to exact data. KRONROD returns vectors of length I X WK WG 0.960491 0.104656 0.00000
0.774597 0.268488 0.555556
0.434244 0.401397 0.00000
0.00000 0.450917 0.888889 Gauss Abscissas
Exact Computed -0.774597 -0.774597
0.00000 -0.00000
0.774597 0.774597 Gauss Weights
Exact Computed 0.555556 0.555556
0.888889 0.888889
0.555556 0.555556 Gauss Kronrod Abscissas
Exact Computed -0.960491 -0.960491
-0.774597 -0.774597
-0.434244 -0.434244
0.00000 -0.00000
0.434244 0.434244
0.774597 0.774597
0.960491 0.960491 Gauss Kronrod Weights
Exact Computed 0.104656 0.104656
0.268488 0.268488
0.401397 0.401397
0.450917 0.450917
0.401397 0.401397
0.268488 0.268488
0.104656 0.104656 TEST02
Request KRONROD to compute the Gauss rule
of order , and the Kronrod extension of
order +=. KRONROD returns vectors of length I X WK WG 0.976560 0.629774E-01 0.00000
0.861136 0.170054 0.347855
0.640286 0.266798 0.00000
0.339981 0.326949 0.652145
0.00000 0.346443 0.00000 KRONROD_PRB:
Normal end of execution. June :: PM
Warning: ieee_inexact is signaling
FORTRAN STOP