function r8vec_multinormal_pdf_test ( )

%*****************************************************************************80
%
%% r8vec_multinormal_pdf_test() tests r8vec_multinormal_pdf().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 August 2015
%
%  Author:
%
%    John Burkardt
%
  initialize ( );

  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'r8vec_multinormal_pdf_test():\n' );
  fprintf ( 1, '  r8vec_multinormal_pdf() evaluates the PDF for the\n' );
  fprintf ( 1, '  multinormal distribution.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The covariance matrix is C.\n' );
  fprintf ( 1, '  The definition uses the inverse of C;\n' );
  fprintf ( 1, '  R8VEC_MULTINORMAL_PDF uses the Cholesky factor.\n' );
  fprintf ( 1, '  Verify that the algorithms are equivalent.\n' );
%
%  Generate a random upper triangular matrix with positive diagonal.
%
  r1 = zeros ( n, n );

  for j = 1 : n
    for i = 1 : j
      r1(i,j) = r8_uniform_01_sample ( );
    end
  end

  r8ut_print ( n, n, r1, '  R1:' );
%
%  Compute a positive definite symmetric covariance matrix C.
%
  c_ge = r1' * r1;

  r8ge_print ( n, n, c_ge, '  C:' );
%
%  Convert to R8PO format.
%
  c = r8ge_to_r8po ( n, c_ge );
%
%  Compute the Cholesky factor R.
%
  r2 = r8po_fa ( n, c );

  r8ut_print ( n, n, r2, '  R2:' );
%
%  Compute the determinant of C.
%
  c_det = r8po_det ( n, r2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Determinant of C = %g\n', c_det );
%
%  Compute the inverse of C.
%
  c_inv = r8po_inverse ( n, r2 );
  r8po_print ( n, c_inv, '  inverse(C):' );
%
%  Compute a random set of means.
%
  mu = zeros ( n, 1 );
  for i = 1 : n
    mu(i) = r8_normal_01_sample ( );
  end

  r8vec_print ( n, mu, '  MU:' );
%
%  Compute X as small variations from MU.
%
  x = zeros ( n, 1 );
  for i = 1 : n
    eps = 0.01 * r8_normal_01_sample ( );
    x(i) = ( 1.0 + eps ) * mu(i);
  end

  r8vec_print ( n, x, '  X:' );
%
%  Compute PDF1 from the function.
%
  pdf1 = r8vec_multinormal_pdf ( n, mu, r2, c_det, x );
%
%  Compute PDF2 from the definition.
%
  y(1:n,1) = x(1:n,1) - mu(1:n,1);

  cx = r8po_mv ( n, c_inv, y );
  xcx = y' * cx;

  pdf2 = 1.0 / sqrt ( ( 2.0 * pi ) ^ n ) ...
    * 1.0 / sqrt ( c_det ) * exp ( - 0.5 * xcx );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF1 = %g\n', pdf1 );
  fprintf ( 1, '  PDF2 = %g\n', pdf2 );

  return
end
