# include <cstdlib>
# include <iomanip>
# include <iostream>
# include <string>
# include <cmath>
# include <ctime>
# include <complex>

using namespace std;

int main ( );
void mxmtst ( double *er, double *fp, double *tm );
void mxm ( double a[], double b[], double c[], int l, int m, int n );
void ffttst ( double *er, double *fp, double *tm );
void cfft2d1 ( int is, int m, int m1, int n, complex <double> x[], 
  complex <double> w[], int ip[] );
void cfft2d2 ( int is, int m, int m1, int n, complex <double> x[], 
  complex <double> w[], int ip[] );
void chotst ( 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 btrtst ( 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 gmttst ( double *er, double *fp, double *tm );
void gmtry ( int nb, int nw, int nwall[], complex <double> proj[], double rmatrx[], 
  complex <double> wall[], double xmax[], complex <double> zcr[] );
void emitst ( double *er, double *fp, double *tm );
void emit ( int nb, int nvm, int nw, double cp[], double dpds[], 
  complex <double> expmz[], complex <double> expz[], complex <double> force[], 
  double gamma[], int nwall[], double ps[], double psi[], 
  complex <double> refpt[], double rhs[], double rmatrx[], double rmom[], 
  complex <double> wall[], complex <double> z[], complex <double> zcr[] );
void vpetst ( 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 );
double r8_abs ( double x );
double r8_max ( double x, double y );
double r8_min ( double x, double y );
void r8vec_copy ( int n, double a1[], double a2[] );
void timestamp ( );
double wtime ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for NAS.
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    Original FORTRAN77 version by David Bailey.
//    C++ version by John Burkardt.
//
{
  double er;
  double er_total;
  double fp;
  double fp_total;
  int i;
  string pn;
  double rt;
  double tm;
  double tm_total;

  er_total = 0.0;
  fp_total = 0.0;
  tm_total = 0.0;

  timestamp ( );
  cout << "\n";
  cout << "NAS:\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "                The NAS kernel benchmark program\n";
  cout << "\n";
  cout << " Program        Error          FP Ops";
  cout << "        Seconds     MFLOPS\n";
  cout << "\n";

  for ( i = 1; i <= 7; i++ )
  {
    if ( i == 1 )
    {
      pn = "MXM     ";
      mxmtst ( &er, &fp, &tm );
    }
    else if ( i == 2 )
    {
      pn = "CFFT2D  ";
      ffttst ( &er, &fp, &tm );
    }
    else if ( i == 3 )
    {
      pn = "CHOLSKY ";
      chotst ( &er, &fp, &tm );
    }
    else if ( i == 4 )
    {
      pn = "BTRIX   ";
      btrtst ( &er, &fp, &tm );
    }
    else if ( i == 5 )
    {
      pn = "GMTRY   ";
      gmttst ( &er, &fp, &tm );
    }
    else if ( i == 6 )
    {
      pn = "EMIT    ";
      emitst ( &er, &fp, &tm );
    }
    else if ( i == 7 )
    {
      pn = "VPENTA  ";
      vpetst ( &er, &fp, &tm );
    }
    rt = 1.0E-06 * fp / tm;
    cout << " "  << setw(8)  << pn
         << "  " << setw(13) << er
         << "  " << setw(13) << fp
         << "  " << setw(10) << tm
         << "  " << setw(10) << rt << "\n";

    er_total = er_total + er;
    fp_total = fp_total + fp;
    tm_total = tm_total + tm;
  }
  pn = "Total   ";
  rt = 1.0E-06 * fp_total / tm_total;
  cout << "\n";
  cout << " "  << setw(8)  << pn
       << "  " << setw(13) << er_total
       << "  " << setw(13) << fp_total
       << "  " << setw(10) << tm_total
       << "  " << setw(10) << rt << "\n";

  cout << "\n";
  cout << "NAS:\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void mxmtst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    MXMTST is the test program for MXM.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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 = r8_abs ( ( c[18+18*l] - ans ) / ans );
  *fp = ( double ) ( 2 * it * l * m * n );

  return;
# undef L
# undef M
# undef N
}
//****************************************************************************80

void mxm ( double a[], double b[], double c[], int l, int m, int n )

//****************************************************************************80
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 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;
}
//****************************************************************************80

void ffttst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    FFTTST is the test program for CFFT2D1 and CFFTD2.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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;
  complex <double> ct;
  double f7;
  int i;
  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;
  complex <double> w1[M];
  complex <double> w2[N];
  complex <double> x[M1*N];
  complex <double> y;

  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] = complex <double> ( t1, t2 );
    }
  }
  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 ( ( real ( 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
}
//****************************************************************************80

void cfft2d1 ( int is, int m, int m1, int n, complex <double> x[], 
  complex <double> w[], int ip[] )

//****************************************************************************80
//
//  Purpose:
//
//    CFFT2D1 performs complex radix 2 FFT''s on the first dimension.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    Original FORTRAN77 version by David Bailey.
//    C++ version by John Burkardt.
//
{
  complex <double> ct;
  complex <double> 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] = complex <double> ( cos ( t ), sin ( t ) );
    }
    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;
}
//****************************************************************************80

void cfft2d2 ( int is, int m, int m1, int n, complex <double> x[], 
  complex <double> w[], int ip[] )

//****************************************************************************80
//
//  Purpose:
//
//    CFFT2D2 performs complex radix 2 FFT''s on the second dimension.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    Original FORTRAN77 version by David Bailey.
//    C++ version by John Burkardt.
//
{
  complex <double> ct;
  complex <double> 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] = complex <double> ( cos ( t ), sin ( t ) );
    }
    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;
}
//****************************************************************************80

void chotst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    CHOTST is the test program for CHOLSKY.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  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;

  for ( k = 0; k <= n; k++ )
  {
    for ( j = -m; j <= 0; j++ )
    {
      for ( i = 0; i <= ida; i++ )
      {
        t = fmod ( f7 * t, 1.0 );
        i1 = i + ( j + m ) * ( ida + 1 ) + k * ( ida + 1 ) * ( m + 1 );
        ax[i1] = t;
      }
    }
  }

  for ( k = 0; k <= n; k++ )
  {
    for ( j = 0; j <= nmat; j++ )
    {
      for ( i = 0; i <= nrhs; i++ )
      {
        t = fmod ( f7 * t, 1.0 );
        i1 = i+j*(nrhs+1)+k*(nrhs+1)*(nmat+1);
        bx[i1] = t;
      }
    }
  }
//
//  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 = r8_abs ( ( b[i1] - ans ) / ans );
  *fp = ( double ) ( it * ( nmat + 1 ) * 4403 );

  return;
# undef IDA
# undef M
# undef N
# undef NMAT
# undef NRHS
}
//****************************************************************************80

void cholsky ( int ida, int nmat, int m, int n, double a[], int nrhs, int idb, 
  double b[] )

//****************************************************************************80
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 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 ( r8_abs ( 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;
}
//****************************************************************************80

void btrtst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    BTRTST is the test program for BTRIX.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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 = r8_abs ( ( 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
}
//****************************************************************************80

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[] )

//****************************************************************************80
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 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
}
//****************************************************************************80

void gmttst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    GMTTST is the test program for GMTRY.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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 lw;
  int nb = NB;
  int nw = NW;
  int nwall[NB];
  complex <double> proj[NW*NB];
  double rmatrx[NW*NB*NW*NB];
  double t1;
  double t2;
  double t30;
  double time1;
  complex <double> wall[NW*NB];
  double xmax[NB];
  complex <double> z1;
  complex <double> zcr[NW*NB];
  complex <double> zi;
  complex <double> zz;

  it = 2;
  ans = -2.57754233214174;
  lw = 2 * nw * nb;
//
//  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] = complex <double> ( t1, t2 );
    }
  }
//
//  Timing.
//
  time1 = wtime ( );

  for ( i = 1; i <= it; i++ )
  {
    gmtry ( nb, nw, nwall, proj, rmatrx, wall, xmax, zcr );
  }

  *tm = wtime ( ) - time1;
//
//  Results.
//
  *er = r8_abs ( ( 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
}
//****************************************************************************80

void gmtry ( int nb, int nw, int nwall[], complex <double> proj[], double rmatrx[], 
  complex <double> wall[], double xmax[], complex <double> zcr[] )

//****************************************************************************80
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    Original FORTRAN77 version by David Bailey.
//    C++ version by John Burkardt.
//
{
  double arcl;
  double dum;
  int i;
  complex <double> I = complex <double> ( 0, 1 );
  int i0;
  int j;
  int j0;
  int k;
  int k1;
  int k2;
  int kp;
  int kron;
  int ks;
  int l;
  int l1;
  int l2;
  int ls;
  int matdim;
  double period;
  double pi;
  double pidp;
  double r0;
  double sig2;
  double sigma;
  double ylimit;
  double ymax;
  double ymin;
  complex <double> z1;
  complex <double> zi;
  complex <double> 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 + abs ( 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 / abs ( 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] = real ( zcr[0+(l-1)*nw] );
    ls = 0;

    for ( k = 1; k <= nwall[l-1]; k++ )
    {
      ymin = r8_min ( ymin, imag ( zcr[k-1+(l-1)*nw] ) );
      ymax = r8_max ( ymax, imag ( zcr[k-1+(l-1)*nw] ) );
      xmax[l-1] = r8_max ( xmax[l-1], real ( zcr[k-1+(l-1)*nw] ) );
      kp = 1 + ( k % nwall[l-1] );

      if ( 0.0 < real ( ( zcr[kp-1+(l-1)*nw] - zcr[k-1+(l-1)*nw] ) * 
        conj ( wall[kp-1+(l-1)*nw] - wall[k-1+(l-1)*nw] ) ) )
      {
        ls = l;
        ks = k;
      }
    }
  }
//
//  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 + ( imag ( 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 = exp ( ( wall[0+(l1-1)*nw] - zcr[j-1+(l2-1)*nw] ) * pidp );
        z1 = z1 - 1.0 / z1;
        dum = sig2 + pow ( real ( z1 ), 2 ) + pow ( imag ( z1 ), 2 );
        for ( i = 2; i <= nwall[l1-1]; i++ )
        {
          zi = exp ( ( 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 ( real ( zz ), 2 ) + pow ( imag ( 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;
}
//****************************************************************************80

void emitst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    EMITST is the test program for EMIT.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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];
  complex <double> dum3;
  complex <double> expmwk;
  complex <double> expmz[NVM];
  complex <double> expwkl;
  complex <double> expz[NVM];
  double f7;
  complex <double> 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];
  complex <double> refpt[NB];
  double rhs[NW*NB];
  double rmatrx[NW*NB*NW*NB];
  double rmom[NB];
  double t1;
  double t2;
  double t30;
  double time1;
  complex <double> uupstr;
  complex <double> wall[NW*NB];
  complex <double> z[NVM];
  complex <double> 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] = complex <double> ( t1, t2 );
      t1 = fmod ( f7 * t2, 1.0 );
      t2 = fmod ( f7 * t1, 1.0 );
      zcr[i+j*nw] = complex <double> ( t1, t2 );
      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] = complex <double> ( t1, t2 );
    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 = r8_abs ( ( 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
}
//****************************************************************************80

void emit ( int nb, int nvm, int nw, double cp[], double dpds[], 
  complex <double> expmz[], complex <double> expz[], complex <double> force[], 
  double gamma[], int nwall[], double ps[], double psi[], 
  complex <double> refpt[], double rhs[], double rmatrx[], double rmom[], 
  complex <double> wall[], complex <double> z[], complex <double> zcr[] )

//****************************************************************************80
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    Original FORTRAN77 version by David Bailey.
//    C++ version by John Burkardt.
//
{
  double chord;
  double cpm;
  double cupst;
  double delt;
  complex <double> dum3;
  complex <double> expmwk;
  complex <double> expwkl;
  int i;
  complex <double> I = complex <double> ( 0, 1 );
  int i0;
  int j;
  int k;
  int k1;
  int k2;
  int l;
  int matdim;
  int nolld;
  int nv;
  double period;
  double pi;
  double pidp;
  double sig2;
  double sps;
  double u0;
  complex <double> uupstr;
  complex <double> zz;

  period = 3.0;
  sig2 = 3.0;
  u0 = 4.0;
  matdim = 500;
  delt = 1.0;
  chord = 5.0;
  pi = 3.141592653589793;
  uupstr = complex <double> ( 3.0, 4.0 );
//
//  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] = exp ( z[i] * pidp );
    expmz[i] = 1.0 / expz[i];
  }
  
  i0 = 0;
  cupst = pow ( real ( uupstr ), 2 ) + pow ( imag ( uupstr ), 2 );

  for ( l = 0; l < nb; l++ )
  {
    for ( k = 0; k < nwall[l]; k++ )
    {
      expwkl = exp ( 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 ( real ( dum3 ), 2 ) + 
          pow ( imag ( dum3 ), 2 ) + sig2 );
        sps = sps + ps[i];
      }
      psi[k] = imag ( 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.
//
  nolld = nv;
  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 ( real ( z[nv] - refpt[l] ), 2 ) +
        pow ( imag ( 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 = r8_max ( 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;
}
//****************************************************************************80

void vpetst ( double *er, double *fp, double *tm )

//****************************************************************************80
//
//  Purpose:
//
//    VPETST is the test program for VPENTA.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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 = r8_abs ( ( 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
}
//****************************************************************************80

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[] )

//****************************************************************************80
//
//  Purpose:
//
//    VPENTA inverts 3 pentadiagonal systems simultaneously.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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;
}
//****************************************************************************80

int i4_max ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MAX returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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;
}
//****************************************************************************80

int i4_min ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MIN returns the smaller of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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;
}
//****************************************************************************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:
//
//    07 May 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

double r8_max ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MAX returns the maximum of two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    07 May 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MAX, the maximum of X and Y.
//
{
  double value;

  if ( y < x )
  {
    value = x;
  } 
  else
  {
    value = y;
  }
  return value;
}
//****************************************************************************80

double r8_min ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MIN returns the minimum of two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    07 May 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MIN, the minimum of X and Y.
//
{
  double value;

  if ( y < x )
  {
    value = y;
  } 
  else
  {
    value = x;
  }
  return value;
}
//****************************************************************************80

void r8vec_copy ( int n, double a1[], double a2[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_COPY copies an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8"s.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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;
}
//****************************************************************************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:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}
//****************************************************************************80

double wtime ( )

//****************************************************************************80
//
//  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 GNU LGPL 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;
}