function [find,vaf,targone,targtwo,targthree, ...
    outpermone,outpermtwo,outpermthree] = ...
    triarobfnd(prox,inperm,kblock)

% TRIAROBFND finds and fits the sum of three 
% anti-Robinson matrices using iterative projection to
% a symmetric proximity matrix in the $L_{2}$-norm 
% based on permutations identified through 
% the use of iterative quadratic assignment.
%
% syntax: [find,vaf,targone,targtwo,targthree,outpermone, ...
%     outpermtwo,outpermthree] = triarobfnd(prox,inperm,kblock)
%
% PROX is the input proximity matrix ($n \times n$ with a zero
% main diagonal and a dissimilarity interpretation);
% INPERM is a given starting permutation of the first $n$ integers;
% FIND is the least-squares optimal matrix (with 
% variance-accounted-for of VAF)
% to PROX and is the sum of the three anti-Robinson matrices
% TARGONE, TARGTWO, and TARGTHREE based on the three row and column
% object orderings given by the ending permutations OUTPERMONE, 
% OUTPERMTWO, and OUTPERMTHREE. KBLOCK defines the block size
% in the use of the iterative quadratic assignment routine.

n = size(prox,1);
proxave = zeros(n,n);

[targ1,vaftarg1,outperm1] = arobfnd(prox,inperm,kblock);

resprox1(outperm1,outperm1) = prox(outperm1,outperm1) - targ1;

[targ2,vaftarg2,outperm2] = arobfnd(resprox1,inperm,kblock);

resprox2(outperm2,outperm2) = resprox1(outperm2,outperm2) - targ2;

[targ3,vaftarg3,outperm3] = arobfnd(resprox2,inperm,kblock);

resprox3(outperm3,outperm3) = resprox2(outperm3,outperm3) - targ3;

find = prox - resprox3;

aveprox = sum(sum(prox))/(n*(n-1));

for i = 1:n
    for j = 1:n
        if(i ~= j)
            proxave(i,j) = aveprox;
        else
            proxave(i,j) = 0.0;
        end
    end
end

diff = sum(sum((prox - find).^2));
denom = sum(sum((prox - proxave).^2));
vaf = 1 - (diff/denom);
            
targone = targ1;
targtwo = targ2;
targthree = targ3;
outpermone = outperm1;
outpermtwo = outperm2;
outpermthree = outperm3;

vafdiff = 1.0;
iteration = 0;

while ((vafdiff >= 1.0e-006) && (iteration <= 100))
    
    vafprev = vaf;
    iteration = iteration + 1;
    
   
    resprox(outpermthree,outpermthree) = prox(outpermthree,outpermthree) - targthree;
    resprox(outpermtwo,outpermtwo) = resprox(outpermtwo,outpermtwo) - targtwo;
    [targone,vafone,outpermone] = arobfnd(resprox,outpermone,kblock);
    
    resprox(outpermthree,outpermthree) = prox(outpermthree,outpermthree) - targthree;
    resprox(outpermone,outpermone) = resprox(outpermone,outpermone) - targone;
    [targtwo,vaftwo,outpermtwo] = arobfnd(resprox,outpermtwo,kblock);
    
    resprox(outpermtwo,outpermtwo) = prox(outpermtwo,outpermtwo) - targtwo;
    resprox(outpermone,outpermone) = resprox(outpermone,outpermone) - targone;
    [targthree,vafthree,outpermthree] = arobfnd(resprox,outpermthree,kblock);
    
    
    find(outpermone,outpermone) = targone;
    find(outpermtwo,outpermtwo) = find(outpermtwo,outpermtwo) + targtwo;
    find(outpermthree,outpermthree) = find(outpermthree,outpermthree) + ...
        + targthree;
    
    diff = sum(sum((prox - find).^2));
    denom = sum(sum((prox - proxave).^2));
    vaf = 1 - (diff/denom);
    
    vafdiff = abs(vaf - vafprev);
end