function colList=colListing(m,n)
mn=m*n;
tm=0.5*m*(m+1);
tm1=0.5*m*(m-1);
tn=0.5*n*(n+1);
tn1=0.5*n*(n-1);
tmn=0.5*mn*(mn+1);
tmn1=0.5*mn*(mn-1); % t(mn-1)

%% part 1: off diagonal entries
jcol=ceil( -0.5+0.5*sqrt(1+8*(1:tmn1)) ) + 1;
icol=(1:tmn1) - 0.5*(jcol-1).*(jcol-2);

kcol=ceil(icol/m);
lcol=ceil(jcol/m);

iicol=icol-m*(kcol-1);
jjcol=jcol-m*(lcol-1);

mtxSmall=[kcol;lcol;iicol;jjcol];
swapkl=logical(iicol>jjcol);
mtxSmall(1:2,swapkl)=[lcol(swapkl); kcol(swapkl)];
mtxSmall(3:4,swapkl)=[jjcol(swapkl); iicol(swapkl)];

% for block indexing
BlkOffDiagLe=logical(kcol<lcol);
BlkOffDiagGe=logical(kcol>lcol);
BlkDiag=logical(kcol==lcol);
mtxSmall1rowOffDiagLe=(n-1)*(mtxSmall(2,BlkOffDiagLe)-1)+...
    mtxSmall(1,BlkOffDiagLe);
mtxSmall1rowOffDiagGe=(n-1)*(mtxSmall(2,BlkOffDiagGe)-1)+...
    mtxSmall(1,BlkOffDiagGe)-1;
mtxSmall1rowDiag=n*(n-1)+mtxSmall(1,BlkDiag);

% for entries within blocks
inBlkOffDiag=logical(iicol~=jjcol);
mtxSmall2rowOffDiag=mtxSmall(3,inBlkOffDiag)+...
    0.5.*(mtxSmall(4,inBlkOffDiag)-2).*(mtxSmall(4,inBlkOffDiag)-1);
mtxSmall2rowDiag=mtxSmall(3,~inBlkOffDiag)+tn1;

% temporary checking
mtxSmall_combined=zeros(2,tmn1);
mtxSmall_combined(1,BlkOffDiagLe)=mtxSmall1rowOffDiagLe;% 5th row for block indexing
mtxSmall_combined(1,BlkOffDiagGe)=mtxSmall1rowOffDiagGe;
mtxSmall_combined(1,BlkDiag)=mtxSmall1rowDiag;
mtxSmall_combined(2,inBlkOffDiag)=mtxSmall2rowOffDiag; % 6th row for in-block indexing
mtxSmall_combined(2,~inBlkOffDiag)=mtxSmall2rowDiag;

%% part 2: diagonal entries
%diagSmall=[kron(1:n,ones(2,4)); repmat(1:m,2,n)];
diagSmall1=kron((n*(n-1)+1):(n*(n-1)+m),ones(1,n));     % blocks
diagSmall2=repmat((tm1+1):tm,1,n);
diagSmall_combined=[diagSmall1; diagSmall2];

%% part3: put everything together
mtx=[mtxSmall_combined , mtxSmall_combined , diagSmall_combined; ...
     ones(1,tmn1), 2*ones(1,tmn1), ones(1,mn)];
 
colList=zeros(1,mn^2);
counter=0;

% off diagonal entries first
for ii=1:tm1,   % the ii-th entry
    temp=find(mtx(2,:)==ii);
    if length(temp)~=2*n^2,
        disp('WARNING: off diagonal, length(find(mtx(2,:)==ii))~=2*n^2??');
        keyboard
    end
    colList((counter+1):2:(counter+2*n^2))=temp(1:n^2);
    colList((counter+2):2:(counter+2*n^2))=temp(n^2+1:end);
    counter=counter+2*n^2;
end

% then the diagonal entries
for ii=tm1+1:tm
    temp=find(mtx(2,:)==ii);
    if length(temp)~=n^2,
        disp('WARNING: on diagonal, length(find(mtx(2,:)==ii))~=n^2??');
        keyboard
    end
    colList((counter+1):2:(counter+n*(n-1)))=temp(1:tn1);
    colList((counter+2):2:(counter+n*(n-1)))=temp(1+tn1:n*(n-1));
    colList((counter+n*(n-1)+1):counter+n^2)=temp((n*(n-1)+1):n^2);
    counter=counter+n^2;
end
 