# include # include # include "sparse_display.h" int main ( ); void test01 ( ); void test02 ( ); void test03 ( ); double *r8mat_uniform_01_new ( int m, int n, int *seed ); double *wathen_ge ( int nx, int ny, int n, int *seed ); int wathen_order ( int nx, int ny ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: MAIN is the main program for SPARSE_DISPLAY_TEST. Discussion: SPARSE_DISPLAY_TEST tests the SPARSE_DISPLAY library. Licensing: This code is distributed under the GNU LGPL license. Modified: 16 September 2014 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "SPARSE_DISPLAY_TEST\n" ); printf ( " C version\n" ); printf ( " Test the SPARSE_DISPLAY library.\n" ); test01 ( ); test02 ( ); test03 ( ); /* Terminate. */ printf ( "\n" ); printf ( "SPARSE_DISPLAY_TEST\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void test01 ( ) /******************************************************************************/ /* Purpose: TEST01 tests SPY_GE for a general storage matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 09 June 2014 Author: John Burkardt */ { double *a; char header[] = "wathen_ge"; int n; int nx; int ny; int seed; printf ( "\n" ); printf ( "TEST01\n" ); printf ( " SPY_GE generates a sparsity plot for a matrix stored\n" ); printf ( " in general (GE) format.\n" ); nx = 5; ny = 5; n = wathen_order ( nx, ny ); seed = 123456789; a = wathen_ge ( nx, ny, n, &seed ); spy_ge ( n, n, a, header ); free ( a ); return; } /******************************************************************************/ void test02 ( ) /******************************************************************************/ /* Purpose: TEST02 tests SPY_FILE in cases where the indices are in a file. Discussion: The files used in this example actually use negative column indices because they were output by DEAL.II and intended to be passed directly into GNUPLOT without any clever commands. So my own "SPY_FILE" is currently set up to deal exactly with such files, and hence, given sensible data will actually show a spy plot that is transposed - still readable and all, but wrong way round. Licensing: This code is distributed under the GNU LGPL license. Modified: 09 June 2014 Author: John Burkardt */ { printf ( "\n" ); printf ( "TEST02\n" ); printf ( " SPY_FILE generates a sparsity plot for a matrix for\n" ); printf ( " which there exists a file containing all the pairs\n" ); printf ( " (I,J) corresponding to nonzero entries.\n" ); spy_file ( "before", "before_data.txt" ); spy_file ( "after", "after_data.txt" ); return; } /******************************************************************************/ void test03 ( ) /******************************************************************************/ /* Purpose: TEST03 tests SPY_GE for a general storage matrix. Discussion: It is not clear to me whether the plot being created is correctly oriented. I might be seeing the transpose of the matrix. One way to check is to set up a moderate sized, highly asymmetric matrix. In this case, I will create a certain upper triangular matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 16 September 2014 Author: John Burkardt */ { double *a; char header[] = "20x30"; int i; int j; int m; int n; printf ( "\n" ); printf ( "TEST03\n" ); printf ( " SPY_GE generates a sparsity plot for a matrix stored\n" ); printf ( " in general (GE) format.\n" ); printf ( " Just to orient ourselves, generate an upper triangular matrix.\n" ); m = 20; n = 30; a = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { if ( j < i || ( ( i - j ) % ( i + 1 ) ) != 0 ) { a[i+j*m] = 0.0; } else { a[i+j*m] = 1.0; } } } spy_ge ( m, n, a, header ); free ( a ); return; } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_UNIFORM_01_NEW fills an R8MAT with pseudorandom values scaled to [0,1]. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the GNU LGPL license. Modified: 30 June 2009 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns. Input/output, int *SEED, the "seed" value. Normally, this value should not be 0, otherwise the output value of SEED will still be 0, and R8_UNIFORM will be 0. On output, SEED has been updated. Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; int k; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10; } } return r; } /******************************************************************************/ double *wathen_ge ( int nx, int ny, int n, int *seed ) /******************************************************************************/ /* Purpose: WATHEN_GE returns the Wathen matrix as a general storage (GE) matrix. Discussion: The Wathen matrix is a finite element matrix which is sparse. The entries of the matrix depend in part on a physical quantity related to density. That density is here assigned random values between 0 and 100. The matrix order N is determined by the input quantities NX and NY, which would usually be the number of elements in the X and Y directions. The value of N is N = 3*NX*NY + 2*NX + 2*NY + 1, The matrix is the consistent mass matrix for a regular NX by NY grid of 8 node serendipity elements. Here is an illustration for NX = 3, NY = 2: 23-24-25-26-27-28-29 | | | | 19 20 21 22 | | | | 12-13-14-15-16-17-18 | | | | 8 9 10 11 | | | | 1--2--3--4--5--6--7 For this example, the total number of nodes is, as expected, N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29 The matrix is symmetric positive definite for any positive values of the density RHO(X,Y). Licensing: This code is distributed under the GNU LGPL license. Modified: 07 June 2014 Author: John Burkardt Reference: Nicholas Higham, Algorithm 694: A Collection of Test Matrices in MATLAB, ACM Transactions on Mathematical Software, Volume 17, Number 3, September 1991, pages 289-305. Andrew Wathen, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA Journal of Numerical Analysis, Volume 7, Number 4, October 1987, pages 449-457. Parameters: Input, int NX, NY, values which determine the size of the matrix. Input, int N, the number of rows and columns. Input/output, int *SEED, the random number seed. Output, double WATHEN_GE[N*N], the matrix. */ { double *a; const double em[8*8] = { 6.0, -6.0, 2.0, -8.0, 3.0, -8.0, 2.0, -6.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, 2.0, -6.0, 6.0, -6.0, 2.0, -8.0, 3.0, -8.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, 3.0, -8.0, 2.0, -6.0, 6.0, -6.0, 2.0, -8.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, 2.0, -8.0, 3.0, -8.0, 2.0, -6.0, 6.0, -6.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 }; int i; int ii; int j; int jj; int kcol; int krow; int node[8]; double *rho; a = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { a[i+j*n] = 0.0; } } rho = r8mat_uniform_01_new ( nx, ny, seed ); for ( j = 0; j < ny; j++ ) { for ( i = 0; i < nx; i++ ) { rho[i+j*nx] = 100.0 * rho[i+j*nx]; } } for ( j = 0; j < nx; j++ ) { for ( i = 0; i < nx; i++ ) { node[0] = 3 * ( j + 1 ) * nx + 2 * ( j + 1 ) + 2 * ( i + 1 ); node[1] = node[0] - 1; node[2] = node[0] - 2; node[3] = ( 3 * ( j + 1 ) - 1 ) * nx + 2 * ( j + 1 ) + ( i + 1 ) - 2; node[4] = ( 3 * ( j + 1 ) - 3 ) * nx + 2 * ( j + 1 ) + 2 * ( i + 1 ) - 4; node[5] = node[4] + 1; node[6] = node[4] + 2; node[7] = node[3] + 1; for ( krow = 0; krow < 8; krow++ ) { ii = node[krow]; for ( kcol = 0; kcol < 8; kcol++ ) { jj = node[kcol]; a[ii+jj*n] = a[ii+jj*n] + rho[i+j*nx] * em[krow+kcol*8]; } } } } free ( rho ); return a; } /******************************************************************************/ int wathen_order ( int nx, int ny ) /******************************************************************************/ /* Purpose: WATHEN_ORDER returns the order of the WATHEN matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 June 2014 Author: John Burkardt Reference: Nicholas Higham, Algorithm 694: A Collection of Test Matrices in MATLAB, ACM Transactions on Mathematical Software, Volume 17, Number 3, September 1991, pages 289-305. Andrew Wathen, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA Journal of Numerical Analysis, Volume 7, 1987, pages 449-457. Parameters: Input, int NX, NY, values which determine the size of A. Output, int WATHEN_ORDER, the order of the matrix, as determined by NX and NY. */ { int n; n = 3 * nx * ny + 2 * nx + 2 * ny + 1; return n; }