# include # include # include # include using namespace std; # include "wathen.hpp" int main ( ); void test01 ( ); void test02 ( ); void test05 ( ); void test06 ( ); void test07 ( ); void test08 ( ); void test10 ( ); void test11 ( ); void test115 ( ); //****************************************************************************80 int main ( ) //****************************************************************************80 // // Purpose: // // MAIN is the main program for WATHEN_PRB. // // Discussion: // // WATHEN_PRB tests the WATHEN library. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 08 June 2014 // // Author: // // John Burkardt // { timestamp ( ); cout << "\n"; cout << "WATHEN_PRB\n"; cout << " C++ version\n"; cout << " Test the WATHEN library.\n"; test01 ( ); test02 ( ); test05 ( ); test06 ( ); test07 ( ); test08 ( ); test10 ( ); test11 ( ); test115 ( ); // // Terminate. // cout << "\n"; cout << "WATHEN_PRB\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 GNU LGPL 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; 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]; info = 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 GNU LGPL license. // // Modified: // // 08 June 2014 // // Author: // // John Burkardt // { double *a; double *b; double e; int i; int info; int *ipvt; int j; int jhi; int jlo; 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]; info = 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 GNU LGPL 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 GNU LGPL 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; 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 ( ); info = 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 GNU LGPL license. // // Modified: // // 08 June 2014 // // Author: // // John Burkardt // { double *a; double *b; double e; int i; int info; int *ipvt; int j; int jhi; int jlo; 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 ( ); 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. // 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/WATHEN_GB. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 08 June 2014 // // Author: // // John Burkardt // { double *a; double *b; double e; int i; int info; int *ipvt; int j; int jhi; int jlo; 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/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 ( ); info = 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 ( ); 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. // 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 GNU LGPL 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 assemble, factor and solve using WATHEN_ST + CG_ST. // // Licensing: // // This code is distributed under the GNU LGPL 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 GNU LGPL license. // // Modified: // // 08 June 2014 // // Author: // // John Burkardt // { double *a; double *b; double e; int i; int j; int jhi; int jlo; int lda; 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; }