# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "wathen_matrix.hpp"

int main ( );
void test01 ( );
void test02 ( );
void test05 ( );
void test06 ( );
void test07 ( );
void test08 ( );
void test10 ( );
void test11 ( );
void test115 ( );
void wathen_xy_test ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    wathen_matrix_test() tests wathen_matrix().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 March 2022
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "wathen_matrix_test():\n";
  cout << "  C++ version\n";
  cout << "  test wathen_matrix()\n";

  test01 ( );
  test02 ( );
  test05 ( );
  test06 ( );
  test07 ( );
  test08 ( );
  test10 ( );
  test11 ( );
  test115 ( );
  wathen_xy_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "wathen_matrix_test()\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test01 assembles, factor and solve using wathen_ge().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int *ipvt;
  int job;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test01():\n";
  cout << "  Assemble, factor and solve a Wathen system\n";
  cout << "  defined by wathen_ge().\n";
  cout << "\n";

  nx = 4;
  ny = 4;
  cout << "  Elements in X direction NX = " << nx << "\n";
  cout << "  Elements in Y direction NY = " << ny << "\n";
  cout << "  Number of elements = " << nx * ny << "\n";
//
//  Compute the number of unknowns.
//
  n = wathen_order ( nx, ny );
  cout << "  Number of nodes N = " << n << "\n";
//
//  Set up a random solution X.
//
  seed = 123456789;
  x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
  seed = 123456789;
  a = wathen_ge ( nx, ny, n, seed );
//
//  Compute the corresponding right hand side B.
//
  b = mv_ge ( n, n, a, x1 );
//
//  Solve the linear system.
//
  ipvt = new int[n];
  dgefa ( a, n, n, ipvt );

  x2 = new double[n];
  for ( i = 0; i < n; i++ )
  {
    x2[i] = b[i];
  }
  job = 0;
  dgesl ( a, n, n, ipvt, x2, job );
//
//  Compute the maximum solution error.
//
  e = r8vec_diff_norm_li ( n, x1, x2 );
  cout << "  Maximum solution error is " << e << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] ipvt;
  delete [] x1;
  delete [] x2;

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

void test02 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test02() assembles, factors and solves using wathen_gb().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int *ipvt;
  int job;
  int lda;
  int md;
  int ml;
  int mu;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test02():\n";
  cout << "  Assemble, factor and solve a Wathen system\n";
  cout << "  using wathen_gb().\n";
  cout << "\n";

  nx = 4;
  ny = 4;
  cout << "  Elements in X direction NX = " << nx << "\n";
  cout << "  Elements in Y direction NY = " << ny << "\n";
  cout << "  Number of elements = " << nx * ny << "\n";
//
//  Compute the number of unknowns.
//
  n = wathen_order ( nx, ny );
  cout << "  Number of nodes N = " << n << "\n";
//
//  Compute the bandwidth.
//
  wathen_bandwidth ( nx, ny, ml, md, mu );
  cout << "  Lower bandwidth ML = " << ml << "\n";
  cout << "  Upper bandwidth MU = " << mu << "\n";
//
//  Set up a random solution X1.
//
  seed = 123456789;
  x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
  seed = 123456789;
  a = wathen_gb ( nx, ny, n, seed );
//
//  Compute the corresponding right hand side B.
//
  b = mv_gb ( n, n, ml, mu, a, x1 );
//
//  Solve the linear system.
//
  lda = 2 * ml + mu + 1;
  ipvt = new int[n];
  dgbfa ( a, lda, n, ml, mu, ipvt );

  x2 = new double[n];
  for ( i = 0; i < n; i++ )
  {
    x2[i] = b[i];
  }
  job = 0;
  dgbsl ( a, lda, n, ml, mu, ipvt, x2, job );
//
//  Compute the maximum solution error.
//
  e = r8vec_diff_norm_li ( n, x1, x2 );
  cout << "  Maximum solution error is " << e << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] ipvt;
  delete [] x1;
  delete [] x2;

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

void test05 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test05() measures the storage needed for the Wathen system.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  int bd1;
  int bd2;
  int bl1;
  int bl2;
  int bu1;
  int bu2;
  int bw1;
  int bw2;
  int n;
  int nx;
  int ny;
  int seed;
  int storage_gb;
  int storage_ge;
  int storage_sparse;
  int test;

  cout << "\n";
  cout << "test05():\n";
  cout << "  For various problem sizes and storage schemes,\n";
  cout << "  measure the storage used for the Wathen system.\n";
  cout << "\n";
  cout << "                                   Predicted  Observed\n";
  cout << "                              GE        Band      Band      ";
  cout << "Band    Sparse\n";
  cout << "    NX  Elements   Nodes   storage     width     width   ";
  cout << "  storage   storage\n";
  cout << "\n";

  nx = 1;
  ny = 1;

  for ( test = 1; test <= 6; test++ )
  {
//
//  Compute the number of unknowns.
//
    n = wathen_order ( nx, ny );
//
//  Predict the bandwidth.
//
    wathen_bandwidth ( nx, ny, bl1, bd1, bu1 );
    bw1 = bl1 + bd1 + bu1;
//
//  Compute the matrix.
//
    seed = 123456789;
    a = wathen_ge ( nx, ny, n, seed );

    storage_ge = n * n;

    bandwidth ( n, n, a, bw2, bl2, bd2, bu2 );
    storage_gb = ( 2 * bl2 + 1 + bu2 ) * n;

    storage_sparse = nonzeros ( n, n, a );
//
//  Report.
//
    cout << setw(6) << nx << "      "
         << setw(4) << nx * ny << "  "
         << setw(6) << n << "  "
         << setw(8) << storage_ge << "  "
         << setw(8) << bw1 << "  "
         << setw(8) << bw2 << "  "
         << setw(8) << storage_gb << "  "
         << setw(8) << storage_sparse << "\n";
//
//  Ready for next iteration.
//
    nx = nx * 2;
    ny = ny * 2;

    delete [] a;
  }

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

void test06 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test06() times wathen_ge() assembly and solution.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int *ipvt;
  int job;
  int n;
  int nx;
  int ny;
  int seed;
  int storage_ge;
  double t0;
  double t1;
  double t2;
  int test;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test06():\n";
  cout << "  For various problem sizes,\n";
  cout << "  time the assembly and factorization of a Wathen system\n";
  cout << "  using the wathen_ge() function.\n";
  cout << "\n";
  cout << 
    "    NX  Elements   Nodes   Storage    Assembly      Factor      Error\n";
  cout << "\n";

  nx = 1;
  ny = 1;

  for ( test = 1; test <= 6; test++ )
  {
//
//  Compute the number of unknowns.
//
    n = wathen_order ( nx, ny );
    storage_ge = n * n;
//
//  Set up a random solution X1.
//
    seed = 123456789;
    x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix, and measure the storage required.
//
    seed = 123456789;

    t0 = cpu_time ( );
    a = wathen_ge ( nx, ny, n, seed );
    t1 = cpu_time ( );
    t1 = t1 - t0;
//
//  Compute the corresponding right hand side B.
//
    b = mv_ge ( n, n, a, x1 );
//
//  Solve the system.
//
    ipvt = new int[n];
    x2 = new double[n];
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    dgefa ( a, n, n, ipvt );
    dgesl ( a, n, n, ipvt, x2, job );
    t2 = cpu_time ( );
    t2 = t2 - t0;
//
//  Compute the maximum solution error.
//
    e = r8vec_diff_norm_li ( n, x1, x2 );
//
//  Report.
//
    cout << setw(6) << nx << "      "
         << setw(4) << nx * ny << "  " 
         << setw(6) << n << "  "
         << setw(8) << storage_ge << "  "
         << setw(10) << t1 << "  "
         << setw(10) << t2 << "  "
         << setw(10) << e << "\n";
//
//  Ready for next iteration.
//
    nx = nx * 2;
    ny = ny * 2;
//
//  Free memory.
//
    delete [] a;
    delete [] b;
    delete [] ipvt;
    delete [] x1;
    delete [] x2;
  }

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

void test07 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test07() times wathen_gb() assembly and solution.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int *ipvt;
  int job;
  int lda;
  int md;
  int ml;
  int mu;
  int n;
  int nx;
  int ny;
  int seed;
  int storage_gb;
  double t0;
  double t1;
  double t2;
  int test;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test07()\n";
  cout << "  For various problem sizes,\n";
  cout << "  time the assembly and factorization of a Wathen system\n";
  cout << "  using the wathen_gb() function.\n";
  cout << "\n";
  cout << "    NX  Elements   Nodes   Storage    Assembly      Factor      Error\n";
  cout << "\n";

  nx = 1;
  ny = 1;

  for ( test = 1; test <= 6; test++ )
  {
//
//  Compute the number of unknowns.
//
    n = wathen_order ( nx, ny );
//
//  Compute the bandwidth.
//
    wathen_bandwidth ( nx, ny, ml, md, mu );
    storage_gb = ( 2 * ml + mu + 1 ) * n;
//
//  Set up a random solution X1.
//
    seed = 123456789;
    x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
    seed = 123456789;
    t0 = cpu_time ( );
    a = wathen_gb ( nx, ny, n, seed );
    t1 = cpu_time ( );
    t1 = t1 - t0;
//
//  Compute the corresponding right hand side B.
//
    b = mv_gb ( n, n, ml, mu, a, x1 );
//
//  Solve the system.
//
    lda = 2 * ml + mu + 1;
    ipvt = new int[n];
    x2 = new double[n];
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    dgbfa ( a, lda, n, ml, mu, ipvt );
    dgbsl ( a, lda, n, ml, mu, ipvt, x2, job );
    t2 = cpu_time ( );
    t2 = t2 - t0;
//
//  Compute the maximum solution error.
//
    e = r8vec_diff_norm_li ( n, x1, x2 );
//
//  Report.
//
    cout << setw(6) << nx << "      "
         << setw(4) << nx * ny << "  " 
         << setw(6) << n << "  "
         << setw(8) << storage_gb << "  "
         << setw(10) << t1 << "  "
         << setw(10) << t2 << "  "
         << setw(10) << e << "\n";
//
//  Ready for next iteration.
//
    nx = nx * 2;
    ny = ny * 2;
//
//  Free memory.
//
    delete [] a;
    delete [] b;
    delete [] ipvt;
    delete [] x1;
    delete [] x2;
  }

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

void test08 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test08() times wathen_ge() and wathen_gb().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int *ipvt;
  int job;
  int lda;
  int md;
  int ml;
  int mu;
  int n;
  int nx;
  int ny;
  int seed;
  int storage_gb;
  int storage_ge;
  double t0;
  double t1;
  double t2;
  int test;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test08()\n";
  cout << "  For various problem sizes,\n";
  cout << "  time the assembly and factorization of a Wathen system\n";
  cout << "  wathen_ge() and wathen_gb().\n";
  cout << "\n";
  cout << "                   NX  Elements   Nodes   Storage    ";
  cout << "  Assembly      Factor      Error\n";

  nx = 1;
  ny = 1;

  for ( test = 1; test <= 6; test++ )
  {
//
//  Compute the number of unknowns.
//
    n = wathen_order ( nx, ny );
    storage_ge = n * n;
//
//  Set up a random solution X1.
//
    seed = 123456789;
    x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
    seed = 123456789;
    t0 = cpu_time ( );
    a = wathen_ge ( nx, ny, n, seed );
    t1 = cpu_time ( );
    t1 = t1 - t0;
//
//  Compute the corresponding right hand side B.
//
    b = mv_ge ( n, n, a, x1 );
//
//  Solve the system.
//
    ipvt = new int[n];
    x2 = new double[n];
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    dgefa ( a, n, n, ipvt );
    dgesl ( a, n, n, ipvt, x2, job );
    t2 = cpu_time ( );
    t2 = t2 - t0;
//
//  Compute the maximum solution error.
//
    e = r8vec_diff_norm_li ( n, x1, x2 );
//
//  Report.
//
    cout << "\n";
    cout << "  WATHEN_GE      "
         << setw(6) << nx << "      "
         << setw(4) << nx * ny << "  " 
         << setw(6) << n << "  "
         << setw(8) << storage_ge << "  "
         << setw(10) << t1 << "  "
         << setw(10) << t2 << "  "
         << setw(10) << e << "\n";
//
//  Free memory.
//
    delete [] a;
    delete [] b;
    delete [] ipvt;
    delete [] x1;
    delete [] x2;
//
//  Compute the bandwidth.
//
    wathen_bandwidth ( nx, ny, ml, md, mu );
    storage_gb = ( 2 * ml + mu + 1 ) * n;
//
//  Set up a random solution X1.
//
    seed = 123456789;
    x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
    seed = 123456789;
    t0 = cpu_time ( );
    a = wathen_gb ( nx, ny, n, seed );
    t1 = cpu_time ( );
    t1 = t1 - t0;
//
//  Compute the corresponding right hand side B.
//
    b = mv_gb ( n, n, ml, mu, a, x1 );
//
//  Solve the system.
//
    lda = 2 * ml + mu + 1;
    ipvt = new int[n];
    x2 = new double[n];
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    dgbfa ( a, lda, n, ml, mu, ipvt );
    dgbsl ( a, lda, n, ml, mu, ipvt, x2, job );
    t2 = cpu_time ( );
    t2 = t2 - t0;
//
//  Compute the maximum solution error.
//
    e = r8vec_diff_norm_li ( n, x1, x2 );
//
//  Report.
//
    cout << "  WATHEN_GB      "
         << setw(6) << nx << "      "
         << setw(4) << nx * ny << "  " 
         << setw(6) << n << "  "
         << setw(8) << storage_gb << "  "
         << setw(10) << t1 << "  "
         << setw(10) << t2 << "  "
         << setw(10) << e << "\n";
//
//  Free memory.
//
    delete [] a;
    delete [] b;
    delete [] ipvt;
    delete [] x1;
    delete [] x2;
//
//  Ready for next iteration.
//
    nx = nx * 2;
    ny = ny * 2;
  }

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

void test10 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test10() assembles, factor and solve using wathen_ge() and cg_ge().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test10():\n";
  cout << "  Assemble, factor and solve a Wathen system\n";
  cout << "  defined by wathen_ge() and cg_ge().\n";
  cout << "\n";

  nx = 1;
  ny = 1;
  cout << "  Elements in X direction NX = " << nx << "\n";
  cout << "  Elements in Y direction NY = " << ny << "\n";
  cout << "  Number of elements = " << nx * ny << "\n";
//
//  Compute the number of unknowns.
//
  n = wathen_order ( nx, ny );
  cout << "  Number of nodes N = " << n << "\n";
//
//  Set up a random solution X.
//
  seed = 123456789;
  x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
  seed = 123456789;
  a = wathen_ge ( nx, ny, n, seed );
//
//  Compute the corresponding right hand side B.
//

  b = mv_ge ( n, n, a, x1 );
//
//  Solve the linear system.
//
  x2 = new double[n];
  for ( i = 0; i < n; i++ )
  {
    x2[i] = 1.0;
  }
  cg_ge ( n, a, b, x2 );
//
//  Compute the maximum solution error.
//
  e = r8vec_diff_norm_li ( n, x1, x2 );
  cout << "  Maximum solution error is " << e << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] x1;
  delete [] x2;

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

void test11 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test11() assembles, factors and solves using wathen_st() + cg_st().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  int *col;
  double e;
  int i;
  int n;
  int nx;
  int ny;
  int nz_num;
  int *row;
  int seed;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test11():\n";
  cout << "  Assemble, factor and solve a Wathen system\n";
  cout << "  defined by wathen_st() and cg_st().\n";
  cout << "\n";

  nx = 1;
  ny = 1;
  cout << "  Elements in X direction NX = " << nx << "\n";
  cout << "  Elements in Y direction NY = " << ny << "\n";
  cout << "  Number of elements = " << nx * ny << "\n";
//
//  Compute the number of unknowns.
//
  n = wathen_order ( nx, ny );
  cout << "  Number of nodes N = " << n << "\n";
//
//  Set up a random solution X1.
//
  seed = 123456789;
  x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix size.
//
  nz_num = wathen_st_size ( nx, ny );
  cout << "  Number of nonzeros NZ_NUM = " << nz_num << "\n";
//
//  Compute the matrix.
//
  seed = 123456789;
  row = new int[nz_num];
  col = new int[nz_num];
  a = wathen_st ( nx, ny, nz_num, seed, row, col );
//
//  Compute the corresponding right hand side B.
//
  b = mv_st ( n, n, nz_num, row, col, a, x1 );
//
//  Solve the linear system.
//
  x2 = new double[n];
  for ( i = 0; i < n; i++ )
  {
    x2[i] = 1.0;
  }
  cg_st ( n, nz_num, row, col, a, b, x2 );
//
//  Compute the maximum solution error.
//
  e = r8vec_diff_norm_li ( n, x1, x2 );
  cout << "  Maximum solution error is " << e << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] col;
  delete [] row;
  delete [] x1;
  delete [] x2;

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

void test115 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test115() assembles, factors and solves using wathen_gb() and cg_gb().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 June 2014
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double e;
  int i;
  int md;
  int ml;
  int mu;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

  cout << "\n";
  cout << "test115():\n";
  cout << "  Assemble, factor and solve a Wathen system\n";
  cout << "  using wathen_gb() and cg_gb().\n";
  cout << "\n";

  nx = 4;
  ny = 4;
  cout << "  Elements in X direction NX = " << nx << "\n";
  cout << "  Elements in Y direction NY = " << ny << "\n";
  cout << "  Number of elements = " << nx * ny << "\n";
//
//  Compute the number of unknowns.
//
  n = wathen_order ( nx, ny );
  cout << "  Number of nodes N = " << n << "\n";
//
//  Compute the bandwidth.
//
  wathen_bandwidth ( nx, ny, ml, md, mu );
  cout << "  Lower bandwidth ML = " << ml << "\n";
  cout << "  Upper bandwidth MU = " << mu << "\n";
//
//  Set up a random solution X1.
//
  seed = 123456789;
  x1 = r8vec_uniform_01_new ( n, seed );
//
//  Compute the matrix.
//
  seed = 123456789;
  a = wathen_gb ( nx, ny, n, seed );
//
//  Compute the corresponding right hand side B.
//
  b = mv_gb ( n, n, ml, mu, a, x1 );
//
//  Solve the linear system.
//
  x2 = new double[n];
  for ( i = 0; i < n; i++ )
  {
    x2[i] = 1.0;
  }
  cg_gb ( n, ml, mu, a, b, x2 );
//
//  Compute the maximum solution error.
//
  e = r8vec_diff_norm_li ( n, x1, x2 );
  cout << "  Maximum solution error is " << e << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] x1;
  delete [] x2;

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

void wathen_xy_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    wathen_xy_test() tests wathen_xy().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 February 2020
//
//  Author:
//
//    John Burkardt
//
{
  int i;
  int j;
  int k;
  int n;
  int nx;
  int ny;
  double *xy;

  cout << "\n";
  cout << "wathen_xy_test():\n";
  cout << "  wathen_xy() returns the (X,Y) coordinates of nodes in the\n";
  cout << "  Wathen finite element system.\n";

  nx = 3;
  ny = 3;
  n = wathen_order ( nx, ny );
  xy = wathen_xy ( nx, ny, n );

  cout << "\n";
  cout << "   k   i   j         x          y\n";
  cout << "\n";

  k = 0;
  for ( j = 0; j <= 2 * ny; j++ )
  {
    if ( ( j % 2 ) == 0 )
    {
      for ( i = 0; i <= 2 * ny; i++ )
      {
        cout << "  " << setw(2) << k
             << "  " << setw(2) << i
             << "  " << setw(2) << j
             << "  " << setw(8) << xy[k]
             << "  " << setw(8) << xy[k+n] << "\n";
        k = k + 1;
      }
    }
    else
    {
      for ( i = 0; i <= ny; i++ )
      {
        cout << "  " << setw(2) << k
             << "  " << setw(2) << i
             << "  " << setw(2) << j
             << "  " << setw(8) << xy[k]
             << "  " << setw(8) << xy[k+n] << "\n";
        k = k + 1;
      }
    }
  }
 
  delete [] xy;

  return;
}