以下例子均来自网络，只是稍作了编辑，方便今后查阅。
子目录
（一） 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(8) :: n3 end subroutine real(8) function func_d(var2d) real(8) :: var2d(2,3) end function real(8) function func_c(n3) real(8) :: n3 end function end interface integer :: n1 real :: n2 real(8) :: n3,n4 real(8) :: var2d(2,3) n1=3 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=0 var2d(1,1)=dble(n1) var2d(1,2)=dble(n2) write(*,*) "var2d(1,:): ",var2d(1,:) write(*,*) "var2d(2,:): ",var2d(2,:) n2=func_d(var2d) write(*,*) "var2d(1,:): ",var2d(1,:) write(*,*) "var2d(2,:): ",var2d(2,:) 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 [][2]); 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[3][2]) { double NN; NN=77.0; printf("%5.2f %5.2f %5.2f ",var2d[0][0],var2d[1][0],var2d[2][0]); printf("%5.2f %5.2f %5.2f ",var2d[0][1],var2d[1][1],var2d[2][1]); var2d[1][0]=NN; printf("%5.2f %5.2f %5.2f ",var2d[0][0],var2d[1][0],var2d[2][0]); printf("%5.2f %5.2f %5.2f ",var2d[0][1],var2d[1][1],var2d[2][1]); 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= 3 n2= 5.00000000 n3= 125.00000000000000 n4= 11.180339887498949 var2d(1,:): 3.0000000000000000 5.0000000000000000 0.0000000000000000 var2d(2,:): 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(1,:): 3.0000000000000000 77.000000000000000 0.0000000000000000 var2d(2,:): 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=90; t.m=0.003; int i; float e = 2.71828; char hello[32]; 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(0); }
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(*,100) HELLO(1:LENGTH),LENGTH,E,T%N,T%M 100 FORMAT(3X,A,2X,I3,4X,F7.5,2X,I3,2X,F7.5) 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 32 2.71828 90 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=3; num_float=5.0; sub_fortran_(&num_int,&num_float,&num_double); num=function_fortran_(&num_double); printf("num_int=%d num_float=%f num_double=%f num=%f",num_int,num_float,num_double,num); return 0; }
sub.f90
subroutine Sub_Fortran(NumInt,NumFloat,NumDouble) implicit none integer :: NumInt real :: NumFloat real(8) :: NumDouble NumDouble=NumFloat**NumInt end subroutine real(8) function Function_Fortran(NumDouble) implicit none real(8) :: 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 返回 1 make: *** [all] 错误 1 [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=3 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 = 0; i < rows; ++i){ for(int j = 0; j < cols; ++j){ std::cout << V[j * rows + i] << " "; } std::cout << std::endl; } std::cout << std::endl; } int main() { size_t N = 3; 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[0], &sz, &sz, &B[0], &sz, &sz, &C[0], &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 0; }
Makefile
all: gfortran -c fortran_matrix_multiply.f90 g++ -c -std=c++11 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++11 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.5/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.5 20150623 (Red Hat 4.8.5-11) (GCC) [She@she-centos7 eg2]$ gcc -v 使用内建 specs。 COLLECT_GCC=gcc COLLECT_LTO_WRAPPER=/usr/libexec/gcc/x86_64-redhat-linux/4.8.5/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.5 20150623 (Red Hat 4.8.5-11) (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 !*****************************************************************************80 ! !! MAIN is the main program for KRONROD_PRB. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 01 September 2011 ! ! 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 ( ) !*****************************************************************************80 ! !! TEST01 tests the code for the odd case N = 3. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 29 November 2010 ! ! 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 = 3 real ( c_double ) eps integer ( c_int ) i integer ( c_int ) i2 real ( c_double ) s real ( c_double ) w1(n+1) real ( c_double ) w2(n+1) real ( c_double ) :: wg(n) = (/ & 0.555555555555555555556D+00, & 0.888888888888888888889D+00, & 0.555555555555555555556D+00 /) real ( c_double ) :: wk(2*n+1) = (/ & 0.104656226026467265194D+00, & 0.268488089868333440729D+00, & 0.401397414775962222905D+00, & 0.450916538658474142345D+00, & 0.401397414775962222905D+00, & 0.268488089868333440729D+00, & 0.104656226026467265194D+00 /) real ( c_double ) x(n+1) real ( c_double ) :: xg(n) = (/ & -0.77459666924148337704D+00, & 0.0D+00, & 0.77459666924148337704D+00 /) real ( c_double ) :: xk(2*n+1) = (/ & -0.96049126870802028342D+00, & -0.77459666924148337704D+00, & -0.43424374934680255800D+00, & 0.0D+00, & 0.43424374934680255800D+00, & 0.77459666924148337704D+00, & 0.96049126870802028342D+00 /) 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+00 call kronrod ( n, eps, x, w1, w2 ) write ( *, '(a)' ) ' ' write ( *, '(a,i2)' ) ' KRONROD returns 3 vectors of length ', n + 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X WK WG' write ( *, '(a)' ) ' ' do i = 1, n + 1 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 = 1, n if ( 2 * i <= n + 1 ) then i2 = 2 * i s = -1.0D+00 else i2 = 2 * ( n + 1 ) - 2 * i s = +1.0D+00 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 = 1, n if ( 2 * i <= n + 1 ) then i2 = 2 * i else i2 = 2 * ( n + 1 ) - 2 * 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 = 1, 2 * n + 1 if ( i <= n + 1 ) then i2 = i s = -1.0D+00 else i2 = 2 * ( n + 1 ) - i s = +1.0D+00 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 = 1, 2 * n + 1 if ( i <= n + 1 ) then i2 = i else i2 = 2 * ( n + 1 ) - i end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6)' ) i, wk(i), w1(i2) end do return end subroutine test02 ( ) !*****************************************************************************80 ! !! TEST02 tests the code for the even case N = 4. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 29 November 2010 ! ! 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 = 4 real ( c_double ) eps integer ( c_int ) i integer ( c_int ) i2 real ( c_double ) s real ( c_double ) w1(n+1) real ( c_double ) w2(n+1) real ( c_double ) x(n+1) 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+00 call kronrod ( n, eps, x, w1, w2 ) write ( *, '(a)' ) ' ' write ( *, '(a,i2)' ) ' KRONROD returns 3 vectors of length ', n + 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X WK WG' write ( *, '(a)' ) ' ' do i = 1, n + 1 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 = 1; } else { ka = 0; } // // Iterative process for the computation of a Kronrod abscissa. // for ( iter = 1; iter <= 50; iter++ ) { b1 = 0.0; b2 = b[m]; yy = 4.0 * (*x) * (*x) - 2.0; d1 = 0.0; if ( even ) { ai = m + m + 1; d2 = ai * b[m]; dif = 2.0; } else { ai = m + 1; d2 = 0.0; dif = 1.0; } for ( k = 1; k <= m; k++ ) { ai = ai - dif; i = m - k + 1; b0 = b1; b1 = b2; d0 = d1; d1 = d2; b2 = yy * b1 - b0 + b[i-1]; if ( !even ) { i = i + 1; } d2 = yy * d1 - d0 + ai * b[i-1]; } 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 == 1 ) { break; } if ( r8_abs ( delta ) <= eps ) { ka = 1; } } // // Catch non-convergence. // if ( ka != 1 ) { cout << " "; cout << "ABWE1 - Fatal error! "; cout << " Iteration limit reached. "; cout << " Last DELTA was " << delta << " "; exit ( 1 ); } // // Computation of the weight. // d0 = 1.0; d1 = *x; ai = 0.0; for ( k = 2; 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 = 1; } else { ka = 0; } // // Iterative process for the computation of a Gaussian abscissa. // for ( iter = 1; iter <= 50; iter++ ) { p0 = 1.0; p1 = *x; pd0 = 0.0; pd1 = 1.0; ai = 0.0; for ( k = 2; 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 == 1 ) { break; } if ( r8_abs ( delta ) <= eps ) { ka = 1; } } // // Catch non-convergence. // if ( ka != 1 ) { cout << " "; cout << "ABWE2 - Fatal error! "; cout << " Iteration limit reached. "; cout << " Last DELTA was " << delta << " "; exit ( 1 ); } // // 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 = 1; k <= m; k++ ) { i = m - k + 1; p0 = p1; p1 = p2; p2 = yy * p1 - p0 + b[i-1]; } 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+1)/2)+1]; tau = new double[(n+1)/2]; m = ( n + 1 ) / 2; even = ( 2 * m == n ); d = 2.0; an = 0.0; for ( k = 1; k <= n; k++ ) { an = an + 1.0; d = d * an / ( an + 0.5 ); } // // Calculation of the Chebyshev coefficients of the orthogonal polynomial. // tau[0] = ( an + 2.0 ) / ( an + an + 3.0 ); b[m-1] = tau[0] - 1.0; ak = an; for ( l = 1; l < m; l++ ) { ak = ak + 2.0; tau[l] = ( ( ak - 1.0 ) * ak - an * ( an + 1.0 ) ) * ( ak + 2.0 ) * tau[l-1] / ( ak * ( ( ak + 3.0 ) * ( ak + 2.0 ) - an * ( an + 1.0 ) ) ); b[m-l-1] = tau[l]; for ( ll = 1; ll <= l; ll++ ) { b[m-l-1] = b[m-l-1] + tau[ll-1] * b[m-l+ll-1]; } } 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 ) ( 2 * n + 1 ); for ( i = 1; 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 = 1; k <= n; k = k + 2 ) { abwe1 ( n, m, eps, coef2, even, b, &xx, w1+k-1 ); w2[k-1] = 0.0; x[k-1] = 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 40 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 << " "; 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++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 [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> (1) breakpoint set at: _IO_marker::abwe1 line: "kronrod.cpp"@101 address: 0x40150F pgdbg> Breakpoint at 0x40150F, function _IO_marker::abwe1, file kronrod.cpp, line 101 #101: ai = m + m + 1; 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++11 kronrod.cpp -o kronrod.o "kronrod.cpp", line 636: 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++11 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 09 June 2017 02:51:39 PM KRONROD_PRB: FORTRAN90 version FORTRAN90 test program calls C++ functions. TEST01 Request KRONROD to compute the Gauss rule of order 3, and the Kronrod extension of order 3+4=7. Compare to exact data. KRONROD returns 3 vectors of length 4 I X WK WG 1 0.960491 0.104656 0.00000 2 0.774597 0.268488 0.555556 3 0.434244 0.401397 0.00000 4 0.00000 0.450917 0.888889 Gauss Abscissas Exact Computed 1 -0.774597 -0.774597 2 0.00000 -0.00000 3 0.774597 0.774597 Gauss Weights Exact Computed 1 0.555556 0.555556 2 0.888889 0.888889 3 0.555556 0.555556 Gauss Kronrod Abscissas Exact Computed 1 -0.960491 -0.960491 2 -0.774597 -0.774597 3 -0.434244 -0.434244 4 0.00000 -0.00000 5 0.434244 0.434244 6 0.774597 0.774597 7 0.960491 0.960491 Gauss Kronrod Weights Exact Computed 1 0.104656 0.104656 2 0.268488 0.268488 3 0.401397 0.401397 4 0.450917 0.450917 5 0.401397 0.401397 6 0.268488 0.268488 7 0.104656 0.104656 TEST02 Request KRONROD to compute the Gauss rule of order 4, and the Kronrod extension of order 4+5=9. KRONROD returns 3 vectors of length 5 I X WK WG 1 0.976560 0.629774E-01 0.00000 2 0.861136 0.170054 0.347855 3 0.640286 0.266798 0.00000 4 0.339981 0.326949 0.652145 5 0.00000 0.346443 0.00000 KRONROD_PRB: Normal end of execution. 09 June 2017 02:51:39 PM Warning: ieee_inexact is signaling FORTRAN STOP