function output = simdo ( qind, qfor, iprod, kdim, input ) %*****************************************************************************80 % %% SIMDO generates multi-indices, simulating nested DO-loops. % % Discussion: % % The loops are assumed to be nested to a depth of K. % % The R-th loop is assumed to have upper limit N(R) and increment Inc(R). % % The total number of executions of the innermost loop is % % N = product ( 1 <= R <= K ) N(R). % % Let these executions be indexed by the single integer J, which % we call the index subscript. % % Each value of J corresponds to a particular set of loop indices, % which we call the subscript vector I(J). % % This routine can start with J and find I(J), or determine % J from I(J). % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 26 July 2008 % % Author: % % M O'Flaherty, G MacKenzie % % Reference: % % M O'Flaherty, G MacKenzie, % Algorithm AS 172: % Direct Simulation of Nested Fortran DO-LOOPS, % Applied Statistics, % Volume 31, Number 1, 1982, pages 71-74. % % Parameters: % % Input, logical QIND. % TRUE to convert an index subscript J to the subscript vector I(J). % FALSE to convert the subscript vector I(J) to the index subscript J. % % Input, logical QFOR, % TRUE if conversion is required in standard Fortran subscripting order, % FALSE otherwise. % % Input, integer IPROD(KDIM), contains the partial products. % If QFOR is FALSE, then % IPROD(S) = product ( 1 <= R <= S ) N(R). % If QFOR is TRUE, then % IPROD(S) = product ( 1 <= R <= S ) N(KDIM+1-R). % % Input, integer KDIM, the nesting depth of the loops. % % Input/output, integer JSUB. % If QIND is TRUE, then JSUB is an input quantity, an index subscript J % to be converted into the subscript vector I(J). % If QIND is FALSE, then JSUB is an output quantity, the index subscript J % corresponding to the subscript vector I(J). % % Input/output, integer IVEC(KDIM). % if QIND is TRUE, then IVEC is an output quantity, the subscript vector I(J) % corresponding to the index subscript J. % If QIND is FALSE, then IVEC is an input quantity, a subscript vector I(J) % for which the corresponding index subscript J is to be computed. % ifault = 0; % % Index subscript to subscript vector conversion. % if ( qind ) jsub = input; if ( iprod(kdim) < jsub ) ifault = 1; fprintf ( 1, '\n' ); fprintf ( 1, 'SIMDO - Fatal error!\n' ); fprintf ( 1, ' JSUB exceeds limit.\n' ); error ( 'SIMDO - Fatal error!' ); end itempv = jsub - 1; for i = 1 : kdim - 1 ik = kdim - i; ivec(i) = floor ( itempv / iprod(ik) ); itempv = itempv - iprod(ik) * ivec(i); ivec(i) = ivec(i) + 1; end ivec(kdim) = itempv + 1; if ( qfor ) ivec(1:kdim) = ivec(kdim:-1:1); end output = ivec; % % Subscript vector to index subscript conversion. % else ivec = input; if ( ~qfor ) ivec(1:kdim) = ivec(kdim:-1:1); end if ( iprod(1) < ivec(1) ) ifault = 2; fprintf ( 1, '\n' ); fprintf ( 1, 'SIMDO - Fatal error!\n' ); fprintf ( 1, ' Illegal entry in IVEC.\n' ); error ( 'SIMDO - Fatal error!' ); end for i = 2 : kdim if ( iprod(i) / iprod(i-1) < ivec(i) ) ifault = 2; fprintf ( 1, '\n' ); fprintf ( 1, 'SIMDO - Fatal error!\n' ); fprintf ( 1, ' Illegal entry in IVEC.\n' ); fprintf ( 1, ' IPROD(%d) = %d\n', i, iprod(i) ); fprintf ( 1, ' IPROD(%d) = %d\n', i-1, iprod(i-1) ); fprintf ( 1, ' IVEC(%d) = %d\n', i, ivec(i) ); error ( 'SIMDO - Fatal error!' ); end end jsub = ivec(1); for i = 2 : kdim jsub = jsub + ( ivec(i) - 1 ) * iprod(i-1); end % % As a courtesy to the caller, UNREVERSE the IVEC vector % if you reversed it. Doesn't matter in MATLAB since the % vector is input only. % % if ( ~qfor ) % ivec(1:kdim) = ivec(kdim:1:-1) % end output = jsub; end return end