# include # include # include # include # include # include int main ( ); void mxm_test ( double *er, double *fp, double *tm ); void mxm ( double a[], double b[], double c[], int l, int m, int n ); void cfft2d_test ( double *er, double *fp, double *tm ); void cfft2d1 ( int is, int m, int m1, int n, double complex x[], double complex w[], int ip[] ); void cfft2d2 ( int is, int m, int m1, int n, double complex x[], double complex w[], int ip[] ); void cholsky_test ( double *er, double *fp, double *tm ); void cholsky ( int ida, int nmat, int m, int n, double a[], int nrhs, int idb, double b[] ); void btrix_test ( double *er, double *fp, double *tm ); void btrix ( int js, int je, int ls, int le, int k, int jd, int kd, int ld, int md, double a[], double b[], double c[], double s[] ); void gmtry_test ( double *er, double *fp, double *tm ); void gmtry ( int nb, int nw, int nwall[], double complex proj[], double rmatrx[], double complex wall[], double xmax[], double complex zcr[] ); void emit_test ( double *er, double *fp, double *tm ); void emit ( int nb, int nvm, int nw, double cp[], double dpds[], double complex expmz[], double complex expz[], double complex force[], double gamma[], int nwall[], double ps[], double psi[], double complex refpt[], double rhs[], double rmatrx[], double rmom[], double complex wall[], double complex z[], double complex zcr[] ); void vpenta_test ( double *er, double *fp, double *tm ); void vpenta ( int jl, int ju, int kl, int ku, int nja, int njb, double a[], double b[], double c[], double d[], double e[], double f[], double x[], double y[] ); int i4_max ( int i1, int i2 ); int i4_min ( int i1, int i2 ); void r8vec_copy ( int n, double a1[], double a2[] ); void timestamp ( ); double wtime ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: nas() executes the NAS kernel benchmark. Discussion: This is a version of the NAS kernel benchmark program, whose original version was created by David Bailey, dated 17 December 1984. Each of the tests begins by filling arrays with pseudorandom values generated by the recursion: x(n+1) = 5^7 * x(n) (mod 2^30) This recursion will generate 2^28 (approx. 268 million) numbers before repeating. For this scheme to work properly, the hardware multiply operation must be correct to 47 bits of precision. Licensing: This code is distributed under the MIT license. Modified: 05 July 2015 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double er; double er_total; double fp; double fp_total; int i; char pn[9]; double rt; double tm; double tm_total; er_total = 0.0; fp_total = 0.0; tm_total = 0.0; timestamp ( ); printf ( "\n" ); printf ( "NAS:\n" ); printf ( " C version\n" ); printf ( "\n" ); printf ( " The NAS kernel benchmark program\n" ); printf ( "\n" ); printf ( " Program Error FP Ops" ); printf ( " Seconds MFLOPS\n" ); printf ( "\n" ); for ( i = 1; i <= 7; i++ ) { if ( i == 1 ) { strcpy ( pn, "BTRIX " ); btrix_test ( &er, &fp, &tm ); } else if ( i == 2 ) { strcpy ( pn, "CFFT2D " ); cfft2d_test ( &er, &fp, &tm ); } else if ( i == 3 ) { strcpy ( pn, "CHOLSKY " ); cholsky_test ( &er, &fp, &tm ); } else if ( i == 4 ) { strcpy ( pn, "EMIT " ); emit_test ( &er, &fp, &tm ); } else if ( i == 5 ) { strcpy ( pn, "GMTRY " ); gmtry_test ( &er, &fp, &tm ); } else if ( i == 6 ) { strcpy ( pn, "MXM " ); mxm_test ( &er, &fp, &tm ); } else if ( i == 7 ) { strcpy ( pn, "VPENTA " ); vpenta_test ( &er, &fp, &tm ); } rt = 1.0E-06 * fp / tm; printf ( " %s %13.4g %13.4g %10.4f %10.2f\n", pn, er, fp, tm, rt ); er_total = er_total + er; fp_total = fp_total + fp; tm_total = tm_total + tm; } strcpy ( pn, "Total " ); rt = 1.0E-06 * fp_total / tm_total; printf ( "\n" ); printf ( " %s %13.4g %13.4g %10.4f %10.2f\n", pn, er_total, fp_total, tm_total, rt ); printf ( "\n" ); printf ( "NAS:\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void mxm_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: MXM_TEST tests MXM. Licensing: This code is distributed under the MIT license. Modified: 07 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define L 256 # define M 128 # define N 64 double a[L*M]; double ans; double b[M*N]; double c[L*N]; double f7; int i; int ii; int it; int j; int l = L; int m = M; int n = N; double t; double t30; double time1; it = 100; ans = 35.2026179738722; /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t = f7 / t30; for ( j = 0; j < m; j++ ) { for ( i = 0; i < l; i++ ) { t = fmod ( f7 * t, 1.0 ); a[i+j*l] = t; } } for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { t = fmod ( f7 * t, 1.0 ); b[i+j*m] = t; } } /* Timing. */ time1 = wtime ( ); for ( ii = 1; ii <= it; ii++ ) { mxm ( a, b, c, l, m, n ); } *tm = wtime ( ) - time1; /* Results. */ *er = fabs ( ( c[18+18*l] - ans ) / ans ); *fp = ( double ) ( 2 * it * l * m * n ); return; # undef L # undef M # undef N } /******************************************************************************/ void mxm ( double a[], double b[], double c[], int l, int m, int n ) /******************************************************************************/ /* Purpose: MXM computes the matrix product C = A * B. Discussion: The function uses 4-way unrolled loops to carry out matrix multiplication. M must be a multiple of 4. Licensing: This code is distributed under the MIT license. Modified: 07 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { int i; int j; int k; for ( k = 0; k < n; k++ ) { for ( i = 0; i < l; i++ ) { c[i+k*l] = 0.0; } } for ( j = 0; j < m; j = j + 4 ) { for ( k = 0; k < n; k++ ) { for ( i = 0; i < l; i++ ) { c[i+k*l] = c[i+k*l] + a[i+ j *l] * b[j +k*m] + a[i+(j+1)*l] * b[j+1+k*m] + a[i+(j+2)*l] * b[j+2+k*m] + a[i+(j+3)*l] * b[j+3+k*m]; } } } return; } /******************************************************************************/ void cfft2d_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: CFFT2D_TEST is the test program for CFFT2D1 and CFFTD2. Licensing: This code is distributed under the MIT license. Modified: 08 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define M 128 # define M1 128 # define N 256 double ans; double f7; int i; // // IP must be dimensions 2*MAX(M,N) // int ip[2*N]; int it; int j; int k; int m = M; int m1 = M1; int n = N; double rmn; double t1; double t2; double t30; double time1; double complex w1[M]; double complex w2[N]; double complex x[M1*N]; it = 100; ans = 0.894799941219277; rmn = 1.0 / ( double ) ( m * n ); /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t2 = f7 / t30; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { t1 = fmod ( f7 * t2, 1.0 ); t2 = fmod ( f7 * t1, 1.0 ); x[i+j*m1] = t1 + t2 * I; } } cfft2d1 ( 0, m, m1, n, x, w1, ip ); cfft2d2 ( 0, m, m1, n, x, w2, ip ); /* Timing. */ time1 = wtime ( ); for ( k = 1; k <= it; k++ ) { for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { x[i+j*m1] = rmn * x[i+j*m1]; } } cfft2d1 ( 1, m, m1, n, x, w1, ip ); cfft2d2 ( 1, m, m1, n, x, w2, ip ); cfft2d2 ( -1, m, m1, n, x, w2, ip ); cfft2d1 ( -1, m, m1, n, x, w1, ip ); } *tm = wtime ( ) - time1; /* Results. */ *er = abs ( ( creal ( x[18+18*m1] ) - ans ) / ans ); *fp = ( double ) ( it * m * n ) * ( 2.0 + 10.0 * log ( ( double ) ( m * n ) ) / log ( 2.0 ) ); return; # undef M # undef M1 # undef N } /******************************************************************************/ void cfft2d1 ( int is, int m, int m1, int n, double complex x[], double complex w[], int ip[] ) /******************************************************************************/ /* Purpose: CFFT2D1 performs complex radix 2 FFT''s on the first dimension. Licensing: This code is distributed under the MIT license. Modified: 11 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double complex ct; double complex cx; int i; int i1; int i2; int ii; int im; int j; int k; int l; int m2; double pi; double t; pi = 3.141592653589793; /* If IS = 0 then initialize only. */ m2 = m / 2; if ( is == 0 ) { for ( i = 0; i < m2; i++ ) { t = 2.0 * pi * ( double ) ( i ) / ( double ) ( m ); w[i] = cos ( t ) + sin ( t ) * I; } return; } /* Perform forward or backward FFT''s according to IS = 1 or -1. */ for ( i = 0; i < m; i++ ) { ip[0+i*2] = i + 1; } l = 1; i1 = 1; for ( ; ; ) { i2 = 3 - i1; for ( j = l; j <= m2; j = j + l ) { cx = w[j-l+1-1]; if ( is < 0 ) { cx = conj ( cx ); } for ( i = j - l + 1; i <= j; i++ ) { ii = ip[i1-1+(i-1)*2]; ip[i2-1+(i+j-l-1)*2] = ii; im = ip[i1-1+(i+m2-1)*2]; ip[i2-1+(i+j-1)*2] = im; for ( k = 1; k <= n; k++ ) { ct = x[ii-1+(k-1)*m1] - x[im-1+(k-1)*m1]; x[ii-1+(k-1)*m1] = x[ii-1+(k-1)*m1] + x[im-1+(k-1)*m1]; x[im-1+(k-1)*m1] = ct * cx; } } } l = 2 * l; i1 = i2; if ( m2 < l ) { break; } } for ( i = 1; i <= m; i++ ) { ii = ip[i1-1+(i-1)*2]; if ( i < ii ) { for ( k = 1; k <= n; k++ ) { ct = x[i-1+(k-1)*m1]; x[i-1+(k-1)*m1] = x[ii-1+(k-1)*m1]; x[ii-1+(k-1)*m1] = ct; } } } return; } /******************************************************************************/ void cfft2d2 ( int is, int m, int m1, int n, double complex x[], double complex w[], int ip[] ) /******************************************************************************/ /* Purpose: CFFT2D2 performs complex radix 2 FFT''s on the second dimension. Licensing: This code is distributed under the MIT license. Modified: 08 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double complex ct; double complex cx; int i; int i1; int i2; int ii; int im; int j; int k; int l; int n2; double pi; double t; pi = 3.141592653589793; /* If IS = 0, then initialize only. */ n2 = n / 2; if ( is == 0 ) { for ( i = 0; i < n2; i++ ) { t = 2.0 * pi * ( double ) ( i ) / ( double ) ( n ); w[i] = cos ( t ) + sin ( t ) * I; } return; } /* Perform forward or backward FFT''s according to IS = 1 or -1. */ for ( i = 0; i < n; i++ ) { ip[0+i*2] = i + 1; } l = 1; i1 = 1; for ( ; ; ) { i2 = 3 - i1; for ( j = l; j <= n2; j = j + l ) { cx = w[j-l+1-1]; if ( is < 0 ) { cx = conj ( cx ); } for ( i = j - l + 1; i <= j; i++ ) { ii = ip[i1-1+(i-1)*2]; ip[i2-1+(i+j-l-1)*2] = ii; im = ip[i1-1+(i+n2-1)*2]; ip[i2-1+(i+j-1)*2] = im; for ( k = 1; k <= m; k++ ) { ct = x[k-1+(ii-1)*m1] - x[k-1+(im-1)*m1]; x[k-1+(ii-1)*m1] = x[k-1+(ii-1)*m1] + x[k-1+(im-1)*m1]; x[k-1+(im-1)*m1] = ct * cx; } } } l = 2 * l; i1 = i2; if ( n2 < l ) { break; } } for ( i = 1; i <= n; i++ ) { ii = ip[i1-1+(i-1)*2]; if ( i < ii ) { for ( k = 1; k <= m; k++ ) { ct = x[k-1+(i-1)*m1]; x[k-1+(i-1)*m1] = x[k-1+(ii-1)*m1]; x[k-1+(ii-1)*m1] = ct; } } } return; } /******************************************************************************/ void cholsky_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: CHOLSKY_TEST tests CHOLSKY. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define IDA 250 # define M 4 # define N 40 # define NMAT 250 # define NRHS 3 double a[(IDA+1)*(M+1)*(N+1)]; double ans; double ax[(IDA+1)*(M+1)*(N+1)]; double b[(NRHS+1)*(NMAT+1)*(N+1)]; double bx[(NRHS+1)*(NMAT+1)*(N+1)]; double f7; int i; int i1; int ida = IDA; int it; int j; int k; int la; int lb; int m = M; int n = N; int nmat = NMAT; int nrhs = NRHS; double t; double t30; double time1; it = 200; ans = 5177.88531774562; la = ( ida + 1 ) * ( m + 1 ) * ( n + 1 ); lb = ( nrhs + 1 ) * ( nmat + 1 ) * ( n + 1 ); /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t = f7 / t30; i1 = 0; for ( k = 0; k <= n; k++ ) { for ( j = -m; j <= 0; j++ ) { for ( i = 0; i <= ida; i++ ) { t = fmod ( f7 * t, 1.0 ); ax[i1] = t; i1 = i1 + 1; } } } i1 = 0; for ( k = 0; k <= n; k++ ) { for ( j = 0; j <= nmat; j++ ) { for ( i = 0; i <= nrhs; i++ ) { t = fmod ( f7 * t, 1.0 ); bx[i1] = t; i1 = i1 + 1; } } } /* Timing. */ time1 = wtime ( ); for ( j = 1; j <= it; j++ ) { r8vec_copy ( la, ax, a ); r8vec_copy ( lb, bx, b ); cholsky ( ida, nmat, m, n, a, nrhs, ida, b ); } *tm = wtime ( ) - time1; /* Results. */ i1 = 1+19*(nrhs+1)+19*(nrhs+1)*(nmat+1); *er = fabs ( ( b[i1] - ans ) / ans ); *fp = ( double ) ( it * ( nmat + 1 ) * 4403 ); return; # undef IDA # undef M # undef N # undef NMAT # undef NRHS } /******************************************************************************/ void cholsky ( int ida, int nmat, int m, int n, double a[], int nrhs, int idb, double b[] ) /******************************************************************************/ /* Purpose: CHOLSKY carries out Cholesky decomposition and back substitution. Discussion: The Cholesky decomposition is performed on a set of input matrices which are provided as a single three-dimensional array. Licensing: This code is distributed under the MIT license. Modified: 07 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double eps; double epss[257]; int i; int i0; int i1; int i2; int i3; int j; int jj; int k; int l; eps = 1.0E-13; /* Cholesky decomposition. */ for ( j = 0; j <= n; j++ ) { i0 = i4_max ( -m, -j ); /* Off diagonal elements. */ for ( i = i0; i <= -1; i++ ) { for ( jj = i0 - i; jj <= -1; jj++ ) { for ( l = 0; l <= nmat; l++ ) { i1 = l + (i+m) * (ida+1) + j * (ida+1)*(m+1); i2 = l + (jj+m) * (ida+1) + (i+j) * (ida+1)*(m+1); i3 = l + (i+jj+m) * (ida+1) + j * (ida+1)*(m+1); a[i1] = a[i1] - a[i2] * a[i3]; } } for ( l = 0; l <= nmat; l++ ) { i1 = l + (i+m) * (ida+1) + j * (ida+1)*(m+1); i2 = l + m * (ida+1) + (i+j) * (ida+1)*(m+1); a[i1] = a[i1] * a[i2]; } } /* Store inverse of diagonal elements. */ for ( l = 0; l <= nmat; l++ ) { i1 = l + + m * (ida+1) + j * (ida+1)*(m+1); epss[l] = eps * a[i1]; } for ( jj = i0; jj <= -1; jj++ ) { for ( l = 0; l <= nmat; l++ ) { i1 = l + m * (ida+1) + j * (ida+1)*(m+1); i2 = l + (jj+m) * (ida+1) + j * (ida+1)*(m+1); a[i1] = a[i1] - pow ( a[i2], 2 ); } } for ( l = 0; l <= nmat; l++ ) { i1 = l + m * (ida+1) + j * (ida+1)*(m+1); a[i1] = 1.0 / sqrt ( fabs ( epss[l] + a[i1] ) ); } } /* Solution. */ for ( i = 0; i <= nrhs; i++ ) { for ( k = 0; k <= n; k++ ) { for ( l = 0; l <= nmat; l++ ) { i1 = i+l*(nrhs+1)+k*(nrhs+1)*(idb+1); i2 = l + m * (ida+1) + k * (ida+1)*(m+1); b[i1] = b[i1] * a[i2]; } for ( jj = 1; jj <= i4_min ( m, n - k ); jj++ ) { for ( l = 0; l <= nmat; l++ ) { i1 = i+l*(nrhs+1)+(k+jj)*(nrhs+1)*(idb+1); i2 = l + (-jj+m) * (ida+1) + (k+jj) * (ida+1)*(m+1); i3 = i+l*(nrhs+1)+k*(nrhs+1)*(idb+1); b[i1] = b[i1] - a[i2] * b[i3]; } } } for ( k = n; 0 <= k; k-- ) { for ( l = 0; l <= nmat; l++ ) { i1 = i + l*(nrhs+1) + k*(nrhs+1)*(idb+1); i2 = l + m * (ida+1) + k * (ida+1)*(m+1); b[i1] = b[i1] * a[i2]; } for ( jj = 1; jj <= i4_min ( m, k ); jj++ ) { for ( l = 0; l <= nmat; l++ ) { i1 = i+l*(nrhs+1)+(k-jj)*(nrhs+1)*(idb+1); i2 = l + (-jj+m) * (ida+1) + k * (ida+1)*(m+1); i3 = i+l*(nrhs+1)+k*(nrhs+1)*(idb+1); b[i1] = b[i1] - a[i2] * b[i3]; } } } } return; } /******************************************************************************/ void btrix_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: BTRIX_TEST tests BTRIX. Licensing: This code is distributed under the MIT license. Modified: 08 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define JD 30 # define KD 30 # define LD 30 # define MD 30 double a[5*5*MD*MD]; double ans; double b[5*5*MD*MD]; double bx[5*5*MD*MD]; double c[5*5*MD*MD]; double f7; int i; int ii; int it; int j; int jd = JD; int je; int js; int k; int kd = KD; int l; int ld = LD; int le; int ls; int md = MD; int nb; int ns; double s[JD*KD*LD*5]; double sx[JD*KD*LD*5]; double t; double t30; double time1; js = 1; je = 28; ls = 1; le = 28; it = 20; ans = -0.286282658663962; nb = 25 * md * md; ns = jd * kd * ld * 5; /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t = f7 / t30; for ( l = 0; l < md; l++ ) { for ( k = 0; k < md; k++ ) { for ( j = 0; j < 5; j++ ) { for ( i = 0; i < 5; i++ ) { t = fmod ( f7 * t, 1.0 ); a[i+j*5+k*5*5+l*5*5*md] = t; t = fmod ( f7 * t, 1.0 ); bx[i+j*5+k*5*5+l*5*5*md] = t; t = fmod ( f7 * t, 1.0 ); c[i+j*5+k*5*5+l*5*5*md] = t; } } } } for ( l = 0; l < 5; l++ ) { for ( k = 0; k < ld; k++ ) { for ( j = 0; j < kd; j++ ) { for ( i = 0; i < jd; i++ ) { t = fmod ( f7 * t, 1.0 ); sx[i+j*jd+k*jd*kd+l*jd*kd*ld] = t; } } } } /* Timing. */ time1 = wtime ( ); for ( ii = 1; ii <= it; ii++ ) { r8vec_copy ( ns, sx, s ); for ( k = 1; k <= kd; k++ ) { r8vec_copy ( nb, bx, b ); btrix ( js, je, ls, le, k, jd, kd, ld, md, a, b, c, s ); } } *tm = wtime ( ) - time1; /* Results. */ *er = fabs ( ( s[18+18*jd+18*jd*kd+0*jd*kd*ld] - ans ) / ans ); *fp = ( double ) ( it * md * ( le - 1 ) * 19165 ); return; # undef JD # undef KD # undef LD # undef MD } /******************************************************************************/ void btrix ( int js, int je, int ls, int le, int k, int jd, int kd, int ld, int md, double a[], double b[], double c[], double s[] ) /******************************************************************************/ /* Purpose: BTRIX is a block tridiagonal solver in one direction. Discussion: The array has four dimensions. The routine solves along the "J" index. Licensing: This code is distributed under the MIT license. Modified: 08 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define MD 30 double c1; double c2; double c3; double c4; double c5; double d1; double d2; double d3; double d4; double d5; int j; int jem1; int l; double l11[MD]; double l21[MD]; double l31[MD]; double l41[MD]; double l51[MD]; double l22[MD]; double l32[MD]; double l33[MD]; double l42[MD]; double l43[MD]; double l44[MD]; double l52[MD]; double l53[MD]; double l54[MD]; double l55[MD]; int m; int n; double u12[MD]; double u13[MD]; double u14[MD]; double u15[MD]; double u23[MD]; double u24[MD]; double u25[MD]; double u34[MD]; double u35[MD]; double u45[MD]; /* Part 1. Forward block sweep. */ for ( j = js; j <= je; j++ ) { /* Step 1. Construct L(I) in B. */ if ( j != js ) { for ( m = 0; m < 5; m++ ) { for ( n = 0; n < 5; n++ ) { for ( l = ls; l <= le; l++ ) { b[m+n*5+j*5*5+l*5*5*md] = b[m+n*5+j*5*5+l*5*5*md] - a[m+0*5+j*5*5+l*5*5*md] * b[0+n*5+(j-1)*5*5+l*5*5*md] - a[m+1*5+j*5*5+l*5*5*md] * b[1+n*5+(j-1)*5*5+l*5*5*md] - a[m+2*5+j*5*5+l*5*5*md] * b[2+n*5+(j-1)*5*5+l*5*5*md] - a[m+3*5+j*5*5+l*5*5*md] * b[3+n*5+(j-1)*5*5+l*5*5*md] - a[m+4*5+j*5*5+l*5*5*md] * b[4+n*5+(j-1)*5*5+l*5*5*md]; } } } } /* Step 2. Compute L inverse. A. Decompose L(I) into L and U. */ for ( l = ls; l <= le; l++ ) { l11[l] = 1.0 / b[0+0*5+j*5*5+l*5*5*md]; u12[l] = b[0+1*5+j*5*5+l*5*5*md] * l11[l]; u13[l] = b[0+2*5+j*5*5+l*5*5*md] * l11[l]; u14[l] = b[0+3*5+j*5*5+l*5*5*md] * l11[l]; u15[l] = b[0+4*5+j*5*5+l*5*5*md] * l11[l]; l21[l] = b[1+0*5+j*5*5+l*5*5*md]; l22[l] = 1.0 / ( b[1+1*5+j*5*5+l*5*5*md] - l21[l] * u12[l] ); u23[l] = ( b[1+2*5+j*5*5+l*5*5*md] - l21[l] * u13[l] ) * l22[l]; u24[l] = ( b[1+3*5+j*5*5+l*5*5*md] - l21[l] * u14[l] ) * l22[l]; u25[l] = ( b[1+4*5+j*5*5+l*5*5*md] - l21[l] * u15[l] ) * l22[l]; l31[l] = b[2+0*5+j*5*5+l*5*5*md]; l32[l] = b[2+1*5+j*5*5+l*5*5*md] - l31[l] * u12[l]; l33[l] = 1.0 / ( b[2+2*5+j*5*5+l*5*5*md] - l31[l] * u13[l] - l32[l] * u23[l] ); u34[l] = ( b[2+3*5+j*5*5+l*5*5*md] - l31[l] * u14[l] - l32[l] * u24[l] ) * l33[l]; u35[l] = ( b[2+4*5+j*5*5+l*5*5*md] - l31[l] * u15[l] - l32[l] * u25[l] ) * l33[l]; } for ( l = ls; l <= le; l++ ) { l41[l] = b[3+0*5+j*5*5+l*5*5*md]; l42[l] = b[3+1*5+j*5*5+l*5*5*md] - l41[l] * u12[l]; l43[l] = b[3+2*5+j*5*5+l*5*5*md] - l41[l] * u13[l] - l42[l] * u23[l]; l44[l] = 1.0 / ( b[3+3*5+j*5*5+l*5*5*md] - l41[l] * u14[l] - l42[l] * u24[l] - l43[l] * u34[l] ); u45[l] = ( b[3+4*5+j*5*5+l*5*5*md] - l41[l] * u15[l] - l42[l] * u25[l] - l43[l] * u35[l] ) * l44[l]; l51[l] = b[4+0*5+j*5*5+l*5*5*md]; l52[l] = b[4+1*5+j*5*5+l*5*5*md] - l51[l] * u12[l]; l53[l] = b[4+2*5+j*5*5+l*5*5*md] - l51[l] * u13[l] - l52[l] * u23[l]; l54[l] = b[4+3*5+j*5*5+l*5*5*md] - l51[l] * u14[l] - l52[l] * u24[l] - l53[l] * u34[l]; l55[l] = 1.0 / ( b[4+4*5+j*5*5+l*5*5*md] - l51[l] * u15[l] - l52[l] * u25[l] - l53[l] * u35[l] - l54[l] * u45[l] ); } /* Step 3. Solve for intermediate vector. A. Construct the right hand side. */ if ( j != js ) { for ( m = 0; m < 5; m++ ) { for ( l = ls; l <= le; l++ ) { s[j+k*jd+l*jd*kd+m*jd*kd*ld] = s[j+k*jd+l*jd*kd+m*jd*kd*ld] - a[m+0*5+j*5*5+l*5*5*md] * s[j-1+k*jd+l*jd*kd+0*jd*kd*ld] - a[m+1*5+j*5*5+l*5*5*md] * s[j-1+k*jd+l*jd*kd+1*jd*kd*ld] - a[m+2*5+j*5*5+l*5*5*md] * s[j-1+k*jd+l*jd*kd+2*jd*kd*ld] - a[m+3*5+j*5*5+l*5*5*md] * s[j-1+k*jd+l*jd*kd+3*jd*kd*ld] - a[m+4*5+j*5*5+l*5*5*md] * s[j-1+k*jd+l*jd*kd+4*jd*kd*ld]; } } } /* B. Intermediate vector. Forward substitution. */ for ( l = ls; l <= le; l++ ) { d1 = s[j+k*jd+l*jd*kd+0*jd*kd*ld] * l11[l]; d2 = ( s[j+k*jd+l*jd*kd+1*jd*kd*ld] - l21[l] * d1 ) * l22[l]; d3 = ( s[j+k*jd+l*jd*kd+2*jd*kd*ld] - l31[l] * d1 - l32[l] * d2 ) * l33[l]; d4 = ( s[j+k*jd+l*jd*kd+3*jd*kd*ld] - l41[l] * d1 - l42[l] * d2 - l43[l] * d3 ) * l44[l]; d5 = ( s[j+k*jd+l*jd*kd+4*jd*kd*ld] - l51[l] * d1 - l52[l] * d2 - l53[l] * d3 - l54[l] * d4 ) * l55[l]; /* Backward substitution. */ s[j+k*jd+l*jd*kd+4*jd*kd*ld] = d5; s[j+k*jd+l*jd*kd+3*jd*kd*ld] = d4 - u45[l] * d5; s[j+k*jd+l*jd*kd+2*jd*kd*ld] = d3 - u34[l] * s[j+k*jd+l*jd*kd+3*jd*kd*ld] - u35[l] * d5; s[j+k*jd+l*jd*kd+1*jd*kd*ld] = d2 - u23[l] * s[j+k*jd+l*jd*kd+2*jd*kd*ld] - u24[l] * s[j+k*jd+l*jd*kd+3*jd*kd*ld] - u25[l] * d5; s[j+k*jd+l*jd*kd+0*jd*kd*ld] = d1 - u12[l] * s[j+k*jd+l*jd*kd+1*jd*kd*ld] - u13[l] * s[j+k*jd+l*jd*kd+2*jd*kd*ld] - u14[l] * s[j+k*jd+l*jd*kd+3*jd*kd*ld] - u15[l] * d5; } /* Step 4. Construct U(I) = inverse(L(I))*C(I+1) by columns and store in B. */ if ( j != je ) { for ( n = 0; n < 5; n++ ) { for ( l = ls; l <= le; l++ ) { /* Forward substitution. */ c1 = c[0+n*5+j*5*5+l*5*5*md] * l11[l]; c2 = ( c[1+n*5+j*5*5+l*5*5*md] - l21[l] * c1 ) * l22[l]; c3 = ( c[2+n*5+j*5*5+l*5*5*md] - l31[l] * c1 - l32[l] * c2 ) * l33[l]; c4 = ( c[3+n*5+j*5*5+l*5*5*md] - l41[l] * c1 - l42[l] * c2 - l43[l] * c3 ) * l44[l]; c5 = ( c[4+n*5+j*5*5+l*5*5*md] - l51[l] * c1 - l52[l] * c2 - l53[l] * c3 - l54[l] * c4 ) * l55[l]; /* Backward substitution. */ b[4+n*5+j*5*5+l*5*5*md] = c5; b[3+n*5+j*5*5+l*5*5*md] = c4 - u45[l] * c5; b[2+n*5+j*5*5+l*5*5*md] = c3 - u34[l] * b[3+n*5+j*5*5+l*5*5*md] - u35[l] * c5; b[1+n*5+j*5*5+l*5*5*md] = c2 - u23[l] * b[2+n*5+j*5*5+l*5*5*md] - u24[l] * b[3+n*5+j*5*5+l*5*5*md] - u25[l] * c5; b[0+n*5+j*5*5+l*5*5*md] = c1 - u12[l] * b[1+n*5+j*5*5+l*5*5*md] - u13[l] * b[2+n*5+j*5*5+l*5*5*md] - u14[l] * b[3+n*5+j*5*5+l*5*5*md] - u15[l] * c5; } } } } /* Part 2. Backward block sweep. */ jem1 = je - 1; for ( j = jem1; js <= j; j-- ) { for ( m = 0; m < 5; m++ ) { for ( l = ls; l <= le; l++ ) { s[j+k*jd+l*jd*kd+m*jd*kd*ld] = s[j+k*jd+l*jd*kd+m*jd*kd*ld] - b[m+0*5+j*5*5+l*5*5*md] * s[j+1+k*jd+l*jd*kd+0*jd*kd*ld] - b[m+1*5+j*5*5+l*5*5*md] * s[j+1+k*jd+l*jd*kd+1*jd*kd*ld] - b[m+2*5+j*5*5+l*5*5*md] * s[j+1+k*jd+l*jd*kd+2*jd*kd*ld] - b[m+3*5+j*5*5+l*5*5*md] * s[j+1+k*jd+l*jd*kd+3*jd*kd*ld] - b[m+4*5+j*5*5+l*5*5*md] * s[j+1+k*jd+l*jd*kd+4*jd*kd*ld]; } } } return; # undef MD } /******************************************************************************/ void gmtry_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: GMTRY_TEST tests GMTRY. Licensing: This code is distributed under the MIT license. Modified: 09 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define NB 5 # define NW 100 double ans; double f7; int i; int it; int j; int nb = NB; int nw = NW; int nwall[NB]; double complex proj[NW*NB]; double rmatrx[NW*NB*NW*NB]; double t1; double t2; double t30; double time1; double complex wall[NW*NB]; double xmax[NB]; double complex zcr[NW*NB]; it = 2; ans = -2.57754233214174; /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t2 = f7 / t30; for ( j = 1; j <= nb; j++ ) { nwall[j-1] = nw; } for ( j = 1; j <= nb; j++ ) { for ( i = 1; i <= nw; i++ ) { t1 = fmod ( f7 * t2, 1.0 ); t2 = fmod ( f7 * t1, 1.0 ); wall[i-1+(j-1)*NW] = t1 + t2 * I; } } /* Timing. */ time1 = wtime ( ); for ( i = 1; i <= it; i++ ) { gmtry ( nb, nw, nwall, proj, rmatrx, wall, xmax, zcr ); } *tm = wtime ( ) - time1; /* Results. */ *er = fabs ( ( rmatrx[18+18*nw*nb] - ans ) / ans ); *fp = ( double ) ( it ) * ( ( double ) ( 120 * ( nb * nw * nb * nw ) ) + 0.666 * ( double ) ( nb * nw * nb * nw * nb * nw ) ); return; # undef NB # undef NW } /******************************************************************************/ void gmtry ( int nb, int nw, int nwall[], double complex proj[], double rmatrx[], double complex wall[], double xmax[], double complex zcr[] ) /******************************************************************************/ /* Purpose: GMTRY computes solid-related arrays. Discussion: This function was extracted from a vortex method program. It sets up arrays needed for the computation, and performs Gauss elimination on the matrix of wall-influence coefficients. Licensing: This code is distributed under the MIT license. Modified: 09 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double arcl; double dum; int i; int i0; int j; int j0; int k; int k1; int k2; int kron; int l; int l1; int l2; int matdim; double period; double pi; double pidp; double r0; double sig2; double sigma; double ylimit; double ymax; double ymin; double complex z1; double complex zi; double complex zz; pi = 3.141592653589793; period = 3.0; /* Compute arclength. */ matdim = 0; arcl = 0.0; ymin = 1.0E+30; ymax = -1.0E+30; pidp = pi / period; for ( l = 1; l <= nb; l++ ) { matdim = matdim + nwall[l-1]; for ( k = 1; k <= nwall[l-1]; k++ ) { k2 = 1 + ( k % nwall[l-1] ); arcl = arcl + cabs ( wall[k-1+(l-1)*nw] - wall[k2-1+(l-1)*nw] ); } } /* Compute core radius. */ r0 = 0.5 * arcl / ( double ) ( matdim ); sigma = r0 / 2.0; /* Define creation points. */ for ( l = 1; l <= nb; l++ ) { for ( k = 1; k <= nwall[l-1]; k++ ) { k1 = 1 + ( k + nwall[l-1] - 2 ) % ( nwall[l-1] ); k2 = 1 + k % nwall[l-1]; zz = wall[k1-1+(l-1)*nw] - wall[k2-1+(l-1)*nw]; zcr[k-1+(l-1)*nw] = wall[k-1+(l-1)*nw] + r0 * I / cabs ( zz ) * zz; } /* Check that wall and creation points are not crossed due to too sharp a concave kink or an error in defining the body. Also find highest, lowest and right-most point. */ xmax[l-1] = creal ( zcr[0+(l-1)*nw] ); for ( k = 1; k <= nwall[l-1]; k++ ) { ymin = fmin ( ymin, cimag ( zcr[k-1+(l-1)*nw] ) ); ymax = fmax ( ymax, cimag ( zcr[k-1+(l-1)*nw] ) ); xmax[l-1] = fmax ( xmax[l-1], creal ( zcr[k-1+(l-1)*nw] ) ); } } /* The "main period" will be between ylimit and ylimit + period. */ ylimit = ( ymin - period + ymax ) / 2.0; /* Project creation points into main period. This is technical. */ for ( l = 1; l <= nb; l++ ) { for ( k = 1; k <= nwall[l-1]; k++ ) { proj[k-1+(l-1)*nw] = zcr[k-1+(l-1)*nw] - period * I * ( ( int ) ( 5.0 + ( cimag ( zcr[k-1+(l-1)*nw] ) - ylimit ) / period ) - 5.0 ); } } /* Compute matrix. */ sig2 = pow ( 2.0 * pidp * sigma, 2 ); i0 = 0; for ( l1 = 1; l1 <= nb; l1++ ) { j0 = 0; for ( l2 = 1; l2 <= nb; l2++ ) { if ( l1 == l2 ) { kron = 1; } else { kron = 0; } for ( j = 1; j <= nwall[l2-1]; j++ ) { rmatrx[i0+(j0+j-1)*nw*nb] = kron; z1 = cexp ( ( wall[0+(l1-1)*nw] - zcr[j-1+(l2-1)*nw] ) * pidp ); z1 = z1 - 1.0 / z1; dum = sig2 + pow ( creal ( z1 ), 2 ) + pow ( cimag ( z1 ), 2 ); for ( i = 2; i <= nwall[l1-1]; i++ ) { zi = cexp ( ( wall[i-1+(l1-1)*nw] - zcr[j-1+(l2-1)*nw] ) * pidp ); zz = zi - 1.0 / zi; rmatrx[i0+i-1+(j0+j-1)*nw*nb] = -0.25 / pi * log ( dum / ( sig2 + pow ( creal ( zz ), 2 ) + pow ( cimag ( zz ), 2 ) ) ); } } j0 = j0 + nwall[l2-1]; } i0 = i0 + nwall[l1-1]; } /* Gauss elimination. */ for ( i = 1; i <= matdim; i++ ) { rmatrx[i-1+(i-1)*nw*nb] = 1.0 / rmatrx[i-1+(i-1)*nw*nb]; for ( j = i + 1; j <= matdim; j++ ) { rmatrx[j-1+(i-1)*nw*nb] = rmatrx[j-1+(i-1)*nw*nb] * rmatrx[i-1+(i-1)*nw*nb]; for ( k = i + 1; k <= matdim; k++ ) { rmatrx[j-1+(k-1)*nw*nb] = rmatrx[j-1+(k-1)*nw*nb] - rmatrx[j-1+(i-1)*nw*nb] * rmatrx[i-1+(k-1)*nw*nb]; } } } return; } /******************************************************************************/ void emit_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: EMIT_TEST tests EMIT. Licensing: This code is distributed under the MIT license. Modified: 11 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define NB 5 # define NV 1000 # define NVM 1500 # define NW 100 double ans; double cp[NW*NB]; double dpds[NW*NB]; double complex expmz[NVM]; double complex expz[NVM]; double f7; double complex force[NB]; double gamma[NVM]; int i; int it; int j; int nb = NB; int nv = NV; int nvm = NVM; int nw = NW; int nwall[NB]; double ps[NVM]; double psi[NW]; double complex refpt[NB]; double rhs[NW*NB]; double rmatrx[NW*NB*NW*NB]; double rmom[NB]; double t1; double t2; double t30; double time1; double complex wall[NW*NB]; double complex z[NVM]; double complex zcr[NW*NB]; it = 10; ans = 6.0088546832072; /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t2 = f7 / t30; for ( j = 0; j < nb; j++ ) { nwall[j] = nw; refpt[j] = 0.0; force[j] = 0.0; rmom[j] = 0.0; for ( i = 0; i < nw; i++ ) { t1 = fmod ( f7 * t2, 1.0 ); t2 = fmod ( f7 * t1, 1.0 ); wall[i+j*nw] = t1 + t2 * I; t1 = fmod ( f7 * t2, 1.0 ); t2 = fmod ( f7 * t1, 1.0 ); zcr[i+j*nw] = t1 + t2 * I; dpds[i+j*nw] = 0.0; } } for ( j = 0; j < nw * nb; j++ ) { rmatrx[j+j*nw*nb] = 1.0; for ( i = 0; i < j; i++ ) { t2 = fmod ( f7 * t2, 1.0 ); rmatrx[i+j*nw*nb] = 0.001 * t2; rmatrx[j+i*nw*nb] = 0.001 * t2; } } for ( i = 0; i < nvm; i++ ) { t1 = fmod ( f7 * t2, 1.0 ); t2 = fmod ( f7 * t1, 1.0 ); z[i] = t1 + t2 * I; t2 = fmod ( f7 * t2, 1.0 ); gamma[i] = t2; } /* Timing. */ time1 = wtime ( ); for ( i = 1; i <= it; i++ ) { emit ( nb, nvm, nw, cp, dpds, expmz, expz, force, gamma, nwall, ps, psi, refpt, rhs, rmatrx, rmom, wall, z, zcr ); } *tm = wtime ( ) - time1; /* Results. */ *er = fabs ( ( rhs[18] - ans ) / ans ); *fp = ( double ) ( it * ( 56 * nv + nb * nw * ( 97 + 44 * nv + 2 * nb * nw ) ) ); return; # undef NB # undef NV # undef NVM # undef NW } /******************************************************************************/ void emit ( int nb, int nvm, int nw, double cp[], double dpds[], double complex expmz[], double complex expz[], double complex force[], double gamma[], int nwall[], double ps[], double psi[], double complex refpt[], double rhs[], double rmatrx[], double rmom[], double complex wall[], double complex z[], double complex zcr[] ) /******************************************************************************/ /* Purpose: EMIT creates new vortices according to certain boundary conditions. Discussion: This function was extracted from a vortex method program. It emits new vortices to satisfy the boundary condition. It also finishes computing pressure, forces, and other quantities. Licensing: This code is distributed under the MIT license. Modified: 10 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { double chord; double cpm; double cupst; double delt; double complex dum3; double complex expmwk; double complex expwkl; int i; int i0; int j; int k; int k1; int k2; int l; int matdim; int nv; double period; double pi; double pidp; double sig2; double sps; double u0; double complex uupstr; period = 3.0; sig2 = 3.0; u0 = 4.0; matdim = 500; delt = 1.0; chord = 5.0; pi = 3.141592653589793; uupstr = 3.0 + 4.0 * I; /* Store exp(z(i)) and exp(-z(i)) to reduce work in inner loop. Note that the NV used here is a variable, whereas the NV in the calling program is a constant. They are separate quantities. */ nv = 1000; pidp = pi / period; for ( i = 0; i < nv; i++ ) { expz[i] = cexp ( z[i] * pidp ); expmz[i] = 1.0 / expz[i]; } i0 = 0; cupst = pow ( creal ( uupstr ), 2 ) + pow ( cimag ( uupstr ), 2 ); for ( l = 0; l < nb; l++ ) { for ( k = 0; k < nwall[l]; k++ ) { expwkl = cexp ( wall[k+l*nw] * pidp ); expmwk = 1.0 / expwkl; sps = 0.0; for ( i = 0; i < nv; i++ ) { dum3 = expz[i] * expmwk - expwkl * expmz[i]; ps[i] = gamma[i] * log ( pow ( creal ( dum3 ), 2 ) + pow ( cimag ( dum3 ), 2 ) + sig2 ); sps = sps + ps[i]; } psi[k] = cimag ( wall[k+l*nw] * conj ( uupstr + u0 * I ) ) - sps * 0.25 / pi; } /* Compute the right-hand side. */ for ( k = 0; k < nwall[l]; k++ ) { rhs[i0+k] = psi[k] - psi[0]; } i0 = i0 + nwall[l]; } /* Solve the system. */ for ( i = 0; i < matdim; i++ ) { for ( j = i + 1; j < matdim; j++ ) { rhs[j] = rhs[j] - rmatrx[j+i*nw*nb] * rhs[i]; } } for ( i = matdim - 1; 0 <= i; i-- ) { rhs[i] = rmatrx[i+i*nw*nb] * rhs[i]; for ( j = 0; j < i; j++ ) { rhs[j] = rhs[j] - rmatrx[j+i*nw*nb] * rhs[i]; } } /* Create new vortices. */ i0 = 0; for ( l = 0; l < nb; l++ ) { for ( k = 0; k < nwall[l]; k++ ) { /* Put the new vortex at the end of the array. */ z[nv] = zcr[k+l*nw]; gamma[nv] = rhs[i0+k]; /* Record the gain of linear and angular momentum. */ force[l] = force[l] + gamma[nv] * z[nv]; rmom[l] = rmom[l] + gamma[nv] * ( pow ( creal ( z[nv] - refpt[l] ), 2 ) + pow ( cimag ( z[nv] - refpt[l] ), 2 ) ); dpds[k+l*nw] = dpds[k+l*nw] - gamma[nv]; nv = nv + 1; } /* Filter and integrate pressure gradient to get pressure. */ cp[0+l*nw] = 0.0; cpm = -1.0E+30; for ( k = 1; k < nwall[l]; k++ ) { k1 = k % nwall[l]; k2 = ( k + nwall[l] - 3 ) % nwall[l]; cp[k+l*nw] = cp[k-1+l*nw] + ( 3.0 * ( dpds[k+l*nw] + dpds[k-1+l*nw] ) + dpds[k1+l*nw] + dpds[k2+l*nw] ) / ( 4.0 * delt * cupst ); cpm = fmax ( cpm, cp[k+l*nw] ); } /* Normalize the pressure. */ for ( k = 0; k < nwall[l]; k++ ) { cp[k+l*nw] = cp[k+l*nw] - cpm; } /* Finish computing force and moment, as time rate of change of linear and angular momentum. */ force[l] = force[l] * 2.0 * I / ( delt * chord * cupst ); rmom[l] = rmom[l] * 2.0 / ( delt * chord * chord * cupst ); i0 = i0 + nwall[l]; } return; } /******************************************************************************/ void vpenta_test ( double *er, double *fp, double *tm ) /******************************************************************************/ /* Purpose: VPENA_TEST tests VPENTA. Licensing: This code is distributed under the MIT license. Modified: 09 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { # define JL 0 # define JU 127 # define KL 0 # define KU 127 # define NJA 128 # define NJB 128 double a[NJA*NJB]; double ans; double b[NJA*NJB]; double c[NJA*NJB]; double d[NJA*NJB]; double e[NJA*NJB]; double f[NJA*NJB*3]; double f7; double fx[NJA*NJB*3]; int i; int it; int j; int jl = JL; int ju = JU; int k; int kl = KL; int ku = KU; int lf; int nja = NJA; int njb = NJB; double t; double t30; double time1; double x[NJA*NJB]; double y[NJA*NJB]; it = 400; ans = -0.354649411858726; lf = nja * njb * 3; /* Random initialization. */ f7 = 78125.0; t30 = 1073741824.0; t = f7 / t30; for ( j = kl; j <= ku; j++ ) { for ( i = jl; i <= ju; i++ ) { t = fmod ( f7 * t, 1.0 ); a[i+j*nja] = t; t = fmod ( f7 * t, 1.0 ); b[i+j*nja] = t; t = fmod ( f7 * t, 1.0 ); c[i+j*nja] = t; t = fmod ( f7 * t, 1.0 ); d[i+j*nja] = t; t = fmod ( f7 * t, 1.0 ); e[i+j*nja] = t; for ( k = 0; k < 3; k++ ) { t = fmod ( f7 * t, 1.0 ); fx[i+j*nja+k*nja*njb] = t; } } } /* Timing. */ time1 = wtime ( ); for ( i = 1; i <= it; i++ ) { r8vec_copy ( lf, fx, f ); vpenta ( jl, ju, kl, ku, nja, njb, a, b, c, d, e, f, x, y ); } *tm = wtime ( ) - time1; /* Results. */ *er = fabs ( ( f[18+18*nja+0*nja*njb] - ans ) / ans ); *fp = ( double ) ( it * ku * ( 40 * ku - 53 ) ); return; # undef JL # undef JU # undef KL # undef KU # undef NJA # undef NJB } /******************************************************************************/ void vpenta ( int jl, int ju, int kl, int ku, int nja, int njb, double a[], double b[], double c[], double d[], double e[], double f[], double x[], double y[] ) /******************************************************************************/ /* Purpose: VPENTA inverts 3 pentadiagonal systems simultaneously. Licensing: This code is distributed under the MIT license. Modified: 10 November 2010 Author: Original FORTRAN77 version by David Bailey. C version by John Burkardt. */ { int j; int jx; int k; double rld; double rld1; double rld2; double rldi; /* Start forward generation process and sweep. */ j = jl; for ( k = kl; k <= ku; k++ ) { rld = c[j+k*nja]; rldi = 1.0 / rld; f[j+k*nja+0*nja*njb] = f[j+k*nja+0*nja*njb] * rldi; f[j+k*nja+1*nja*njb] = f[j+k*nja+1*nja*njb] * rldi; f[j+k*nja+2*nja*njb] = f[j+k*nja+2*nja*njb] * rldi; x[j+k*nja] = d[j+k*nja] * rldi; y[j+k*nja] = e[j+k*nja] * rldi; } j = jl + 1; for ( k = kl; k <= ku; k++ ) { rld1 = b[j+k*nja]; rld = c[j+k*nja] - rld1 * x[j-1+k*nja]; rldi = 1.0 / rld; f[j+k*nja+0*nja*njb] = ( f[j+k*nja+0*nja*njb] - rld1 * f[j-1+k*nja+0*nja*njb] ) * rldi; f[j+k*nja+1*nja*njb] = ( f[j+k*nja+1*nja*njb] - rld1 * f[j-1+k*nja+1*nja*njb] ) * rldi; f[j+k*nja+2*nja*njb] = ( f[j+k*nja+2*nja*njb] - rld1 * f[j-1+k*nja+2*nja*njb] ) * rldi; x[j+k*nja] = ( d[j+k*nja] - rld1 * y[j-1+k*nja] ) * rldi; y[j+k*nja] = e[j+k*nja] * rldi; } for ( j = jl + 2; j <= ju - 2; j++ ) { for ( k = kl; k <= ku; k++ ) { rld2 = a[j+k*nja]; rld1 = b[j+k*nja] - rld2 * x[j-2+k*nja]; rld = c[j+k*nja] - ( rld2 * y[j-2+k*nja] + rld1 * x[j-1+k*nja] ); rldi = 1.0 / rld; f[j+k*nja+0*nja*njb] = ( f[j+k*nja+0*nja*njb] - rld2 * f[j-2+k*nja+0*nja*njb] - rld1 * f[j-1+k*nja+0*nja*njb] ) * rldi; f[j+k*nja+1*nja*njb] = ( f[j+k*nja+1*nja*njb] - rld2 * f[j-2+k*nja+1*nja*njb] - rld1 * f[j-1+k*nja+1*nja*njb] ) * rldi; f[j+k*nja+2*nja*njb] = ( f[j+k*nja+2*nja*njb] - rld2 * f[j-2+k*nja+2*nja*njb] - rld1 * f[j-1+k*nja+2*nja*njb] ) * rldi; x[j+k*nja] = ( d[j+k*nja] - rld1 * y[j-1+k*nja] ) * rldi; y[j+k*nja] = e[j+k*nja] * rldi; } } j = ju - 1; for ( k = kl; k <= ku; k++ ) { rld2 = a[j+k*nja]; rld1 = b[j+k*nja] - rld2 * x[j-2+k*nja]; rld = c[j+k*nja] - ( rld2 * y[j-2+k*nja] + rld1 * x[j-1+k*nja] ); rldi = 1.0 / rld;; f[j+k*nja+0*nja*njb] = ( f[j+k*nja+0*nja*njb] - rld2 * f[j-2+k*nja+0*nja*njb] - rld1 * f[j-1+k*nja+0*nja*njb] ) * rldi; f[j+k*nja+1*nja*njb] = ( f[j+k*nja+1*nja*njb] - rld2 * f[j-2+k*nja+1*nja*njb] - rld1 * f[j-1+k*nja+1*nja*njb] ) * rldi; f[j+k*nja+2*nja*njb] = ( f[j+k*nja+2*nja*njb] - rld2 * f[j-2+k*nja+2*nja*njb] - rld1 * f[j-1+k*nja+2*nja*njb] ) * rldi; x[j+k*nja] = ( d[j+k*nja] - rld1 * y[j-1+k*nja] ) * rldi; } j = ju; for ( k = kl; k <= ku; k++ ) { rld2 = a[j+k*nja]; rld1 = b[j+k*nja] - rld2 * x[j-2+k*nja]; rld = c[j+k*nja] - ( rld2 * y[j-2+k*nja] + rld1 * x[j-1+k*nja] ); rldi = 1.0 / rld; f[j+k*nja+0*nja*njb] = ( f[j+k*nja+0*nja*njb] - rld2 * f[j-2+k*nja+0*nja*njb] - rld1 * f[j-1+k*nja+0*nja*njb] ) * rldi; f[j+k*nja+1*nja*njb] = ( f[j+k*nja+1*nja*njb] - rld2 * f[j-2+k*nja+1*nja*njb] - rld1 * f[j-1+k*nja+1*nja*njb] ) * rldi; f[j+k*nja+2*nja*njb] = ( f[j+k*nja+2*nja*njb] - rld2 * f[j-2+k*nja+2*nja*njb] - rld1 * f[j-1+k*nja+2*nja*njb] ) * rldi; } /* Back sweep solution. */ for ( k = kl; k <= ku; k++ ) { f[ju+k*nja+0*nja*njb] = f[ju+k*nja+0*nja*njb]; f[ju+k*nja+1*nja*njb] = f[ju+k*nja+1*nja*njb]; f[ju+k*nja+2*nja*njb] = f[ju+k*nja+2*nja*njb]; f[ju-1+k*nja+0*nja*njb] = f[ju-1+k*nja+0*nja*njb] - x[ju-1+k*nja] * f[ju+k*nja+0*nja*njb]; f[ju-1+k*nja+1*nja*njb] = f[ju-1+k*nja+1*nja*njb] - x[ju-1+k*nja] * f[ju+k*nja+1*nja*njb]; f[ju-1+k*nja+2*nja*njb] = f[ju-1+k*nja+2*nja*njb] - x[ju-1+k*nja] * f[ju+k*nja+2*nja*njb]; } for ( j = 2; j <= ju - jl; j++ ) { jx = ju - j; for ( k = kl; k <= ku; k++ ) { f[jx+k*nja+0*nja*njb] = f[jx+k*nja+0*nja*njb] - x[jx+k*nja] * f[jx+1+k*nja+0*nja*njb] - y[jx+k*nja] * f[jx+2+k*nja+0*nja*njb]; f[jx+k*nja+1*nja*njb] = f[jx+k*nja+1*nja*njb] - x[jx+k*nja] * f[jx+1+k*nja+1*nja*njb] - y[jx+k*nja] * f[jx+2+k*nja+1*nja*njb]; f[jx+k*nja+2*nja*njb] = f[jx+k*nja+2*nja*njb] - x[jx+k*nja] * f[jx+1+k*nja+2*nja*njb] - y[jx+k*nja] * f[jx+2+k*nja+2*nja*njb]; } } return; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY copies an R8VEC. Discussion: An R8VEC is a vector of R8"s. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Input, double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* 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 MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ double wtime ( ) /******************************************************************************/ /* Purpose: WTIME reports the elapsed wallclock time. Discussion: The reliability of this function depends in part on the value of CLOCKS_PER_SECOND. Licensing: This code is distributed under the MIT license. Modified: 27 April 2009 Author: John Burkardt Parameters: Output, double WTIME, the a reading of the wall clock timer, in seconds. */ { double value; value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC; return value; }