subroutine r8mat_cholesky_factor ( n, a, c, flag )

!*****************************************************************************80
!
!! r8mat_cholesky_factor() computes the Cholesky factor of a symmetric matrix.
!
!  Discussion:
!
!    The matrix must be symmetric and positive semidefinite.
!
!    For a positive semidefinite symmetric matrix A, the Cholesky factorization
!    is a lower triangular matrix L such that:
!
!      A = L * L'
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    08 April 2009
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N matrix.
!
!  Output:
!
!    real ( kind = rk8 ) C(N,N), the lower triangular Cholesky factor.
!
!    integer FLAG:
!    0, no error occurred.
!    1, the matrix is not positive definite.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) c(n,n)
  integer flag
  integer i
  integer j
  real ( kind = rk8 ) sum2

  flag = 0

  c(1:n,1:n) = a(1:n,1:n)

  do j = 1, n

    c(1:j-1,j) = 0.0D+00

    do i = j, n

      sum2 = c(j,i) - dot_product ( c(j,1:j-1), c(i,1:j-1) )

      if ( i == j ) then
        if ( sum2 <= 0.0D+00 ) then
          flag = 1
          return
        else
          c(i,j) = sqrt ( sum2 )
        end if
      else
        if ( c(j,j) /= 0.0D+00 ) then
          c(i,j) = sum2 / c(j,j)
        else
          c(i,j) = 0.0D+00
        end if
      end if

    end do

  end do

  return
end
subroutine r8mat_cholesky_solve ( n, l, b, x )

!*****************************************************************************80
!
!! r8mat_cholesky_solve() solves a Cholesky factored linear system A * x = b.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!    This routine works with the lower triangular Cholesky factor A = L * L'.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) L(N,N), the N by N lower Cholesky factor of the
!    system matrix A.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) b(n)
  real ( kind = rk8 ) l(n,n)
  real ( kind = rk8 ) x(n)
!
!  Solve L * y = b.
!
  call r8mat_l_solve ( n, l, b, x )
!
!  Solve L' * x = y.
!
  call r8mat_lt_solve ( n, l, x, x )

  return
end
subroutine r8mat_l_solve ( n, a, b, x )

!*****************************************************************************80
!
!! r8mat_l_solve() solves a lower triangular linear system.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) b(n)
  integer i
  real ( kind = rk8 ) x(n)
!
!  Solve L * x = b.
!
  do i = 1, n
    x(i) = ( b(i) - dot_product ( a(i,1:i-1), x(1:i-1) ) ) / a(i,i)
  end do

  return
end
subroutine r8mat_lt_solve ( n, a, b, x )

!*****************************************************************************80
!
!! r8mat_lt_solve() solves a transposed lower triangular linear system.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!    Given the lower triangular matrix A, the linear system to be solved is:
!
!      A' * x = b
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) b(n)
  integer i
  real ( kind = rk8 ) x(n)
!
!  Solve L'*x = b.
!
  do i = n, 1, -1
    x(i) = ( b(i) - dot_product ( x(i+1:n), a(i+1:n,i) ) ) / a(i,i)
  end do

  return
end
function r8vec_norm_affine ( n, v0, v1 )

!*****************************************************************************80
!
!! r8vec_norm_affine() returns the affine norm of an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    The affine vector L2 norm is defined as:
!
!      R8VEC_NORM_AFFINE(V0,V1)
!        = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 October 2010
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the vectors.
!
!    real ( kind = rk8 ) V0(N), the base vector.
!
!    real ( kind = rk8 ) V1(N), the vector whose affine norm is desired.
!
!  Output:
!
!    real ( kind = rk8 ) R8VEC_NORM_AFFINE, the L2 norm of V1-V0.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) v0(n)
  real ( kind = rk8 ) v1(n)

  r8vec_norm_affine = sqrt ( sum ( ( v0(1:n) - v1(1:n) ) ** 2 ) )

  return
end
function triangle_num ( n )

!*****************************************************************************80
!
!! triangle_num() returns the N-th triangular number.
!
!  Discussion:
!
!    The N-th triangular number T(N) is formed by the sum of the first
!    N integers:
!
!      T(N) = sum ( 1 <= I <= N ) I
!
!    By convention, T(0) = 0.
!
!    T(N) can be computed quickly by the formula:
!
!      T(N) = ( N * ( N + 1 ) ) / 2
!
!  First Values:
!
!     0
!     1
!     3
!     6
!    10
!    15
!    21
!    28
!    36
!    45
!    55
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    11 August 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the index of the desired number, 
!    which must be at least 0.
!
!    Output, integer TRIANGLE_NUM, the N-th triangular number.
!
  implicit none

  integer n
  integer triangle_num

  triangle_num = ( n * ( n + 1 ) ) / 2

  return
end
subroutine vandermonde_interp_2d_matrix ( n, m, x, y, a )

!*****************************************************************************80
!
!! vandermonde_interp_2d_matrix() computes a Vandermonde 2D interpolation matrix.
!
!  Discussion:
!
!    We assume the approximating function has the form of a polynomial
!    in X and Y of total degree M.
!
!      p(x,y) = c00 
!             + c10 * x                + c01 * y
!             + c20 * x^2   + c11 * xy + c02 * y^2
!             + ...
!             + cm0 * x^(m) + ...      + c0m * y^m.
!
!    If we let T(K) = the K-th triangular number 
!            = sum ( 1 <= I <= K ) I
!    then the number of coefficients in the above polynomial is T(M+1).
!
!    We have n data locations (x(i),y(i)) and values z(i) to approximate:
!
!      p(x(i),y(i)) = z(i)
!
!    and we assume that N = T(M+1).
!
!    This can be cast as an NxN linear system for the polynomial
!    coefficients:
!
!      [ 1 x1 y1  x1^2 ... y1^m ] [ c00 ] = [  z1 ]
!      [ 1 x2 y2  x2^2 ... y2^m ] [ c10 ] = [  z2 ]
!      [ 1 x3 y3  x3^2 ... y3^m ] [ c01 ] = [  z3 ]
!      [ ...................... ] [ ... ] = [ ... ]
!      [ 1 xn yn  xn^2 ... yn^m ] [ c0n ] = [  zn ]
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 September 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data points.  It is necessary 
!    that N = T(M+1), where T(K) is the K-th triangular number.
!
!    Input, integer M, the degree of the polynomial.
!
!    Input, real ( kind = rk ) X(N), Y(N), the data locations.
!
!    Output, real ( kind = rk ) A(N,N), the Vandermonde matrix for X.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a(n,n)
  integer ex
  integer ey
  integer j
  integer m
  integer s
  integer tmp1
  integer triangle_num
  real ( kind = rk ) x(n)
  real ( kind = rk ) y(n)

  tmp1 = triangle_num ( m + 1 )

  if ( n /= tmp1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'VANDERMONDE_INTERP_2D_MATRIX - Fatal error!'
    write ( *, '(a)' ) '  For interpolation, we need N = T(M+1).'
    write ( *, '(a,i6)' ) '  But we have N = ', n
    write ( *, '(a,i6)' ) '  M = ', m
    write ( *, '(a,i6)' ) '  and T(M+1) = ', tmp1
    stop
  end if

  j = 0

  do s = 0, m
    do ex = s, 0, -1
      ey = s - ex
      j = j + 1
      a(1:n,j) = x(1:n) ** ex * y(1:n) ** ey
    end do
  end do
 
  return
end