function [inds,posmul]=ssp21sp(s,siz,fl)
% [inds,posmul]=SSP21SP(s,siz,fl)
%
% Expresses the distance between a linearly indexed grid node and all the
% other grid nodes as indices into the vector that contains the distances
% between the first linear index of the grid and all the other nodes
%
% INPUT:
%
% s The linear index into a two-dimensional array
% siz The size of the array (rows, columns)
% fl 0 Runs a non-redundant full loop with some extra overhead [default]
% 1 Runs a redundant full loop
%
% 2 Runs a short loop for all positions PAST the given one
%
% OUTPUT:
%
% inds The indices that map Dij(1,inds)=Dij(s,:) in the distance
% matrix Dij that contains all pairwise distances between the
% elements s and sp in the grid
% posmul The relevant index positions and their multiplicity
%
% SEE ALSO:
%
% SSPDIST, XXPDIST
%
% Last modified by fjsimons-at-alum.mit.edu, 01/15/2014
defval('fl',1)
% Keep time
t=tic;
% The size of the array
M=siz(1);
N=siz(2);
MN=M*N;
% A handful of special cases
if s==1
inds=1:MN;
posmul=[inds ; ones(size(inds))];
return
end
% Allocate index array
inds=zeros(1,MN);
%-[fl==3]*[s+1]);
if s==MN
inds=MN:-1:1;
posmul=[inds ; ones(size(inds))];
return
end
% Only one remaining case after contraction of loops below is when max/min
% don't commute which is when there is a NaN from 0/0, on the last grid row,
% i.e. when a=-1 and d=0 it fails to fill out to the left since we picked
% the max-min order to favor the right side... so prepare for that
if rem(s-1,M)+1==M
inds(s-M:-M:2)=[1:round([s-M-1]/M)]*M+1;
end
if s>1 && s0 && 1==3
disp(sprintf('%i out of %i are repeated hits, time taken %6.4fs',...
redfax,MN,toc(t)))
end
end
% Check that every element has been hit as it should
diferm(sum(~~inds),MN)
% If second output is requested
if nargout==2
pos=unique(inds);
posmul=hist(inds,pos);
% And the reconstruction must be happening
diferm(gamini(pos,posmul),sort(inds))
% Return both the position and its multiplicity
posmul=[pos; posmul];
end