admin 管理员组文章数量: 887021
2024年1月24日发(作者：专业数据恢复机构)
#include
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
//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
for( auto &elem : V){ elem = distr(eng); } }
//Print matrix V(rows, cols) storage in column-major format void print_matrix(vector
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
//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
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.167265194D+00, & 0.268488729D+00, & 0.462222905D+00, & 0.454142345D+00, & 0.462222905D+00, & 0.268488729D+00, & 0.167265194D+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.96028342D+00, & -0.77459666924148337704D+00, & -0.43424374934680255800D+00, &
-0.43424374934680255800D+00, & 0.0D+00, & 0.43424374934680255800D+00, & 0.77459666924148337704D+00, & 0.96028342D+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
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 returnendsubroutine 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)' ) ' ' 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 returnend
C++函数 #include
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,
// 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
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 << "n"; cout << "ABWE1 - Fatal error!n"; cout << " Iteration limit reached.n"; cout << " Last DELTA was " << delta << "n"; exit ( 1 ); }//// Computation of the weight.// d0 = 1.0; d1 = *x; ai = 0.0;
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;}//****************************************************************************80void 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.
// 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 << "n"; cout << "ABWE2 - Fatal error!n"; cout << " Iteration limit reached.n"; cout << " Last DELTA was " << delta << "n"; 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;}//****************************************************************************80void 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.
// 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 );
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
//****************************************************************************80double 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;}//****************************************************************************80void 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
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 PMWarning: ieee_inexact is signalingFORTRAN STOP
版权声明：本文标题：Fortran与CC++混合编程示例 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.freenas.com.cn/jishu/1706044584h499486.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论