# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "wathen_matrix.h"

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

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  Purpose:

    wathen_matrix_test() tests wathen_matrix().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    18 February 2020

  Author:

    John Burkardt
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "wathen_matrix_test()\n" );
  printf ( "  C version\n" );
  printf ( "  test wathen_matrix().\n" );

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

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

void test01 ( )

/******************************************************************************/
/*
  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 info;
  int *ipvt;
  int job;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

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

  nx = 4;
  ny = 4;
  printf ( "  Elements in X direction NX = %d\n", nx );
  printf ( "  Elements in Y direction NY = %d\n", ny );
  printf ( "  Number of elements = %d\n", nx * ny );
/*
  Compute the number of unknowns.
*/
  n = wathen_order ( nx, ny );
  printf ( "  Number of nodes N = %d\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 = ( int * ) malloc ( n * sizeof ( int ) );
  info = dgefa ( a, n, n, ipvt );
  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "test01 - Fatal error!\n" );
    printf ( "  DGEFA returned nonzero info.\n" );
    return;
  }

  x2 = ( double * ) malloc ( n * sizeof ( double ) );
  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 );
  printf ( "  Maximum solution error is %g\n", e );
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( ipvt );
  free ( x1 );
  free ( x2 );

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

void test02 ( )

/******************************************************************************/
/*
  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 info;
  int *ipvt;
  int job;
  int lda;
  int md;
  int ml;
  int mu;
  int n;
  int nx;
  int ny;
  int seed;
  double *x1;
  double *x2;

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

  nx = 4;
  ny = 4;
  printf ( "  Elements in X direction NX = %d\n", nx );
  printf ( "  Elements in Y direction NY = %d\n", ny );
  printf ( "  Number of elements = %d\n", nx * ny );
/*
  Compute the number of unknowns.
*/
  n = wathen_order ( nx, ny );
  printf ( "  Number of nodes N = %d\n", n );
/*
  Compute the bandwidth.
*/
  wathen_bandwidth ( nx, ny, &ml, &md, &mu );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );
/*
  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 = ( int * ) malloc ( n * sizeof ( int ) );
  info = dgbfa ( a, lda, n, ml, mu, ipvt );
  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "test04 - Fatal error!\n" );
    printf ( "  DGBFA returned nonzero info.\n" );
    return;
  }
  x2 = ( double * ) malloc ( n * sizeof ( double ) );
  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 );
  printf ( "  Maximum solution error is %g\n", e );
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( ipvt );
  free ( x1 );
  free ( x2 );

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

void test05 ( )

/******************************************************************************/
/*
  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;

  printf ( "\n" );
  printf ( "test05\n" );
  printf ( "  For various problem sizes and storage schemes,\n" );
  printf ( "  measure the storage used for the Wathen system.\n" );
  printf ( "\n" );
  printf ( "                                   Predicted  Observed\n" );
  printf ( "                              GE        Band      Band      " );
  printf ( "Band    Sparse\n" );
  printf ( "    NX  Elements   Nodes   storage     width     width   " );
  printf ( "  storage   storage\n" );
  printf ( "\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.
*/
    printf ( "  %4d      %4d  %6d  %8d  %8d  %8d  %8d  %8d\n",
      nx, nx * ny, n, storage_ge, bw1, bw2, storage_gb, storage_sparse );
/*
  Ready for next iteration.
*/
    nx = nx * 2;
    ny = ny * 2;

    free ( a );
  }

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

void test06 ( )

/******************************************************************************/
/*
  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 info;
  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;

  printf ( "\n" );
  printf ( "test06():\n" );
  printf ( "  For various problem sizes,\n" );
  printf ( "  time the assembly and factorization of a Wathen system\n" );
  printf ( "  using the wathen_ge() function.\n" );
  printf ( "\n" );
  printf ( 
    "    NX  Elements   Nodes   Storage    Assembly      Factor      Error\n" );
  printf ( "\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 = ( int * ) malloc ( n * sizeof ( int ) );
    x2 = ( double * ) malloc ( n * sizeof ( double ) );
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    info = dgefa ( a, n, n, ipvt );
    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "test06 - Fatal error!\n" );
      printf ( "  DGEFA returned nonzero info.\n" );
      return;
    }
    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.
*/
    printf ( "  %4d      %4d  %6d  %8d  %10.2e  %10.2e  %10.2e\n",
      nx, nx * ny, n, storage_ge, t1, t2, e );
/*
  Ready for next iteration.
*/
    nx = nx * 2;
    ny = ny * 2;
/*
  Free memory.
*/
    free ( a );
    free ( b );
    free ( ipvt );
    free ( x1 );
    free ( x2 );
  }

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

void test07 ( )

/******************************************************************************/
/*
  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 info;
  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;

  printf ( "\n" );
  printf ( "test07():\n" );
  printf ( "  For various problem sizes,\n" );
  printf ( "  time the assembly and factorization of a Wathen system\n" );
  printf ( "  using the wathen_gb() function.\n" );
  printf ( "\n" );
  printf ( "    NX  Elements   Nodes   Storage    Assembly      Factor      Error\n" );
  printf ( "\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 = ( int * ) malloc ( n * sizeof ( int ) );
    x2 = ( double * ) malloc ( n * sizeof ( double ) );
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    info = dgbfa ( a, lda, n, ml, mu, ipvt );
    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "test07 - Fatal error!\n" );
      printf ( "  DGBFA returned nonzero info.\n" );
      return;
    }
    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.
*/
    printf ( "  %4d      %4d  %6d  %8d  %10.2e  %10.2e  %10.2e\n",
      nx, nx * ny, n, storage_gb, t1, t2, e );
/*
  Ready for next iteration.
*/
    nx = nx * 2;
    ny = ny * 2;
/*
  Free memory.
*/
    free ( a );
    free ( b );
    free ( ipvt );
    free ( x1 );
    free ( x2 );
  }

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

void test08 ( )

/******************************************************************************/
/*
  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 info;
  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;

  printf ( "\n" );
  printf ( "test08():\n" );
  printf ( "  For various problem sizes,\n" );
  printf ( "  time the assembly and factorization of a Wathen system\n" );
  printf ( "  wathen_ge() and wathen_gb().\n" );
  printf ( "\n" );
  printf ( "                   NX  Elements   Nodes   Storage    " );
  printf ( "  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 = ( int * ) malloc ( n * sizeof ( int ) );
    x2 = ( double * ) malloc ( n * sizeof ( double ) );
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    info = dgefa ( a, n, n, ipvt );
    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "wathen_GE - Fatal error.\n" );
      printf ( "  DGEFA returned nonzero info.\n" );
      return;
    }
    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.
*/
    printf ( "\n" );
    printf ( 
      "  wathen_GE      %4d      %4d  %6d  %8d  %10.2e  %10.2e  %10.2e\n",
      nx, nx * ny, n, storage_ge, t1, t2, e );
/*
  Free memory.
*/
    free ( a );
    free ( b );
    free ( ipvt );
    free ( x1 );
    free ( 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 = ( int * ) malloc ( n * sizeof ( int ) );
    x2 = ( double * ) malloc ( n * sizeof ( double ) );
    for ( i = 0; i < n; i++ )
    {
      x2[i] = b[i];
    }
    job = 0;

    t0 = cpu_time ( );
    info = 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.
*/
    printf ( 
      "  wathen_GB      %4d      %4d  %6d  %8d  %10.2e  %10.2e  %10.2e\n",
      nx, nx * ny, n, storage_gb, t1, t2, e );
/*
  Free memory.
*/
    free ( a );
    free ( b );
    free ( ipvt );
    free ( x1 );
    free ( x2 );
/*
  Ready for next iteration.
*/
    nx = nx * 2;
    ny = ny * 2;
  }

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

void test10 ( )

/******************************************************************************/
/*
  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;

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

  nx = 1;
  ny = 1;
  printf ( "  Elements in X direction NX = %d\n", nx );
  printf ( "  Elements in Y direction NY = %d\n", ny );
  printf ( "  Number of elements = %d\n", nx * ny );
/*
  Compute the number of unknowns.
*/
  n = wathen_order ( nx, ny );
  printf ( "  Number of nodes N = %d\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 = ( double * ) malloc ( n * sizeof ( double ) );
  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 );
  printf ( "  Maximum solution error is %g\n", e ); 
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( x1 );
  free ( x2 );

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

void test11 ( )

/******************************************************************************/
/*
  Purpose:

    test11() assemble, factor and solve 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;

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

  nx = 1;
  ny = 1;
  printf ( "  Elements in X direction NX = %d\n", nx );
  printf ( "  Elements in Y direction NY = %d\n", ny );
  printf ( "  Number of elements = %d\n", nx * ny );
/*
  Compute the number of unknowns.
*/
  n = wathen_order ( nx, ny );
  printf ( "  Number of nodes N = %d\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 );
  printf ( "  Number of nonzeros NZ_NUM = %d\n", nz_num );
/*
  Compute the matrix.
*/
  seed = 123456789;
  row = ( int * ) malloc ( nz_num * sizeof ( int ) );
  col = ( int * ) malloc ( nz_num * sizeof ( int ) );
  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 = ( double * ) malloc ( n * sizeof ( double ) );
  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 );
  printf ( "  Maximum solution error is %g\n", e );
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( col );
  free ( row );
  free ( x1 );
  free ( x2 );

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

void test115 ( )

/******************************************************************************/
/*
  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;

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

  nx = 4;
  ny = 4;
  printf ( "  Elements in X direction NX = %d\n", nx );
  printf ( "  Elements in Y direction NY = %d\n", ny );
  printf ( "  Number of elements = %d\n", nx * ny );
/*
  Compute the number of unknowns.
*/
  n = wathen_order ( nx, ny );
  printf ( "  Number of nodes N = %d\n", n );
/*
  Compute the bandwidth.
*/
  wathen_bandwidth ( nx, ny, &ml, &md, &mu );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );
/*
  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 = ( double * ) malloc ( n * sizeof ( double ) );
  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 );
  printf ( "  Maximum solution error is %g\n", e );
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( x1 );
  free ( x2 );

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

void wathen_xy_test ( )

/******************************************************************************/
/*
  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;

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

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

  printf ( "\n" );
  printf ( "   k   i   j         x          y\n" );
  printf ( "\n" );

  k = 0;
  for ( j = 0; j <= 2 * ny; j++ )
  {
    if ( ( j % 2 ) == 0 )
    {
      for ( i = 0; i <= 2 * ny; i++ )
      {
        printf ( "  %2d  %2d  %2d  %8.4f  %8.4f\n", k, i, j, xy[k], xy[k+n] );
        k = k + 1;
      }
    }
    else
    {
      for ( i = 0; i <= ny; i++ )
      {
        printf ( "  %2d  %2d  %2d  %8.4f  %8.4f\n", k, i, j, xy[k], xy[k+n] );
        k = k + 1;
      }
    }
  }
 
  free ( xy );

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

double cpu_time ( )

/******************************************************************************/
/*
  Purpose:

    cpu_time() returns the current reading on the CPU clock.

  Discussion:

    The CPU time measurements available through this routine are often
    not very accurate.  In some cases, the accuracy is no better than
    a hundredth of a second.  

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    06 June 2005

  Author:

    John Burkardt

  Parameters:

    Output, double CPU_TIME, the current reading of the CPU clock, in seconds.
*/
{
  double value;

  value = ( double ) clock ( ) 
        / ( double ) CLOCKS_PER_SEC;

  return value;
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

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

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

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}