function cuda_loop ( blocks, threads, n )

%*****************************************************************************80
%
%% CUDA_LOOP simulates the behavior of a CUDA loop.
%
%  Discussion:
%
%    A CUDA kernel "kernel()" is invoked by a command of the form
%     
%      kernel << blocks, threads >> ( args )
%
%    where blocks and threads are each vectors of up to 3 values,
%    listing the number of blocks and number of threads to be used.
%
%    If a problem involves N tasks, then tasks are allotted to 
%    specific CUDA processes in an organized fashion.  Some processes
%    may get no tasks, one task, or multiple tasks.  
%
%    Each process is given variables that can be used to determine
%    the tasks to be performed:
%
%      gridDim.x, gridDim.y, gridDim.z: the block dimensions as
%      given by the user in "blocks";
%
%      blockDim.x, blockDim.y, blockDim.z: the thread dimensions as
%      given by the user in "threads";
%
%      blockIdx.x, blockIdx.y, blockId.z: the block indices for this process.
%
%      threadIdx.x, threadIdx.y, threadIdx.z: the thread indices for this process.
%
%    Essentially, a process can determine its linear index K by:
%
%      K = threadIdx.x
%        +  blockdim.x  * threadIdx.y
%        +  blockDim.x  *  blockDim.y  * threadIdx.z
%        +  blockDim.x  *  blockDim.y  *  blockDim.z  * blockIdx.x
%        +  blockDim.x  *  blockDim.y  *  blockDim.z  *  gridDim.x  * blockIdx.y
%        +  blockDim.x  *  blockDim.y  *  blockDim.z  *  gridDim.x  *  gridDim.y  * blockIdx.z
%
%    Set task T = K.
%
%    while ( T < N )
%      carry out task T;
%      T = T + blockDim.x * blockDim.y * blockDim.z * gridDim.x * gridDim.y * gridDim.z.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    29 March 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer BLOCKS[3], the CUDA block values.  These 
%    should be nonnegative.  Typically, the third entry is 1.  Generally, 
%    the first two values cannot be greater than 35,535.
%
%    Input, integer THREADS[3], the CUDA thread values.  These 
%    should be nonnegative.  Typically, there is a maximum value imposed on 
%    these quantities, which depends on the GPU model.
%
%    Input, integer N, the number of tasks to be carried out.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CUDA_LOOP:\n' );
  fprintf ( 1, '  Simulate the assignment of N tasks to the blocks\n' );
  fprintf ( 1, '  and threads of a GPU using CUDA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of tasks is %d\n', n );
  fprintf ( 1, '  BLOCKS:  {%4d,%4d,%4d}\n', blocks(1), blocks(2), blocks(3) );
  fprintf ( 1, '  THREADS: {%4d,%4d,%4d}\n', threads(1), threads(2), threads(3) );

  k1 = 0;

  blockDimx = threads(1);
  blockDimy = threads(2);
  blockDimz = threads(3);

  gridDimx = blocks(1);
  gridDimy = blocks(2);
  gridDimz = blocks(3);

  chunk = blocks(2) * blocks(1) * threads(3) * threads(2) * threads(1);
  fprintf ( 1, '  Total threads = %d\n', chunk );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Process   Process (bx,by,bz) (tx,ty,tz)  Tasks...\n' );
  fprintf ( 1, '  Increment Formula\n' );
  fprintf ( 1, '\n' );

  for blockIdz = 0 : gridDimz - 1
    for blockIdy = 0 : gridDimy - 1
      for blockIdx = 0 : gridDimx - 1
        for threadIdz = 0 : blockDimz - 1
          for threadIdy = 0 : blockDimy - 1
            for threadIdx = 0 : blockDimx - 1
              t = k1;
              k2 = ...
                threadIdx ...
                + blockDimx * threadIdy ...
                + blockDimx * blockDimy * threadIdz ...
                + blockDimx * blockDimy * blockDimz * blockIdx ...
                + blockDimx * blockDimy * blockDimz * gridDimx * blockIdy ...
                + blockDimx * blockDimy * blockDimz * gridDimx * gridDimy * blockIdz;

              fprintf ( 1, '  %7d  %7d: (%2d,%2d,%2d) (%2d,%2d,%2d', ...
                k1, k2, blockIdx, blockIdy, blockIdz, threadIdx, threadIdy, threadIdz );
              while ( t < n )
                fprintf ( 1, '%3d', t );
                t = t + chunk;
              end
              fprintf ( 1, '\n' );
              k1 = k1 + 1;
            end
          end
        end
      end
    end
  end

  return
end

