function nas ( )

%*****************************************************************************80
%
%% 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:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  te = 0.0;
  tf = 0.0;
  tt = 0.0;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NAS:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                The NAS kernel benchmark program\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ' Program        Error          FP Ops' );
  fprintf ( 1, '        Seconds     MFLOPS\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 7

    if ( i == 1 )
      pn = 'MXM     ';
      [ er(1), fp(1), tm(1) ] = mxmtst ( );
    elseif ( i == 2 )
      pn = 'CFFT2D  ';
      [ er(2), fp(2), tm(2) ] = ffttst ( );
    elseif ( i == 3 )
      pn = 'CHOLSKY ';
      [ er(3), fp(3), tm(3) ] = chotst ( );
    elseif ( i == 4 )
      pn = 'BTRIX   ';
      [ er(4), fp(4), tm(4) ] = btrtst ( );
    elseif ( i == 5 )
      pn = 'GMTRY   ';
      [ er(5), fp(5), tm(5) ] = gmttst ( );
    elseif ( i == 6 )
      pn = 'EMIT    ';
      [ er(6), fp(6), tm(6) ] = emitst ( );
    elseif ( i == 7 )
      pn = 'VPENTA  ';
      [ er(7), fp(7), tm(7) ] = vpetst ( );
    end

    rt(i) = 1.0D-06 * fp(i) / tm(i);
    fprintf ( 1, ' %s  %13.4e  %13.4e  %10.4e  %10.2e\n', ...
      pn, er(i), fp(i), tm(i), rt(i) );

    te = te + er(i);
    tf = tf + fp(i);
    tt = tt + tm(i);

  end

  er(8) = te;
  fp(8) = tf;
  tm(8) = tt;
  rt(8) = 1.0E-06 * tf / tt;
  i = 8;
  pn = 'Total   ';
  fprintf ( 1, '\n' );
  fprintf ( 1, ' %s  %13.4e  %13.4e  %10.4e  %10.2e\n', ...
    pn, er(i), fp(i), tm(i), rt(i) );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NAS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ er, fp, tm ] = mxmtst ( )

%*****************************************************************************80
%
%% MXMTST is the test program for MXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  l = 256;
  m = 128;
  n = 64;

  it = 100;
  ans = 35.2026179738722;
%
%  Random initialization.
%
  f7 = 78125.0;
  t30 = 1073741824.0;
  t = f7 / t30;

  a = zeros(l,m);
  for j = 1 : m
    for i = 1 : l
      t = mod ( f7 * t, 1.0 );
      a(i,j) = t;
    end
  end

  b = zeros(m,n);
  for j = 1 : n
    for i = 1 : m
      t = mod ( f7 * t, 1.0 );
      b(i,j) = t;
    end
  end
%
%  Timing.
%
  tic;

  for ii = 1 : it
    c = mxm ( a, b, l, m, n );
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( c(19,19) - ans ) / ans );
  fp = 2 * it * l * m * n;

  return
end
function c = mxm ( a, b, l, m, n )

%*****************************************************************************80
%
%% 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:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  c = a * b;

  return
end
function [ er, fp, tm ] = ffttst ( )

%*****************************************************************************80
%
%% FFTTST is the test program for CFFT2D1 and CFFTD2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  m = 128;
  m1 = 128;
  n = 256;

  I = sqrt ( - 1 );
  it = 100;
  ans = 0.894799941219277;
  rmn = 1.0 / ( m * n );
%
%  Random initialization.
%
  f7 = 78125.0;
  t30 = 1073741824.0;
  t2 = f7 / t30;

  for j = 1 : n
    for i = 1 : m
      t1 = mod ( f7 * t2, 1.0 );
      t2 = mod ( f7 * t1, 1.0 );
      x(i,j) = t1 + I * t2;
    end
  end

  w1 = [];
  [ x, w1 ] = cfft2d1 ( 0, m, m1, n, x, w1 );
  w2 = [];
  [ x, w2 ] = cfft2d2 ( 0, m, m1, n, x, w2 );
%
%  Timing.
%
  tic;

  for k = 1 : it

    x(1:m,1:n) = rmn * x(1:m,1:n);

    [ x, w1 ] = cfft2d1 ( 1, m, m1, n, x, w1 );
    [ x, w2 ] = cfft2d2 ( 1, m, m1, n, x, w2 );
    [ x, w2 ] = cfft2d2 ( -1, m, m1, n, x, w2 );
    [ x, w1 ] = cfft2d1 ( -1, m, m1, n, x, w1 );

  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( ( x(19,19) ) - ans ) / ans );
  fp = it * m * n * ( 2.0 + 10.0 * log ( m * n ) / log ( 2.0 ) );

  return
end
function [ x, w ] = cfft2d1 ( is, m, m1, n, x, w )

%*****************************************************************************80
%
%% CFFT2D1 performs complex radix 2 FFT's on the first dimension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  I = sqrt ( - 1 );
%
%  If IS = 0 then initialize only.
%
  m2 = floor ( m / 2 );

  if ( is == 0 )
    w = zeros ( m2, 1 );
    for i = 1 : m2
      t = 2.0 * pi * ( i - 1 ) / m;
      w(i) = cos ( t ) + I * sin ( t );
    end
    return
  end
%
%  Perform forward or backward FFT's according to IS = 1 or -1.
%
  ip = zeros(2,m);

  for i = 1 : m
    ip(1,i) = i;
  end

  l = 1;
  i1 = 1;

  while ( 1 )

    i2 = 3 - i1;

    for j = l : l : m2

      cx = w(j-l+1);
      if ( is < 0 )
        cx = conj ( cx );
      end

      for i = j - l + 1 : j
        ii = ip(i1,i);
        ip(i2,i+j-l) = ii;
        im = ip(i1,i+m2);
        ip(i2,i+j) = im;
        for k = 1 : n
          ct = x(ii,k) - x(im,k);
          x(ii,k) = x(ii,k) + x(im,k);
          x(im,k) = ct * cx;
        end
      end

    end

    l = 2 * l;
    i1 = i2;

    if ( m2 < l )
      break
    end

  end

  for i = 1 : m
    ii = ip(i1,i);
    if ( i < ii )
      for k = 1 : n
        ct = x(i,k);
        x(i,k) = x(ii,k);
        x(ii,k) = ct;
      end
    end
  end

  return
end
function [ x, w ] = cfft2d2 ( is, m, m1, n, x, w )

%*****************************************************************************80
%
%% CFFT2D2 performs complex radix 2 FFT's on the second dimension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  I = sqrt ( - 1 );
%
%  If IS = 0, then initialize only.
%
  n2 = floor ( n / 2 );

  if ( is == 0 )
    w = zeros ( n2, 1 );
    for i = 1 : n2
      t = 2.0 * pi * ( i - 1 ) / n;
      w(i) = cos ( t ) + I * sin ( t );
    end
    return
  end
%
%  Perform forward or backward FFT's according to IS = 1 or -1.
%
  ip = zeros(2,n);

  for i = 1 : n
    ip(1,i) = i;
  end

  l = 1;
  i1 = 1;

  while ( 1 )

    i2 = 3 - i1;

    for j = l : l : n2

      cx = w(j-l+1);
      if ( is < 0 )
        cx = conj ( cx );
      end

      for i = j - l + 1 : j
        ii = ip(i1,i);
        ip(i2,i+j-l) = ii;
        im = ip(i1,i+n2);
        ip(i2,i+j) = im;
        for k = 1 : m
          ct = x(k,ii) - x(k,im);
          x(k,ii) = x(k,ii) + x(k,im);
          x(k,im) = ct * cx;
        end
      end
    end

    l = 2 * l;
    i1 = i2;

    if ( n2 < l )
      break
    end

  end

  for i = 1 : n
    ii = ip(i1,i);
    if ( i < ii )
      for k = 1 : m
        ct = x(k,i);
        x(k,i) = x(k,ii);
        x(k,ii) = ct;
      end
    end
  end

  return
end
function [ er, fp, tm ] = chotst ( )

%*****************************************************************************80
%
%% CHOTST is the test program for CHOLSKY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  ida = 250;
  m = 4;
  n = 40;
  nmat = 250;
  nrhs = 3;

  it = 200;
  it = 1;
  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;

  ax = zeros(ida+1,m+1,n+1);

  for k = 0 : n
    for j = - m : 0
      for i = 0 : ida
        t = mod ( f7 * t, 1.0 );
        ax(i+1,j+m+1,k+1) = t;
      end
    end
  end

  bx = zeros(nrhs+1,nmat+1,n+1);

  for k = 0 : n
    for j = 0 : nmat
      for i = 0 : nrhs
        t = mod ( f7 * t, 1.0 );
        bx(i+1,j+1,k+1) = t;
      end
    end
  end
%
%  Timing.
%
  
  tic;

  for j = 1 : it
    a = ax;
    b = bx;
    b = cholsky ( ida, nmat, m, n, a, nrhs, ida, b );
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( b(1+1,19+1,19+1) - ans ) / ans );
  fp = it * ( nmat + 1 ) * 4403;

  return
end
function b = cholsky ( ida, nmat, m, n, a, nrhs, idb, b )

%*****************************************************************************80
%
%% 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:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  eps = 1.0E-13;
%
%  Cholesky decomposition.
%
  for j = 0 : n

    i0 = max ( -m, -j );
%
%  Off diagonal elements.
%
    for i = i0 : -1
      for jj = i0 - i : -1
        for l = 0 : nmat
          a(l+1,i+m+1,j+1) = a(l+1,i+m+1,j+1) ...
            - a(l+1,jj+m+1,i+j+1) * a(l+1,i+jj+m+1,j+1);
        end
      end
      for l = 0 : nmat
        a(l+1,i+m+1,j+1) = a(l+1,i+m+1,j+1) * a(l+1,0+m+1,i+j+1);
      end
    end
%
%  Store inverse of diagonal elements.
%
    for l = 0 : nmat
      epss(l+1) = eps * a(l+1,0+m+1,j+1);
    end

    for jj = i0 : -1
      for l = 0 : nmat
        a(l+1,0+m+1,j+1) = a(l+1,0+m+1,j+1) - a(l+1,jj+m+1,j+1).^2;
      end
    end

    for l = 0 : nmat
      a(l+1,0+m+1,j+1) = 1.0 / sqrt ( abs ( epss(l+1) + a(l+1,0+m+1,j+1) ) );
    end

  end
%
%  Solution.
%
  for i = 0 : nrhs

    for k = 0 : n

      for l = 0 : nmat
        b(i+1,l+1,k+1) = b(i+1,l+1,k+1) * a(l+1,0+m+1,k+1);
      end

      for jj = 1 : min ( m, n - k )
        for l = 0 : nmat
          b(i+1,l+1,k+jj+1) = b(i+1,l+1,k+jj+1) ...
            - a(l+1,-jj+m+1,k+jj+1) * b(i+1,l+1,k+1);
        end
      end

    end

    for k = n : -1 : 0

      for l = 0 : nmat
        b(i+1,l+1,k+1) = b(i+1,l+1,k+1) * a(l+1,0+m+1,k+1);
      end

      for jj = 1 : min ( m, k )
        for l = 0 : nmat
          b(i+1,l+1,k-jj+1) = b(i+1,l+1,k-jj+1) ...
            - a(l+1,-jj+m+1,k+1) * b(i+1,l+1,k+1);
        end
      end

    end

  end

  return
end
function [ er, fp, tm ] = btrtst ( )

%*****************************************************************************80
%
%% BTRTST is the test program for BTRIX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  jd = 30;
  kd = 30;
  ld = 30;
  md = 30;

  js = 2;
  je = 29;
  ls = 2;
  le = 29;
  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;

  a = zeros(5,5,md,md);
  b = zeros(5,5,md,md);
  c = zeros(5,5,md,md);

  for l = 1 : md
    for k = 1 : md
      for j = 1 : 5
        for i = 1 : 5
          t = mod ( f7 * t, 1.0 );
          a(i,j,k,l) = t;
          t = mod ( f7 * t, 1.0 );
          bx(i,j,k,l) = t;
          t = mod ( f7 * t, 1.0 );
          c(i,j,k,l) = t;
        end
      end
    end
  end

  sx = zeros(jd,kd,ld,5);

  for l = 1 : 5
    for k = 1 : ld
      for j = 1 : kd
        for i = 1 : jd
          t = mod ( f7 * t, 1.0 );
          sx(i,j,k,l) = t;
        end
      end
    end
  end
%
%  Timing.
%
  tic;

  for ii = 1 : it
    s(1:jd,1:kd,1:ld,1:5) = sx(1:jd,1:kd,1:ld,1:5);
    for k = 1 : kd
      b(1:5,1:5,1:md,1:md) = bx(1:5,1:5,1:md,1:md);
      s = btrix ( js, je, ls, le, k, jd, kd, ld, md, a, b, c, s );
    end
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( s(19,19,19,1) - ans ) / ans );
  fp = it * md * ( le - 1 ) * 19165;

  return
end
function s = btrix ( js, je, ls, le, k, jd, kd, ld, md, a, b, c, s )

%*****************************************************************************80
%
%% 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:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%

%
%  Part 1.  Forward block sweep.
%
  for j = js : je
%
%  Step 1.  Construct L(I) in B.
%
    if ( j ~= js )
      for m = 1 : 5
        for n = 1 : 5
          for l = ls : le
            b(m,n,j,l) = b(m,n,j,l)       ...
              - a(m,1,j,l) * b(1,n,j-1,l) ...
              - a(m,2,j,l) * b(2,n,j-1,l) ...
              - a(m,3,j,l) * b(3,n,j-1,l) ...
              - a(m,4,j,l) * b(4,n,j-1,l) ...
              - a(m,5,j,l) * b(5,n,j-1,l);
          end
        end
      end
    end
%
%  Step 2.  Compute L inverse.
%
%  A.  Decompose L(I) into L and U.
%
    for l = ls : le
      l11(l) = 1.0 / b(1,1,j,l);
      u12(l) = b(1,2,j,l) * l11(l);
      u13(l) = b(1,3,j,l) * l11(l);
      u14(l) = b(1,4,j,l) * l11(l);
      u15(l) = b(1,5,j,l) * l11(l);
      l21(l) = b(2,1,j,l);
      l22(l) = 1.0 / ( b(2,2,j,l) - l21(l) * u12(l) );
      u23(l) = ( b(2,3,j,l) - l21(l) * u13(l) ) * l22(l);
      u24(l) = ( b(2,4,j,l) - l21(l) * u14(l) ) * l22(l);
      u25(l) = ( b(2,5,j,l) - l21(l) * u15(l) ) * l22(l);
      l31(l) = b(3,1,j,l);
      l32(l) = b(3,2,j,l) - l31(l) * u12(l);
      l33(l) = 1.0 / ( b(3,3,j,l) - l31(l) * u13(l) - l32(l) * u23(l) );
      u34(l)  = ( b(3,4,j,l) - l31(l) * u14(l) - l32(l) * u24(l) ) * l33(l);
      u35(l)  = ( b(3,5,j,l) - l31(l) * u15(l) - l32(l) * u25(l) ) * l33(l);
    end

    for l = ls : le
      l41(l) = b(4,1,j,l);
      l42(l) = b(4,2,j,l) - l41(l) * u12(l);
      l43(l) = b(4,3,j,l) - l41(l) * u13(l) - l42(l) * u23(l);
      l44(l) = 1.0 / ( b(4,4,j,l) - l41(l) * u14(l) ...
        - l42(l) * u24(l) - l43(l) * u34(l) );
      u45(l) = ( b(4,5,j,l) - l41(l) * u15(l) - l42(l) * u25(l) ...
        - l43(l) * u35(l) ) * l44(l);
      l51(l) = b(5,1,j,l);
      l52(l) = b(5,2,j,l) - l51(l) * u12(l);
      l53(l) = b(5,3,j,l) - l51(l) * u13(l) - l52(l) * u23(l);
      l54(l) = b(5,4,j,l) - l51(l) * u14(l) - l52(l) * u24(l) - l53(l) * u34(l);
      l55(l) = 1.0 / ( b(5,5,j,l) - l51(l) * u15(l) - l52(l) * u25(l) ...
        - l53(l) * u35(l) - l54(l) * u45(l) );
    end
%
%  Step 3.  Solve for intermediate vector.
%
%  A.  Construct the right hand side.
%
    if ( j ~= js )
      for m = 1 : 5
        for l = ls : le
          s(j,k,l,m) = s(j,k,l,m)       ...
            - a(m,1,j,l) * s(j-1,k,l,1) ...
            - a(m,2,j,l) * s(j-1,k,l,2) ...
            - a(m,3,j,l) * s(j-1,k,l,3) ...
            - a(m,4,j,l) * s(j-1,k,l,4) ... 
            - a(m,5,j,l) * s(j-1,k,l,5);
        end
      end
    end
%
%  B. Intermediate vector.
%
%  Forward substitution.
%
    for l = ls : le
      d1 = s(j,k,l,1) * l11(l);
      d2 = ( s(j,k,l,2) - l21(l) * d1 ) * l22(l);
      d3 = ( s(j,k,l,3) - l31(l) * d1 - l32(l) * d2 ) * l33(l);
      d4 = ( s(j,k,l,4) - l41(l) * d1 - l42(l) * d2 - l43(l) * d3 ) * l44(l);
      d5 = ( s(j,k,l,5) - l51(l) * d1 - l52(l) * d2 - l53(l) * d3 ...
        - l54(l) * d4 ) * l55(l);
%
%  Backward substitution.
%
      s(j,k,l,5) = d5;
      s(j,k,l,4) = d4 - u45(l) * d5;
      s(j,k,l,3) = d3 - u34(l) * s(j,k,l,4) - u35(l) * d5;
      s(j,k,l,2) = d2 - u23(l) * s(j,k,l,3) - u24(l) * s(j,k,l,4) - u25(l) * d5;
      s(j,k,l,1) = d1 - u12(l) * s(j,k,l,2) - u13(l) * s(j,k,l,3) ...
        - u14(l) * s(j,k,l,4) - u15(l) * d5;
    end
%
%  Step 4.  Construct U(I) = inverse(L(I))*C(I+1) by columns and store in B.
%
    if ( j ~= je )

      for n = 1 : 5
        for l = ls : le
%
%  Forward substitution.
%
          c1 = c(1,n,j,l) * l11(l);
          c2 = ( c(2,n,j,l) - l21(l) * c1 ) * l22(l);
          c3 = ( c(3,n,j,l) - l31(l) * c1 - l32(l) * c2 ) * l33(l);
          c4 = ( c(4,n,j,l) - l41(l) * c1 - l42(l) * c2 - l43(l) * c3 ) * l44(l);
          c5 = ( c(5,n,j,l) - l51(l) * c1 - l52(l) * c2 ...
            - l53(l) * c3 - l54(l) * c4 ) * l55(l);
%
%  Backward substitution.
%
          b(5,n,j,l) = c5;
          b(4,n,j,l) = c4 - u45(l) * c5;
          b(3,n,j,l) = c3 - u34(l) * b(4,n,j,l) - u35(l) * c5;
          b(2,n,j,l) = c2 - u23(l) * b(3,n,j,l) - u24(l) * b(4,n,j,l) ...
            - u25(l) * c5;
          b(1,n,j,l) = c1 - u12(l) * b(2,n,j,l) - u13(l) * b(3,n,j,l) ...
            - u14(l) * b(4,n,j,l) - u15(l) * c5;
        end
      end

    end

  end
%
%  Part 2.  Backward block sweep.
%
  jem1 = je - 1;

  for j = jem1 : -1 : js
    for m = 1 : 5
      for l = ls : le
        s(j,k,l,m) = s(j,k,l,m)       ...
          - b(m,1,j,l) * s(j+1,k,l,1) ...
          - b(m,2,j,l) * s(j+1,k,l,2) ...
          - b(m,3,j,l) * s(j+1,k,l,3) ...
          - b(m,4,j,l) * s(j+1,k,l,4) ... 
          - b(m,5,j,l) * s(j+1,k,l,5);
      end
    end
  end

  return
end
function [ er, fp, tm ] = gmttst ( )

%*****************************************************************************80
%
%% GMTTST is the test program for GMTRY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  nb = 5;
  nw = 100;

  I = sqrt ( - 1 );

  it = 2;
  ans = -2.57754233214174;
  lw = 2 * nw * nb;
%
%  Random initialization.
%
  f7 = 78125.0;
  t30 = 1073741824.0;
  t2 = f7 / t30;

  nwall(1:nb) = nw;

  wall = zeros(nw,nb);
  for j = 1 : nb
    for i = 1 : nw
      t1 = mod ( f7 * t2, 1.0 );
      t2 = mod ( f7 * t1, 1.0 );
      wall(i,j) = t1 + t2 * I;
    end
  end
%
%  Timing.
%
  tic;

  for i = 1 : it
    rmatrx = gmtry ( nb, nw, nwall, wall );
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( rmatrx(19,19) - ans ) / ans );
  fp = it * ( 120 * ( nb * nw )^2 + 0.666 * ( nb * nw )^3 );

  return
end
function rmatrx = gmtry ( nb, nw, nwall, wall )

%*****************************************************************************80
%
%% 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:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  I = sqrt ( - 1 );
  period = 3.0;
%
%  Compute arclength.
%
  matdim = 0;
  arcl = 0.0;
  ymin = 1.0E+30;
  ymax = -1.0E+30;
  pidp = pi / period;

  for l = 1 : nb
    matdim = matdim + nwall(l);
    for k = 1 : nwall(l)
      arcl = arcl + abs ( wall(k,l) - wall(1+mod(k,nwall(l)),l) );
    end
  end
%
%  Compute core radius.
%
  r0 = 0.5 * arcl / matdim;
  sigma = r0 / 2.0;
%
%  Define creation points.
%
  xmax = zeros(nb,1);
  zcr = zeros(nw,nb);

  for l = 1 : nb

    for k = 1 : nwall(l)
      zz = wall(1+mod(k+nwall(l)-2,nwall(l)),l) - wall(1+mod(k,nwall(l)),l);
      zcr(k,l) = wall(k,l) + r0 * I / abs ( zz ) * zz;
    end
%
%  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) = real ( zcr(1,l) );
    ls = 0;

    for k = 1 : nwall(l)

      ymin = min ( ymin, imag ( zcr(k,l) ) );
      ymax = max ( ymax, imag ( zcr(k,l) ) );
      xmax(l) = max ( xmax(l), real ( zcr(k,l) ) );
      kp = 1 + mod ( k, nwall(l) );

      if ( 0.0 < real ( ( zcr(kp,l) - zcr(k,l) ) * ...
          conj ( wall(kp,l) - wall(k,l) ) ) )
        ls = l;
        ks = k;
      end

    end

  end
%
%  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.
%
  proj = zeros(nw,nb);

  for l = 1 : nb
    for k = 1 : nwall(l)
      proj(k,l) = zcr(k,l) - I * period * ...
        ( floor ( 5.0 + ( imag ( zcr(k,l) ) - ylimit ) / period ) ...
        - 5.0 );
    end
  end
%
%  Compute matrix.
%
  rmatrx = zeros(nw*nb,nw*nb);
  sig2 = ( 2.0 * pidp * sigma )^2;
  i0 = 0;

  for l1 = 1 : nb

    j0 = 0;

    for l2 = 1 : nb

      if ( l1 == l2 )
        kron = 1;
      else
        kron = 0;
      end

      for j = 1 : nwall(l2)
        rmatrx(i0+1,j0+j) = kron;
        z1 = exp ( ( wall(1,l1) - zcr(j,l2) ) * pidp );
        z1 = z1 - 1.0 / z1;
        dum = sig2 + real ( z1 )^2 + imag ( z1 )^2;
        for i = 2 : nwall(l1)
          zi = exp ( ( wall(i,l1) - zcr(j,l2) ) * pidp );
          zz = zi - 1.0 / zi;
          rmatrx(i0+i,j0+j) = -0.25 / pi * log ( dum / ...
            ( sig2 + real ( zz )^2 + imag ( zz )^2 ) );
        end
      end

      j0 = j0 + nwall(l2);

    end

    i0 = i0 + nwall(l1);

  end
%
%  Gauss elimination.
%
  for i = 1 : matdim
    rmatrx(i,i) = 1.0 / rmatrx(i,i);
    for j = i + 1 : matdim
      rmatrx(j,i) = rmatrx(j,i) * rmatrx(i,i);
      for k = i + 1 : matdim
        rmatrx(j,k) = rmatrx(j,k) - rmatrx(j,i) * rmatrx(i,k);
      end
    end
  end

  return
end
function [ er, fp, tm ] = emitst ( )

%*****************************************************************************80
%
%% EMITST is the test program for EMIT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  nb = 5;
  nv = 1000;
  nvm = 1500;
  nw = 100;

  I = sqrt ( - 1 );
  it = 10;
  ans = 6.0088546832072;
%
%  Random initialization.
%
  f7 = 78125.0;
  t30 = 1073741824.0;
  t2 = f7 / t30;

  dpds = zeros(nw,nb);
  force = zeros(nb,1);
  nwall = zeros(nb,1);
  refpt = zeros(nb,1);
  rmom = zeros(nb,1);
  wall = zeros(nw,nb);
  zcr = zeros(nw,nb);

  for j = 1 : nb
    nwall(j) = nw;
    for i = 1 : nw
      t1 = mod ( f7 * t2, 1.0 );
      t2 = mod ( f7 * t1, 1.0 );
      wall(i,j) = t1 + t2 * I;
      t1 = mod ( f7 * t2, 1.0 );
      t2 = mod ( f7 * t1, 1.0 );
      zcr(i,j) = t1 + t2 * I;
      dpds(i,j) = 0.0;
    end
  end

  rmatrx = zeros ( nw * nb, nw * nb );
  for j = 1 : nw * nb
    rmatrx(j,j) = 1.0;
    for i = 1 : j - 1
      t2 = mod ( f7 * t2, 1.0 );
      rmatrx(i,j) = 0.001 * t2;
      rmatrx(j,i) = 0.001 * t2;
    end
  end

  gamma = zeros ( nvm, 1 );
  z = zeros ( nvm, 1 );

  for i = 1 : nvm
    t1 = mod ( f7 * t2, 1.0 );
    t2 = mod ( f7 * t1, 1.0 );
    z(i) = t1 + t2 * I;
    t2 = mod ( f7 * t2, 1.0 );
    gamma(i) = t2;
  end
%
%  Timing.
%
  tic;

  for i = 1 : it
    rhs = emit ( nb, nvm, nw, dpds, force, gamma, nwall, refpt, rmatrx, ...
      rmom, wall, z, zcr );
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( rhs(19) - ans ) / ans );
  fp = it * ( 56 * nv + nb * nw * ( 97 + 44 * nv + 2 * nb * nw ) );

  return
end
function rhs = emit ( nb, nvm, nw, dpds, force, gamma, nwall, refpt, rmatrx, ...
  rmom, wall, z, zcr )

%*****************************************************************************80
%
%% 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:
%
%    04 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  I = sqrt ( - 1 );
  period = 3.0;
  sig2 = 3.0;
  u0 = 4.0;
  matdim = 500;
  delt = 1.0;
  chord = 5.0;
  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;

  expz = zeros(nvm,1);
  expmz = zeros(nvm,1);

  for i = 1 : nv
    expz(i) = exp ( z(i) * pidp );
    expmz(i) = 1.0 / expz(i);
  end

  i0 = 0;
  cupst = ( real ( uupstr ) )^2 + ( imag ( uupstr ) )^2;

  ps = zeros(nvm,1);
  psi = zeros(nw,1);

  for l = 1 : nb

    for k = 1 : nwall(l)

      expwkl = exp ( wall(k,l) * pidp );
      expmwk = 1.0 / expwkl;

      sps = 0.0;
      for i = 1 : nv
        dum3 = expz(i) * expmwk - expwkl * expmz(i);
        ps(i) = gamma(i) * log ( real ( dum3 )^2 + imag ( dum3 )^2 + sig2 );
        sps = sps + ps(i);
      end

      psi(k) = imag ( wall(k,l) * conj ( uupstr + u0 * I ) ) ...
        - sps * 0.25 / pi;
    end
%
%  Compute the right-hand side.
%
    for k = 1 : nwall(l)
      rhs(i0+k) = psi(k) - psi(1);
    end

    i0 = i0 + nwall(l);

  end
%
%  Solve the system.
%
  for i = 1 : matdim
    for j = i + 1 : matdim
      rhs(j) = rhs(j) - rmatrx(j,i) * rhs(i);
    end
  end

  for i = matdim : -1 : 1
    rhs(i) = rmatrx(i,i) * rhs(i);
    for j = 1 : i - 1
      rhs(j) = rhs(j) - rmatrx(j,i) * rhs(i);
    end
  end
%
%  Create new vortices.
%
  nolld = nv;
  i0 = 0;

  cp = zeros(nw,nb);

  for l = 1 : nb

    for k = 1 : nwall(l)
%
%  Put the new vortex at the end of the array.
%
      nv = nv + 1;
      z(nv) = zcr(k,l);
      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) * ( real ( z(nv) - refpt(l) )^2 ...
        + imag ( z(nv) - refpt(l) )^2 );
      dpds(k,l) = dpds(k,l) - gamma(nv);

    end
%
%  Filter and integrate pressure gradient to get pressure.
%
    cp(1,l) = 0.0;
    cpm = -1.0E+30;

    for k = 2 : nwall(l)
      cp(k,l) = cp(k-1,l) + ( 3.0 * ( dpds(k,l) + dpds(k-1,l) ) ...
        + dpds(1+mod(k,nwall(l)),l) ...
        + dpds(1+mod(k+nwall(l)-3,nwall(l)),l) ) / ( 4.0 * delt * cupst );
      cpm = max ( cpm, cp(k,l) );
    end
%
%  Normalize the pressure.
%
    for k = 1 : nwall(l)
      cp(k,l) = cp(k,l) - cpm;
    end
%
%  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^2 * cupst );

    i0 = i0 + nwall(l);

  end

  return
end
function [ er, fp, tm ] = vpetst ( )

%*****************************************************************************80
%
%% VPETST is the test program for VPENTA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  jl = 1;
  ju = 128;
  kl = 1;
  ku = 128;
  nja = 128;
  njb = 128;

  it = 400;
  ans = -0.354649411858726;
  lf = nja * njb * 3;
%
%  Random initialization.
%
  f7 = 78125.0;
  t30 = 1073741824.0;
  t = f7 / t30;

  a = zeros(nja,njb);
  b = zeros(nja,njb);
  c = zeros(nja,njb);
  d = zeros(nja,njb);
  e = zeros(nja,njb);
  fx = zeros(nja,njb,3);

  for j = kl : ku
    for i = jl : ju
      t = mod ( f7 * t, 1.0 );
      a(i,j) = t;
      t = mod ( f7 * t, 1.0 );
      b(i,j) = t;
      t = mod ( f7 * t, 1.0 );
      c(i,j) = t;
      t = mod ( f7 * t, 1.0 );
      d(i,j) = t;
      t = mod ( f7 * t, 1.0 );
      e(i,j) = t;
      for k = 1 : 3
        t = mod ( f7 * t, 1.0 );
        fx(i,j,k) = t;
      end
    end
  end
%
%  Timing.
%
  tic;

  for i = 1 : it
    f(1:nja,1:njb,1:3) = fx(1:nja,1:njb,1:3);
    f = vpenta ( jl, ju, kl, ku, nja, njb, a, b, c, d, e, f );
  end

  tm = toc;
%
%  Results.
%
  er = abs ( ( f(19,19,1) - ans ) / ans );
  fp = it * ku * ( 40 * ku - 53 );

  return
end
function f = vpenta ( jl, ju, kl, ku, nja, njb, a, b, c, d, e, f )

%*****************************************************************************80
%
%% VPENTA inverts 3 pentadiagonal systems simultaneously.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2010
%
%  Author:
%
%    Original FORTRAN77 version by David Bailey.
%    MATLAB version by John Burkardt.
%
  x = zeros(nja,njb);
  y = zeros(nja,njb);
%
%  Start forward generation process and sweep.
%
  j = jl;
  for k = kl : ku
    rld = c(j,k);
    rldi = 1.0 / rld;
    f(j,k,1) = f(j,k,1) * rldi;
    f(j,k,2) = f(j,k,2) * rldi;
    f(j,k,3) = f(j,k,3) * rldi;
    x(j,k) = d(j,k) * rldi;
    y(j,k) = e(j,k) * rldi;
  end

  j = jl + 1;
  for k = kl : ku
    rld1 = b(j,k);
    rld = c(j,k) - rld1 * x(j-1,k);
    rldi = 1.0 / rld;
    f(j,k,1) = ( f(j,k,1) - rld1 * f(j-1,k,1) ) * rldi;
    f(j,k,2) = ( f(j,k,2) - rld1 * f(j-1,k,2) ) * rldi;
    f(j,k,3) = ( f(j,k,3) - rld1 * f(j-1,k,3) ) * rldi;
    x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi;
    y(j,k) = e(j,k) * rldi;
  end

  for j = jl + 2 : ju - 2
    for k = kl : ku
      rld2 = a(j,k);
      rld1 = b(j,k) - rld2 * x(j-2,k);
      rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) );
      rldi = 1.0 / rld;
      f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) - rld1 * f(j-1,k,1) ) * rldi;
      f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) - rld1 * f(j-1,k,2) ) * rldi;
      f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) - rld1 * f(j-1,k,3) ) * rldi;
      x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi;
      y(j,k) = e(j,k) * rldi;
    end
  end

  j = ju - 1;
  for k = kl : ku
    rld2 = a(j,k);
    rld1 = b(j,k) - rld2 * x(j-2,k);
    rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) );
    rldi = 1.0 / rld;
    f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) - rld1 * f(j-1,k,1) ) * rldi;
    f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) - rld1 * f(j-1,k,2) ) * rldi;
    f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) - rld1 * f(j-1,k,3) ) * rldi;
    x(j,k) = ( d(j,k) - rld1 * y(j-1,k) ) * rldi;
  end

  j = ju;
  for k = kl : ku
    rld2 = a(j,k);
    rld1 = b(j,k) - rld2 * x(j-2,k);
    rld = c(j,k) - ( rld2 * y(j-2,k) + rld1 * x(j-1,k) );
    rldi = 1.0 / rld;
    f(j,k,1) = ( f(j,k,1) - rld2 * f(j-2,k,1) - rld1 * f(j-1,k,1) ) * rldi;
    f(j,k,2) = ( f(j,k,2) - rld2 * f(j-2,k,2) - rld1 * f(j-1,k,2) ) * rldi;
    f(j,k,3) = ( f(j,k,3) - rld2 * f(j-2,k,3) - rld1 * f(j-1,k,3) ) * rldi;
  end
%
%  Back sweep solution.
%
  for k = kl : ku
    f(ju,k,1) = f(ju,k,1);
    f(ju,k,2) = f(ju,k,2);
    f(ju,k,3) = f(ju,k,3);
    f(ju-1,k,1) = f(ju-1,k,1) - x(ju-1,k) * f(ju,k,1);
    f(ju-1,k,2) = f(ju-1,k,2) - x(ju-1,k) * f(ju,k,2);
    f(ju-1,k,3) = f(ju-1,k,3) - x(ju-1,k) * f(ju,k,3);
  end

  for j = 2 : ju - jl
    jx = ju - j;
    for k = kl : ku
      f(jx,k,1) = f(jx,k,1) - x(jx,k) * f(jx+1,k,1) - y(jx,k) * f(jx+2,k,1);
      f(jx,k,2) = f(jx,k,2) - x(jx,k) * f(jx+1,k,2) - y(jx,k) * f(jx+2,k,2);
      f(jx,k,3) = f(jx,k,3) - x(jx,k) * f(jx+1,k,3) - y(jx,k) * f(jx+2,k,3);
    end
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
