function r_tiling = polyomino_multihedral_tiling_print ( r_shape, p_num, ...
  p_shapes, d, x, label );

%*****************************************************************************80
%
%% polyomino_multihedral_tiling_print prints a multihedral tiling.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 June 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer R_SHAPE(R_M,R_N), a matrix of 0's and 1's representing the 
%    region.
%
%    Input, integer P_NUM, the number of polyomino shapes to be used.
%
%    Input, integer P_SHAPES(M_R,N_R,P_NUM), a binary matrix describing the 
%    polyominoes.  The matrix is "top-left tight", that is, the I-th polyomino 
%    is stored in P_SHAPES(1:M_R,1:N_R,I) with a 1 in the first row, a 1 in the 
%    first column, but with trailing rows and columns of zeros allowed.
%
%    Input, integer D(P_NUM), tells us how many copies of each polyomino
%    we must use in the tiling.
%
%    Input, integer X(N), a solution to the problem, consisting of a binary 
%    value (0 or 1) for every possible placement of a polyomino variant onto R.
%
%    Input, string LABEL, a label for the printout, such as 
%    'Solution #3 of the 2x4 Example'.
%

%
%  Get the linear system satisfied by A*x=b.
%
  [ a, b, parent ] = polyomino_multihedral_matrix ( r_shape, p_num, p_shapes, d );
%
%  Get the size of the linear system.
%
  [ a_m, a_n ] = size ( a );
%
%  If X(J) is 1, replace it by J.  
%
  xv = (1:a_n)' .* x;
%
%  The first A_M - P_NUM - 1 rows of A are equations about covering 
%  each cell of the region.
%
%  Multiplying this matrix times XV gives us 
%
  axv = a(1:a_m-p_num-1,1:a_n) * xv;
%
%  Index the cells in R which are set to 1.
%
  r_index = polyomino_index ( r_shape );
%
%  We need to know the size of R.
%
  [ r_m, r_n ] = size ( r_shape );
%
%  R_SHAPE is binary (0 or 1).
%  R-TILING replaces each 1 by the index of the polyomino variant which
%  covers it.
%
  r_tiling = zeros ( r_m, r_n );

  for i = 1 : r_m
    for j = 1 : r_n
      if ( r_index(i,j) == 0 )
        r_tiling(i,j) = 0;
      else
        r_tiling(i,j) = axv(r_index(i,j));
      end
    end
    
  end
%
%  Print the tiling.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', label );
  fprintf ( 1, '\n' );

  smax = max ( max ( r_tiling ) );
  
  for i = 1 : r_m
    for j = 1 : r_n
      if ( 10 <= smax && smax < 100 )
        if ( r_tiling(i,j) == 0 )
          fprintf ( 1, '   ' );
        else
          fprintf ( 1, ' %2d', r_tiling(i,j) );
        end
      else
        if ( s(i,j) == 0 )
          fprintf ( 1, '  ' );
        else
          fprintf ( 1, ' %d', r_tiling(i,j) );
        end
      end
    end
    fprintf ( 1, '\n' );
  end

  return
end
