function test_inverse ( )

%*****************************************************************************80
%
%% test_inverse() tests the inverse computations.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 October 2022
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test_inverse():\n' );
  fprintf ( 1, '  A = a test matrix;\n' );
  fprintf ( 1, '  B = inverse as returned by test_matrix().\n' );
  fprintf ( 1, '  C = inverse as computed by Matlab''s inv() function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  |A|    = Frobenius norm of A.\n' );
  fprintf ( 1, '  |C|    = Frobenius norm of C.\n' );
  fprintf ( 1, '  |I-AC| = Frobenius norm of I-A*C.\n' );
  fprintf ( 1, '  |I-AB| = Frobenius norm of I-A*B.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Title                    N      ' );
  fprintf ( 1, '|A|      |C|      |I-AC|        |I-AB|\n' );
  fprintf ( 1, '\n' );
%
%  aegerter
%
  title = 'aegerter';
  n = 5;
  a = aegerter_matrix ( n );
  b = aegerter_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  anticirculant
%
  title = 'anticirculant';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = anticirculant_matrix ( n, n, x );
  b = anticirculant_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  antisummation
%
  title = 'antisummation';
  n = 5;
  a = antisummation_matrix ( n, n );
  b = antisummation_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  bab
%
  title = 'bab';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo ) * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = bab_matrix ( n, alpha, beta );
  b = bab_inverse ( n, alpha, beta );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  bauer
%
  title = 'bauer';
  n = 6;
  a = bauer_matrix ( );
  b = bauer_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  bernstein
%
  title = 'bernstein';
  n = 5;
  a = bernstein_matrix ( n );
  b = bernstein_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  bis
%
  title = 'bis';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = bis_matrix ( alpha, beta, n, n );
  b = bis_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  biw
%
  title = 'biw';
  n = 5;
  a = biw_matrix ( n );
  b = biw_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  bodewig
%
  title = 'bodewig';
  n = 4;
  a = bodewig_matrix ( );
  b = bodewig_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  boothroyd
%
  title = 'boothroyd';
  n = 5;
  a = boothroyd_matrix ( n );
  b = boothroyd_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  borderband
%
  title = 'borderband';
  n = 5;
  a = borderband_matrix ( n );
  b = borderband_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  carry
%
  title = 'carry';
  n = 5;
  i4_lo = 2;
  i4_hi = 20;
  key = 123456789;
  rng ( key );
  k = randi ( [ i4_lo, i4_hi ], 1 );
  a = carry_matrix ( n, k );
  b = carry_inverse ( n, k );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  cauchy
%
  title = 'cauchy';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  y = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = cauchy_matrix ( n, x, y );
  b = cauchy_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  cheby_t
%
  title = 'cheby_t';
  n = 5;
  a = cheby_t_matrix ( n );
  b = cheby_t_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  cheby_u
%
  title = 'cheby_u';
  n = 5;
  a = cheby_u_matrix ( n );
  b = cheby_u_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  cheby_van2
%
  title = 'cheby_van2';
  n = 5;
  a = cheby_van2_matrix ( n );
  b = cheby_van2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  cheby_van3
%
  title = 'cheby_van3';
  n = 5;
  a = cheby_van3_matrix ( n );
  b = cheby_van3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  chow
%
  title = 'chow';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = chow_matrix ( alpha, beta, n, n );
  b = chow_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  circulant
%
  title = 'circulant';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = circulant_matrix ( n, n, x );
  b = circulant_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  circulant2
%
  title = 'circulant2';
  n = 5;
  a = circulant2_matrix ( n );
  b = circulant2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  clement1
%  N must be even.
%
  title = 'clement1';
  n = 6;
  a = clement1_matrix ( n );
  b = clement1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  clement2
%
  title = 'clement2';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n - 1, 1 );
  y = r8_lo + ( r8_hi - r8_lo ) * rand ( n - 1, 1 );
  a = clement2_matrix ( n, x, y );
  b = clement2_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  combin
%
  title = 'combin';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = combin_matrix ( alpha, beta, n );
  b = combin_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  companion
%
  title = 'companion';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = companion_matrix ( n, x );
  b = companion_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  complex_i
%
  title = 'complex_i';
  n = 2;
  a = complex_i_matrix ( );
  b = complex_i_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  conex1
%
  title = 'conex1';
  n = 4;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = conex1_matrix ( alpha );
  b = conex1_inverse ( alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  conex2
%
  title = 'conex2';
  n = 3;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = conex2_matrix ( alpha );
  b = conex2_inverse ( alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  conex3
%
  title = 'conex3';
  n = 5;
  a = conex3_matrix ( n );
  b = conex3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  conference
%
  title = 'conference';
  n = 6;
  a = conference_matrix ( n );
  b = conference_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub2
%
  title = 'daub2';
  n = 4;
  a = daub2_matrix ( n );
  b = daub2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub4
%
  title = 'daub4';
  n = 8;
  a = daub4_matrix ( n );
  b = daub4_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub6
%
  title = 'daub6';
  n = 12;
  a = daub6_matrix ( n );
  b = daub6_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub8
%
  title = 'daub8';
  n = 16;
  a = daub8_matrix ( n );
  b = daub8_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub10
%
  title = 'daub10';
  n = 20;
  a = daub10_matrix ( n );
  b = daub10_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  daub12
%
  title = 'daub12';
  n = 24;
  a = daub12_matrix ( n );
  b = daub12_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  defective
%
  title = 'defective';
  n = 2;
  a = defective_matrix ( );
  b = defective_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  diagonal
%
  title = 'diagonal';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = diagonal_matrix ( n, n, x );
  b = diagonal_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  dif1
%
  title = 'dif1';
  n = 6;
  a = dif1_matrix ( n, n );
  b = dif1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  dif2
%
  title = 'dif2';
  n = 5;
  a = dif2_matrix ( n, n );
  b = dif2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  dorr
%
  title = 'dorr';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = dorr_matrix ( alpha, n );
  b = dorr_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  downshift
%
  title = 'downshift';
  n = 5;
  a = downshift_matrix ( n );
  b = downshift_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  drmac_bad
%
  title = 'drmac_bad';
  n = 3;
  a = drmac_bad_matrix ( );
  b = drmac_bad_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  drmac_good
%
  title = 'drmac_good';
  n = 3;
  a = drmac_good_matrix ( );
  b = drmac_good_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  eulerian
%
  title = 'eulerian';
  n = 5;
  a = eulerian_matrix ( n, n );
  b = eulerian_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  exchange
%
  title = 'exchange';
  n = 5;
  a = exchange_matrix ( n, n );
  b = exchange_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  fibonacci2
%
  title = 'fibonacci2';
  n = 5;
  a = fibonacci2_matrix ( n );
  b = fibonacci2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  fibonacci3
%
  title = 'fibonacci3';
  n = 5;
  a = fibonacci3_matrix ( n );
  b = fibonacci3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  fiedler
%
  title = 'fiedler';
  n = 7;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  x = sort ( x );
  a = fiedler_matrix ( n, n, x );
  b = fiedler_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  forsythe
%
  title = 'forsythe';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = forsythe_matrix ( alpha, beta, n );
  b = forsythe_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  fourier_cosine
%
  title = 'fourier_cosine';
  n = 5;
  a = fourier_cosine_matrix ( n );
  b = fourier_cosine_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  fourier_sine
%
  title = 'fourier_sine';
  n = 5;
  a = fourier_sine_matrix ( n );
  b = fourier_sine_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  frank
%
  title = 'frank';
  n = 5;
  a = frank_matrix ( n );
  b = frank_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  gfpp
%
  title = 'gfpp';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = gfpp_matrix ( n, alpha );
  b = gfpp_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  givens
%
  title = 'givens';
  n = 5;
  a = givens_matrix ( n, n );
  b = givens_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  gk323
%
  title = 'gk323';
  n = 5;
  a = gk323_matrix ( n, n );
  b = gk323_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  gk324
%
  title = 'gk324';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n - 1, 1 );
  a = gk324_matrix ( n, n, x );
  b = gk324_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  golub
%
  title = 'golub';
  n = 5;
  key = 123456789;
  a = golub_matrix ( n, key );
  b = golub_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hankel
%
  title = 'hankel';
  n = 5;
  key = 123456789;
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( 2 * n + 1, 1 );
  a = hankel_matrix ( n, x );
  b = hankel_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hankel_n
%
  title = 'hankel_n';
  n = 6;
  a = hankel_n_matrix ( n );
  b = hankel_n_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hanowa
%
  title = 'hanowa';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = hanowa_matrix ( alpha, n );
  b = hanowa_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  harman
%
  title = 'harman';
  n = 8;
  a = harman_matrix ( );
  b = harman_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hartley
%
  title = 'hartley';
  n = 5; 
  a = hartley_matrix ( n );
  b = hartley_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  helmert
%
  title = 'helmert';
  n = 5; 
  a = helmert_matrix ( n );
  b = helmert_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  helmert2
%
  title = 'helmert2';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = helmert2_matrix ( n, x );
  b = helmert2_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hermite
%
  title = 'hermite';
  n = 5; 
  a = hermite_matrix ( n );
  b = hermite_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  herndon
%
  title = 'herndon';
  n = 5; 
  a = herndon_matrix ( n );
  b = herndon_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  hilbert
%
  title = 'hilbert';
  n = 5; 
  a = hilbert_matrix ( n, n );
  b = hilbert_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  householder
%
  title = 'householder';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = householder_matrix ( n, x );
  b = householder_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  identity
%
  title = 'identity';
  n = 5;
  a = identity_matrix ( n, n );
  b = identity_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ill3
%
  title = 'ill3';
  n = 3;
  a = ill3_matrix ( );
  b = ill3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  integration
%
  title = 'integration';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = integration_matrix ( alpha, n );
  b = integration_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  involutory
%
  title = 'involutory';
  n = 5;
  a = involutory_matrix ( n );
  b = involutory_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  jacobi
%
  title = 'jacobi';
  n = 6;
  a = jacobi_matrix ( n, n );
  b = jacobi_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  jordan
%
  title = 'jordan';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = jordan_matrix ( n, n, alpha );
  b = jordan_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  kahan
%
  title = 'kahan';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = kahan_matrix ( alpha, n, n );
  b = kahan_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  kershaw
%
  title = 'kershaw';
  n = 4;
  a = kershaw_matrix ( );
  b = kershaw_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  kershawtri
%
  title = 'kershawtri';
  n = 5;
  x_n = floor ( ( n + 1 ) / 2 );
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( x_n, 1 );
  a = kershawtri_matrix ( n, x );
  b = kershawtri_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  kms
%
  title = 'kms';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = kms_matrix ( alpha, n, n );
  b = kms_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  laguerre
%
  title = 'laguerre';
  n = 5;
  a = laguerre_matrix ( n );
  b = laguerre_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  legendre
%
  title = 'legendre';
  n = 5;
  a = legendre_matrix ( n );
  b = legendre_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  lehmer
%
  title = 'lehmer';
  n = 5;
  a = lehmer_matrix ( n, n );
  b = lehmer_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  lesp
%
  title = 'lesp';
  n = 5;
  a = lesp_matrix ( n, n );
  b = lesp_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  lietzke
%
  title = 'lietzke';
  n = 5;
  a = lietzke_matrix ( n );
  b = lietzke_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  line_adj
%
  title = 'line_adj';
  n = 6;
  a = line_adj_matrix ( n );
  b = line_adj_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  lotkin
%
  title = 'lotkin';
  n = 5;
  a = lotkin_matrix ( n, n );
  b = lotkin_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  maxij
%
  title = 'maxij';
  n = 5;
  a = maxij_matrix ( n, n );
  b = maxij_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  milnes
%
  title = 'milnes';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = milnes_matrix ( n, n, x );
  b = milnes_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  minij
%
  title = 'minij';
  n = 5;
  a = minij_matrix ( n, n );
  b = minij_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  moler1
%
  title = 'moler1';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = moler1_matrix ( alpha, n, n );
  b = moler1_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  moler3
%
  title = 'moler3';
  n = 5;
  a = moler3_matrix ( n, n );
  b = moler3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ortega
%
  title = 'ortega';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  v1 = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  v2 = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  v3 = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = ortega_matrix ( n, v1, v2, v3 );
  b = ortega_inverse ( n, v1, v2, v3 );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  orthogonal_random
%
  title = 'orthogonal_random';
  n = 5;
  key = 123456789;
  a = orthogonal_random_matrix ( n, key );
  b = orthogonal_random_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  orthogonal_symmetric
%
  title = 'orthogonal_symmetric';
  n = 5;
  a = orthogonal_symmetric_matrix ( n );
  b = orthogonal_symmetric_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  oto
%
  title = 'oto';
  n = 5;
  a = oto_matrix ( n, n );
  b = oto_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  parter
%
  title = 'parter';
  n = 5;
  a = parter_matrix ( n, n );
  b = parter_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  pascal1
%
  title = 'pascal1';
  n = 5;
  a = pascal1_matrix ( n );
  b = pascal1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  pascal2
%
  title = 'pascal2';
  n = 5;
  a = pascal2_matrix ( n );
  b = pascal2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  pascal3
%
  title = 'pascal3';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = pascal3_matrix ( n, alpha );
  b = pascal3_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  pei
%
  title = 'pei';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = pei_matrix ( alpha, n );
  b = pei_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  permutation_random
%
  title = 'permutation_random';
  n = 5;
  key = 123456789;
  a = permutation_random_matrix ( n, key );
  b = permutation_random_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  plu
%
  title = 'plu';
  n = 5;
  key = 123456789;
  rng ( key );
  for i = 1 : n
    pivot(i) = randi ( [ i, n ], 1 );
  end
  a = plu_matrix ( n, pivot );
  b = plu_inverse ( n, pivot );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ris
%
  title = 'ris';
  n = 5;
  a = ris_matrix ( n );
  b = ris_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rodman
%
  title = 'rodman';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = rodman_matrix ( n, n, alpha );
  b = rodman_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rutis1
%
  title = 'rutis1';
  n = 4;
  a = rutis1_matrix ( );
  b = rutis1_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rutis2
%
  title = 'rutis2';
  n = 4;
  a = rutis2_matrix ( );
  b = rutis2_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rutis3
%
  title = 'rutis3';
  n = 4;
  a = rutis3_matrix ( );
  b = rutis3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rutis4
%
  title = 'rutis4';
  n = 4;
  a = rutis4_matrix ( n );
  b = rutis4_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  rutis5
%
  title = 'rutis5';
  n = 4;
  a = rutis5_matrix ( );
  b = rutis5_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  schur_block
%
  title = 'schur_block';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x_n = floor ( ( n + 1 ) / 2 );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( x_n, 1 );
  y_n = floor ( n / 2 );
  y = r8_lo + ( r8_hi - r8_lo ) * rand ( y_n, 1 );
  a = schur_block_matrix ( n, x, y );
  b = schur_block_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  spd_random
%
  title = 'spd_random';
  n = 5;
  key = 123456789;
  a = spd_random_matrix ( n, key );
  b = spd_random_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  spline
%
  title = 'spline';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = spline_matrix ( n, x );
  b = spline_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  stirling
%
  title = 'stirling';
  n = 5;
  a = stirling_matrix ( n, n );
  b = stirling_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  summation
%
  title = 'summation';
  n = 5;
  a = summation_matrix ( n, n );
  b = summation_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  sweet1
%
  title = 'sweet1';
  n = 6;
  a = sweet1_matrix ( );
  b = sweet1_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  sweet2
%
  title = 'sweet2';
  n = 6;
  a = sweet2_matrix ( );
  b = sweet2_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  sweet3
%
  title = 'sweet3';
  n = 6;
  a = sweet3_matrix ( );
  b = sweet3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  sweet4
%
  title = 'sweet4';
  n = 13;
  a = sweet4_matrix ( );
  b = sweet4_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  sylvester_kac
%
  title = 'sylvester_kac';
  n = 6;
  a = sylvester_kac_matrix ( n );
  b = sylvester_kac_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  symmetric_random
%
  title = 'symmetric_random';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = symmetric_random_matrix ( n, x, key );
  b = symmetric_random_inverse ( n, x, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  toeplitz
%
  title = 'toeplitz';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( 2 * n - 1, 1 );
  a = toeplitz_matrix ( n, x );
  b = toeplitz_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  tri_upper
%
  title = 'tri_upper';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = tri_upper_matrix ( alpha, n );
  b = tri_upper_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  tris
%
  title = 'tris';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  beta = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  gamma = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = tris_matrix ( n, n, alpha, beta, gamma );
  b = tris_inverse ( n, alpha, beta, gamma );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  triv
%
  title = 'triv';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n - 1, 1 );
  y = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  z = r8_lo + ( r8_hi - r8_lo ) * rand ( n - 1, 1 );
  a = triv_matrix ( n, x, y, z );
  b = triv_inverse ( n, x, y, z );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  triw
%
  title = 'triw';
  n = 5;
  i4_lo = 0;
  i4_hi = n - 1;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  k = randi ( [ i4_lo, i4_hi ], 1 );
  alpha = r8_lo + ( r8_hi - r8_lo )  * rand ( 1, 1 );
  a = triw_matrix ( alpha, k, n );
  b = triw_inverse ( alpha, k, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  unitary_random
%
  title = 'unitary_random';
  n = 5;
  key = 123456789;
  a = unitary_random_matrix ( n, key );
  b = unitary_random_inverse ( n, key );
  c = inv ( a );
  error_ab = c8mat_is_inverse ( n, a, b );
  error_ac = c8mat_is_inverse ( n, a, c );
  norma_frobenius = c8mat_norm_fro ( n, n, a );
  normc_frobenius = c8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  upshift
%
  title = 'upshift';
  n = 5;
  a = upshift_matrix ( n );
  b = upshift_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  vand1
%
  title = 'vand1';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = vand1_matrix ( n, x );
  b = vand1_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  vand2
%
  title = 'vand2';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  key = 123456789;
  rng ( key );
  x = r8_lo + ( r8_hi - r8_lo ) * rand ( n, 1 );
  a = vand2_matrix ( n, x );
  b = vand2_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  wilk03
%
  title = 'wilk03';
  n = 3;
  a = wilk03_matrix ( );
  b = wilk03_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  wilk04
%
  title = 'wilk04';
  n = 4;
  a = wilk04_matrix ( );
  b = wilk04_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  wilk05
%
  title = 'wilk05';
  n = 5;
  a = wilk05_matrix ( );
  b = wilk05_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  wilk21
%
  title = 'wilk21';
  n = 21;
  a = wilk21_matrix ( n );
  b = wilk21_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  wilson
%
  title = 'wilson';
  n = 4;
  a = wilson_matrix ( );
  b = wilson_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = norm ( a, 'fro' );
  normc_frobenius = norm ( c, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );

  return
end
