function [outperm, rawindex, allperms, index, squareprox] = ...
   ordertm(proxtm, targ, inperm, kblock)

% ORDERTM carries out an iterative 
% Quadratic Assignment maximization task using the
% two-mode proximity matrix PROXTM 
% (with entries deviated from the mean proximity) 
% in the upper-right- and lower-left-hand portions of
% a defined square ($n x n$) proximity matrix 
% (called SQUAREPROX with a dissimilarity interpretation)
% with zeros placed elsewhere ($n$ = number of rows + 
% number of columns of PROXTM = nrow + ncol).
%
% syntax: [outperm, rawindex, allperms, index, squareprox] = ...
%   ordertm(proxtm, targ, inperm, kblock)
%
% Three separate local operations are used to permute
% the rows and columns of the square 
% proximity matrix to maximize the cross-product
% index with respect to a square target matrix TARG:
% pairwise interchanges of objects in the 
% permutation defining the row and column 
% order of the square proximity matrix; the insertion of from 1 to 
% KBLOCK (which is less than or equal to $n-1$) consecutive objects
% in the permutation defining the row and column order of the 
% data matrix; the rotation of from 2 to KBLOCK (which is less than
% or equal to $n-1$) consecutive objects in 
% the permutation defining the row and column order of the data 
% matrix. INPERM is the input beginning permutation (a permutation
% of the first $n$ integers).
% PROXTM is the two-mode $nrow x ncol$ input proximity matrix.
% TARG is the $n x n$ input target matrix.
% OUTPERM is the final permutation of SQUAREPROX with the 
% cross-product index RAWINDEX
% with respect to TARG. ALLPERMS is a cell array containing INDEX 
% entries corresponding to all the
% permutations identified in the optimization from ALLPERMS{1} 
% = INPERM to ALLPERMS{INDEX} = OUTPERM.

[nrow ncol] = size(proxtm);
n = nrow + ncol;
prox = zeros(n,n);
aveprox = sum(sum(proxtm))/(nrow*ncol);

for irow = 1:nrow
   for icol = 1:ncol
      
      prox(irow,icol+nrow) = proxtm(irow,icol) - aveprox;
      prox(icol+nrow,irow) = proxtm(irow,icol) - aveprox;
      
   end
end

outperm = inperm;
index = 1;
allperms{index} = inperm;
begindex = sum(sum(prox(inperm,inperm).*targ));

for iterate = 1:2
   
   nchange = 1;
   
while (nchange == 1)
   
   nchange=0;
   
   for k = 1:(n-1)
      for j = (k+1):n
                   
         intrperm = outperm;
             
         intrperm(k) = outperm(j);
         intrperm(j) = outperm(k);
         
         tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
         if(tryindex > (begindex + 1.0e-008))
            nchange = 1;
            begindex = tryindex;
            outperm = intrperm;
            index = index + 1;
            allperms{index} = intrperm;
         end
         
      end
   end
end

rawindex = begindex;
nchange = 1;

while (nchange == 1)
   
   nchange=0;
   
   for k = 1:kblock
      for insertpt = 1:(n+1)
         for nlimlow = 1:(n+1-k)
            
            intrperm = outperm;
            
            if (nlimlow > insertpt)
               
               jtwo = 0;
               for j = insertpt:(insertpt+k-1)
                  intrperm(j) =outperm(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = (insertpt+k):(nlimlow+k-1);
                  intrperm(j) = outperm(insertpt+jone);
                  jone = jone + 1;
               end
               
            elseif ((nlimlow+k) < insertpt)
               
               jtwo = 0;
               for j = (insertpt-k):(insertpt-1)
                  intrperm(j) = outperm(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = nlimlow:(insertpt-k-1)
                  intrperm(j) = outperm(nlimlow+k+jone);
                  jone = jone + 1;
               end
               
            else
               
            end
            
            tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
            if(tryindex > (begindex + 1.0e-008))
               nchange = 1;
               begindex = tryindex;
               outperm = intrperm;
               index = index +1;
               allperms{index} = intrperm;
            end
            
         end
      end
   end
end


rawindex = begindex;
nchange = 1;



while (nchange == 1)
   
   nchange=0;
   
   for k = 2:kblock
      for nlimlow = 1:(n+1-k)
                   
         intrperm = outperm;
            
         for j = 1:k
            intrperm(nlimlow+j-1) = outperm(nlimlow+k-j);
         end
         
         tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
         if(tryindex > (begindex + 1.0e-008))
            nchange = 1;
            begindex = tryindex;
            outperm = intrperm;
            index = index + 1;
            allperms{index} = intrperm;
         end
         
      end
   end
end



rawindex = begindex;
nchange = 1;


end


squareprox = prox;