function cost = pivot_sequence_to_matrix_chain_cost ( n_mats, p, dims )

%*****************************************************************************80
%
%% pivot_sequence_to_matrix_chain_cost() evaluates a particular matrix chain product.
%
%  Discussion:
%
%    There are n_mats matrices to multiply.
%    There are n_mats+1 dimensions to consider.
%    There are n_mats-1 matrix multiplications to carry out.
%
%    The cost, in terms of multiplications of pairs of real numbers,
%    of a multiplying a single pair of matrices A and B,
%    of dimensions LxM and MxN, is L*M*N.
%
%    The cost of multiplying a sequence of N matrices A*B*C*...*Z
%    will in general vary, depending on the order in which the N-1
%    multiplications are carried out.
%
%    This function assumes that a particular multiplication ordering 
%    has been specified in the pivot sequence p(), and returns the
%    corresponding cost.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 June 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer n_mats: the number of matrices to multiply.
%
%    integer p(n_mats-1): the pivot sequence.
%    On the i-th step of the multiplication procedure, there are
%    n-i multiplications to choose from, and we choose the p(i) one.
%    It must be the case that, for each i, 1 <= p(i) <= n - i
%
%    integer dim(n_mats+1): the matrix dimensions.
%    For 1 <= i < n, the i-th matrix has dimension dim(i) x dim(i+1).
%
%  Output:
%
%    integer cost: the matrix multiplication cost.
%
  n_dims = n_mats + 1;
  n_mults = n_mats - 1;

  cost = 0;
  n_dims2 = n_dims;
  dims2(1:n_dims) = dims(1:n_dims);
%
%  Carry out n_mults multiplications.
%  On step I, we carry out multiplication J.
%  Multiplication J involves 
%    dims(j)dims(j+1) * dims(j+1)dims(j+2) => dims(j)dims(j+2)
%  at an additional cost of dim(sj) * dims(j+1) * dims(j+2),
%  while removing dims(j+1) from the list of dimensions.
%
  for i = 1 : n_mults
    j = p(i);
    d1 = dims2(j);
    d2 = dims2(j+1);
    d3 = dims2(j+2);
    cost = cost + d1 * d2 * d3;
    dims2(j+1:n_dims2-1) = dims2(j+2:n_dims2);
    n_dims2 = n_dims2 - 1;
  end

  return
end
