module mexoperators
 interface mxs
    module procedure mxs_r2
    module procedure mxs_c2
    module procedure mxs_i2
    module procedure mxs_r1
    module procedure mxs_c1
    module procedure mxs_i1
    module procedure mxs_r
    module procedure mxs_c
    module procedure mxs_i
 end interface mxs
 interface assignment (=)
    module procedure l_to_i
    module procedure l1_to_i1
    module procedure l2_to_i2
    module procedure l_to_r
    module procedure l1_to_r1
    module procedure l2_to_r2
    module procedure l_to_c
    module procedure l1_to_c1
    module procedure l2_to_c2
 end interface
contains

 function mxs_r2(a) result(out)
  real(8), dimension(:,:) :: a
  real(8) out
  out=a(1,1)
 end function mxs_r2
 function mxs_c2(a) result(out)
  complex(8), dimension(:,:) :: a
  complex(8) out
  out=a(1,1)
 end function mxs_c2
 function mxs_i2(a) result(out)
  integer, dimension(:,:) :: a
  integer out
  out=a(1,1)
 end function mxs_i2
 function mxs_r1(a) result(out)
  real(8), dimension(:) :: a
  real(8) out
  out=a(1)
 end function mxs_r1
 function mxs_c1(a) result(out)
  complex(8), dimension(:) :: a
  complex(8) out
  out=a(1)
 end function mxs_c1
 function mxs_i1(a) result(out)
  integer, dimension(:) :: a
  integer out
  out=a(1)
 end function mxs_i1
 function mxs_r(a) result(out)
  real(8) :: a
  real(8) out
  out=a
 end function mxs_r
 function mxs_c(a) result(out)
  complex(8) :: a
  complex(8) out
  out=a
 end function mxs_c
 function mxs_i(a) result(out)
  integer :: a
  integer out
  out=a
 end function mxs_i

 subroutine l_to_i (numout, login)
  logical, intent(in)  :: login
  integer, intent(out) :: numout
  if (login) then
   numout=1.0
  else
   numout=0.0
  end if
 end subroutine l_to_i
 subroutine l1_to_i1 (numout, login)
  logical, intent(in)  :: login(:)
  integer, intent(out) :: numout(size(login))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l1_to_i1
 subroutine l2_to_i2 (numout, login)
  logical, intent(in)  :: login(:,:)
  integer, intent(out) :: numout(size(login,1),size(login,2))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l2_to_i2
 subroutine l_to_r (numout, login)
  logical, intent(in)  :: login
  real(8), intent(out) :: numout
  if (login) then
   numout=1.0
  else
   numout=0.0
  end if
 end subroutine l_to_r
 subroutine l1_to_r1 (numout, login)
  logical, intent(in)  :: login(:)
  real(8), intent(out) :: numout(size(login))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l1_to_r1
 subroutine l2_to_r2 (numout, login)
  logical, intent(in)  :: login(:,:)
  real(8), intent(out) :: numout(size(login,1),size(login,2))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l2_to_r2
 subroutine l_to_c (numout, login)
  logical, intent(in)  :: login
  complex(8), intent(out) :: numout
  if (login) then
   numout=1.0
  else
   numout=0.0
  end if
 end subroutine l_to_c
 subroutine l1_to_c1 (numout, login)
  logical, intent(in)  :: login(:)
  complex(8), intent(out) :: numout(size(login))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l1_to_c1
 subroutine l2_to_c2 (numout, login)
  logical, intent(in)  :: login(:,:)
  complex(8), intent(out) :: numout(size(login,1),size(login,2))
  where (login) numout=1.0
  where (.not.login) numout=0.0
 end subroutine l2_to_c2

end module mexoperators



module mexcallback
 interface mlcall
  module procedure callback_r
 end interface mlcall

 interface mlcall0
  module procedure callback0_r
 end interface mlcall0
contains
 subroutine callback_r(out,in_1,funstr)
 real(8), dimension(:,:) :: in_1
 integer :: in_1_ptr
 real(8), pointer, dimension(:,:) :: out
 integer :: out_ptr
 character (len=*) :: funstr
 integer lhs(50), rhs(50)
 integer out_m, out_n, in_1_m, in_1_n
  in_1_m=size(in_1,dim=1); in_1_n=size(in_1,dim=2)
  rhs(1)=mxCreateFull(in_1_m,in_1_n,0)
  in_1_ptr=mxGetPr(rhs(1))
  call mxCopyReal8ToPtr(in_1,in_1_ptr,in_1_m*in_1_n)
  call mxSetM(rhs(1),in_1_m);call mxSetN(rhs(1),in_1_n)
  call mexCallMATLAB(1,lhs(1),1,rhs(1:1),funstr)
  out_m=mxGetM(lhs(1));out_n=mxGetN(lhs(1))
  allocate(out(out_m,out_n))
  out_ptr=mxGetPr(lhs(1))
  call mxCopyPtrToReal8(out_ptr,out,out_m*out_n)
 end subroutine callback_r

 subroutine callback0_r(in_1,funstr)
 real(8), dimension(:,:) :: in_1
 integer :: in_1_ptr
 character (len=*) :: funstr
 integer lhs(50), rhs(50)
 integer in_1_m, in_1_n
  in_1_m=size(in_1,dim=1); in_1_n=size(in_1,dim=2)
  rhs(1)=mxCreateFull(in_1_m,in_1_n,0)
  in_1_ptr=mxGetPr(rhs(1))
  call mxCopyReal8ToPtr(in_1,in_1_ptr,in_1_m*in_1_n)
  call mxSetM(rhs(1),in_1_m);call mxSetN(rhs(1),in_1_n)
  call mexCallMATLAB(0,NULL,1,rhs(1:1),funstr)
 end subroutine callback0_r

end module mexcallback



module mexfunctions
contains
subroutine placeholder(m)
integer m
end subroutine placeholder


end module mexfunctions



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Gateway!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mexfunction(nlhs, plhs, nrhs, prhs)
integer(4) plhs(*), prhs(*)
integer(4) nlhs, nrhs      
! input and output pointers
integer(4) :: prox_ptr,inperm_ptr,dumone_ptr
real(8), allocatable :: prox(:,:),inperm(:,:),dumone(:,:)
integer(4) :: find_ptr,vaf_ptr
real(8), allocatable :: find(:,:),vaf(:,:)
!     Any other variables needed
integer prox_m,prox_n,inperm_m,inperm_n,dumone_m,dumone_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 3) then
  print *,'ms_atreefnd requires 3 input arguments';return
elseif (nlhs .ne. 2) then
  print *,'ms_atreefnd requires 2 output arguments';return
endif
!     Get the sizes of all the input variables
prox_m=mxGetm(prhs(1));prox_n=mxGetn(prhs(1))
inperm_m=mxGetm(prhs(2));inperm_n=mxGetn(prhs(2))
dumone_m=mxGetm(prhs(3));dumone_n=mxGetn(prhs(3))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(prox_m,prox_n,0)
plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
find_ptr=mxGetPr(plhs(1))
vaf_ptr=mxGetPr(plhs(2))
!     Copy right hand arguments to local arrays
prox_ptr=mxGetPr(prhs(1))
inperm_ptr=mxGetPr(prhs(2))
dumone_ptr=mxGetPr(prhs(3))
!     Allocate and copy data to arrays
allocate(prox(prox_m,prox_n),inperm(inperm_m,inperm_n),dumone(dumone_m,dumone_n))
allocate(find(prox_m,prox_n),vaf(dumone_m,dumone_n))
call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n)
call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n)
call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n)
!     Do the actual computations in a subroutine
call ms_atreefnd(find,vaf(1,1),prox,inperm,dumone(1,1))
call mxCopyReal8ToPtr(find,find_ptr,prox_m*prox_n)
call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n)
deallocate(find,vaf,prox,inperm,dumone)
return

contains


subroutine ms_atreefnd(find,vaf,prox,inperm,dumone)
!     COMPUTATIONAL SUBROUTINE
!     Before anything, list what modules we will use.
use mexfunctions;   use mexoperators;   use mexcallback
!     size variables
!     Input/Output local mirrors
real(8) prox(:,:)
real(8) inperm(:,:)
real(8) dumone
real(8) find(:,:)
real(8) vaf
!     Matlab function pointers
!     All other local variables
real(8) denom
real(8) diff
real(8) aveprox
real(8), allocatable :: fit(:,:)
real(8), allocatable :: targ(:,:)
real(8) del
real(8) p6
real(8) p5
real(8) p4
real(8) p3
real(8) p2
real(8) p1
integer jfour
integer jthree
integer jtwo
integer jone
real(8) temp
integer j
integer i
integer ifour
integer ithree
integer itwo
integer ione
real(8) indexll
real(8) indxord(1,4)
real(8), allocatable :: proxave(:,:)
real(8) iterate
real(8) cr
real(8) work(657720,1)
real(8) n
!     Fill in vars going in and out
allocate( proxave(size(prox,1),size(prox,2)), targ(size(prox,1),size(prox,2)), fit(s&
      &ize(prox,1),size(prox,2)))
! --- Main computational routine. ---------------------------------!
n = size(prox,1)
work = 0.
find = prox
cr = 1.0
iterate = 0.0
proxave = 0.
indxord = 0.
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
iterate = iterate + 1.0
do ione = 1,nint((n-3))
do itwo = (ione+1),nint((n-2))
do ithree = (itwo+1),nint((n-1))
do ifour = (ithree+1),nint(n)
indxord(1,1) = inperm(1,ione)
indxord(1,2) = inperm(1,itwo)
indxord(1,3) = inperm(1,ithree)
indxord(1,4) = inperm(1,ifour)
do i = 1,3
do j = (i+1),4
if((indxord(1,i) > indxord(1,j))) then
temp = indxord(1,i)
indxord(1,i) = indxord(1,j)
indxord(1,j) = temp
endif
enddo
enddo
jone = indxord(1,1)
jtwo = indxord(1,2)
jthree = indxord(1,3)
jfour = indxord(1,4)
p1 = find(jone,jtwo)
p2 = find(jone,jthree)
p3 = find(jone,jfour)
p4 = find(jtwo,jthree)
p5 = find(jtwo,jfour)
p6 = find(jthree,jfour)
if((iterate <= 100.0)) then
find(jone,jtwo) = find(jone,jtwo) - work(int(indexll)+1,1)
find(jone,jthree) = find(jone,jthree) - work(int(indexll)+2,1)
find(jone,jfour) = find(jone,jfour) - work(int(indexll)+3,1)
find(jtwo,jthree) = find(jtwo,jthree) - work(int(indexll)+4,1)
find(jtwo,jfour) = find(jtwo,jfour) - work(int(indexll)+5,1)
find(jthree,jfour) = find(jthree,jfour) - work(int(indexll)+6,1)
endif
if((((find(jone,jtwo) + find(jthree,jfour)) <=                   (find(jone,jthree) &
      &+ find(jtwo,jfour))) .and.                   ((find(jone,jtwo) + find(jthree,&
      &jfour)) <=                    (find(jone,jfour) + find(jtwo,jthree))))) then
del = (find(jone,jthree) + find(jtwo,jfour) -                   (find(jone,jfour) + &
      &find(jtwo,jthree)))/4.0
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = -del
work(int(indexll)+6,1) = 0.0
find(jone,jthree) = find(jone,jthree) - del
find(jtwo,jfour) = find(jtwo,jfour) - del
find(jone,jfour) = find(jone,jfour) + del
find(jtwo,jthree) = find(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((find(jone,jthree) + find(jtwo,jfour)) <=                   (find(jone,jtwo) +&
      & find(jthree,jfour))) .and.                   ((find(jone,jthree) + find(jtwo&
      &,jfour)) <=                    (find(jone,jfour) + find(jtwo,jthree))))) then
del = (find(jone,jtwo) + find(jthree,jfour) -                   (find(jone,jfour) + &
      &find(jtwo,jthree)))/4.0
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = 0.0
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = 0.0
work(int(indexll)+6,1) = -del
find(jone,jtwo) = find(jone,jtwo) - del
find(jthree,jfour) = find(jthree,jfour) - del
find(jone,jfour) = find(jone,jfour) + del
find(jtwo,jthree) = find(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((find(jone,jfour) + find(jtwo,jthree)) <=                   (find(jone,jtwo) +&
      & find(jthree,jfour))) .and.                   ((find(jone,jfour) + find(jtwo,&
      &jthree)) <=                    (find(jone,jthree) + find(jtwo,jfour))))) then
del = (find(jone,jtwo) + find(jthree,jfour) -                   (find(jone,jthree) +&
      & find(jtwo,jfour)))/4.0
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) =  del
work(int(indexll)+3,1) =  0.0
work(int(indexll)+4,1) =  0.0
work(int(indexll)+5,1) =  del
work(int(indexll)+6,1) = -del
find(jone,jtwo) = find(jone,jtwo) - del
find(jthree,jfour) = find(jthree,jfour) - del
find(jone,jthree) = find(jone,jthree) + del
find(jtwo,jfour) = find(jtwo,jfour) + del
indexll = indexll + 6.0
endif
cr = cr + abs(p1-find(jone,jtwo)) + abs(p2-find(jone,jthree))                   + ab&
      &s(p3-find(jone,jfour)) + abs(p4-find(jtwo,jthree))                   + abs(p5&
      &-find(jtwo,jfour)) + abs(p6-find(jthree,jfour))
enddo
enddo
enddo
enddo
enddo
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
find(jtwo,jone) = find(jone,jtwo)
enddo
enddo
targ = find
n = size(prox,1)
work = 0.
fit = prox
cr = 1.0
proxave = 0.
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
do jone = 1,nint((n-3))
do jtwo = (jone+1),nint((n-2))
do jthree = (jtwo+1),nint((n-1))
do jfour = (jthree+1),nint(n)
p1 = fit(jone,jtwo)
p2 = fit(jone,jthree)
p3 = fit(jone,jfour)
p4 = fit(jtwo,jthree)
p5 = fit(jtwo,jfour)
p6 = fit(jthree,jfour)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1)
fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1)
fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1)
fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1)
fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1)
if ((((targ(jone,jtwo) + targ(jthree,jfour)) <=                   (targ(jone,jthree)&
      & + targ(jtwo,jfour))) .and.                   ((targ(jone,jtwo) + targ(jthree&
      &,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
del = (fit(jone,jthree) + fit(jtwo,jfour) -                   (fit(jone,jfour) + fit&
      &(jtwo,jthree)))/4.0
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = -del
work(int(indexll)+6,1) = 0.0
fit(jone,jthree) = fit(jone,jthree) - del
fit(jtwo,jfour) = fit(jtwo,jfour) - del
fit(jone,jfour) = fit(jone,jfour) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jthree) + targ(jtwo,jfour)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                   ((targ(jone,jthree) + targ(jtwo&
      &,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
del = (fit(jone,jtwo) + fit(jthree,jfour) -                   (fit(jone,jfour) + fit&
      &(jtwo,jthree)))/4.0
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = 0.0
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = 0.0
work(int(indexll)+6,1) = -del
fit(jone,jtwo) = fit(jone,jtwo) - del
fit(jthree,jfour) = fit(jthree,jfour) - del
fit(jone,jfour) = fit(jone,jfour) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jfour) + targ(jtwo,jthree)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                   ((targ(jone,jfour) + targ(jtwo,&
      &jthree)) <=                    (targ(jone,jthree) + targ(jtwo,jfour))))) then
del = (fit(jone,jtwo) + fit(jthree,jfour) -                   (fit(jone,jthree) + fi&
      &t(jtwo,jfour)))/4.0
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) =  del
work(int(indexll)+3,1) =  0.0
work(int(indexll)+4,1) =  0.0
work(int(indexll)+5,1) =  del
work(int(indexll)+6,1) = -del
fit(jone,jtwo) = fit(jone,jtwo) - del
fit(jthree,jfour) = fit(jthree,jfour) - del
fit(jone,jthree) = fit(jone,jthree) + del
fit(jtwo,jfour) = fit(jtwo,jfour) + del
indexll = indexll + 6.0
endif
cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree))                   + abs(&
      &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree))                   + abs(p5-fit&
      &(jtwo,jfour)) + abs(p6-fit(jthree,jfour))
enddo
enddo
enddo
enddo
do jone = 1,nint((n-3))
do jtwo = (jone+1),nint((n-2))
do jthree = (jtwo+1),nint((n-1))
do jfour = (jthree+1),nint(n)
p1 = fit(jone,jtwo)
p2 = fit(jone,jthree)
p3 = fit(jone,jfour)
p4 = fit(jtwo,jthree)
p5 = fit(jtwo,jfour)
p6 = fit(jthree,jfour)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1)
fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1)
fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1)
fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1)
fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1)
if ((((targ(jone,jtwo) + targ(jthree,jfour)) <=                   (targ(jone,jthree)&
      & + targ(jtwo,jfour))) .and.                   ((targ(jone,jtwo) + targ(jthree&
      &,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
if (((fit(jone,jtwo) + fit(jthree,jfour)) <=                   (fit(jone,jthree) + f&
      &it(jtwo,jfour)))) then
del = 0.0
else
del = (fit(jone,jtwo) + fit(jthree,jfour) -                       (fit(jone,jthree) &
      &+ fit(jtwo,jfour)))/4.0
endif
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = del
work(int(indexll)+3,1) =  0.0
work(int(indexll)+4,1) = 0.0
work(int(indexll)+5,1) = del
work(int(indexll)+6,1) = -del
fit(jone,jtwo) = fit(jone,jtwo) - del
fit(jthree,jfour) = fit(jthree,jfour) - del
fit(jone,jthree) = fit(jone,jthree) + del
fit(jtwo,jfour) = fit(jtwo,jfour) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jthree) + targ(jtwo,jfour)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                    ((targ(jone,jthree) + targ(jtw&
      &o,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
if (((fit(jone,jthree) + fit(jtwo,jfour)) <=                   (fit(jone,jtwo) + fit&
      &(jthree,jfour)))) then
del = 0.0
else
del = (fit(jone,jthree) + fit(jtwo,jfour) -                       (fit(jone,jtwo) + &
      &fit(jthree,jfour)))/4.0
endif
work(int(indexll)+1,1) = del
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) =  0.0
work(int(indexll)+4,1) = 0.0
work(int(indexll)+5,1) = -del
work(int(indexll)+6,1) = del
fit(jone,jthree) = fit(jone,jthree) - del
fit(jtwo,jfour) = fit(jtwo,jfour) - del
fit(jone,jtwo) = fit(jone,jtwo) + del
fit(jthree,jfour) = fit(jthree,jfour) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jfour) + targ(jtwo,jthree)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                   ((targ(jone,jfour) + targ(jtwo,&
      &jthree)) <=                    (targ(jone,jthree) + targ(jtwo,jfour))))) then
if (((fit(jone,jfour) + fit(jtwo,jthree)) <=                   (fit(jone,jtwo) + fit&
      &(jthree,jfour)))) then
del = 0.0
else
del = (fit(jone,jfour) + fit(jtwo,jthree) -                       (fit(jone,jtwo) + &
      &fit(jthree,jfour)))/4.0
endif
work(int(indexll)+1,1) = del
work(int(indexll)+2,1) = 0.0
work(int(indexll)+3,1) =  -del
work(int(indexll)+4,1) = -del
work(int(indexll)+5,1) = 0.0
work(int(indexll)+6,1) = del
fit(jone,jfour) = fit(jone,jfour) - del
fit(jtwo,jthree) = fit(jtwo,jthree) - del
fit(jone,jtwo) = fit(jone,jtwo) + del
fit(jthree,jfour) = fit(jthree,jfour) + del
indexll = indexll + 6.0
endif
cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree))                   + abs(&
      &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree))                   + abs(p5-fit&
      &(jtwo,jfour)) + abs(p6-fit(jthree,jfour))
enddo
enddo
enddo
enddo
do jone = 1,nint((n-3))
do jtwo = (jone+1),nint((n-2))
do jthree = (jtwo+1),nint((n-1))
do jfour = (jthree+1),nint(n)
p1 = fit(jone,jtwo)
p2 = fit(jone,jthree)
p3 = fit(jone,jfour)
p4 = fit(jtwo,jthree)
p5 = fit(jtwo,jfour)
p6 = fit(jthree,jfour)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1)
fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1)
fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1)
fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1)
fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1)
if ((((targ(jone,jtwo) + targ(jthree,jfour)) <=                   (targ(jone,jthree)&
      & + targ(jtwo,jfour))) .and.                   ((targ(jone,jtwo) + targ(jthree&
      &,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
if(((fit(jone,jtwo) + fit(jthree,jfour)) <=                   (fit(jone,jfour) + fit&
      &(jtwo,jthree)))) then
del = 0.0
else
del = (fit(jone,jtwo) + fit(jthree,jfour) -                       (fit(jone,jfour) +&
      & fit(jtwo,jthree)))/4.0
endif
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = 0.0
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = 0.0
work(int(indexll)+6,1) = -del
fit(jone,jtwo) = fit(jone,jtwo) - del
fit(jthree,jfour) = fit(jthree,jfour) - del
fit(jone,jfour) = fit(jone,jfour) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jthree) + targ(jtwo,jfour)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                   ((targ(jone,jthree) + targ(jtwo&
      &,jfour)) <=                    (targ(jone,jfour) + targ(jtwo,jthree))))) then
if(((fit(jone,jthree) + fit(jtwo,jfour)) <=                   (fit(jone,jfour) + fit&
      &(jtwo,jthree)))) then
del = 0.0
else
del = (fit(jone,jthree) + fit(jtwo,jfour) -                       (fit(jone,jfour) +&
      & fit(jtwo,jthree)))/4.0
endif
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) =  del
work(int(indexll)+4,1) = del
work(int(indexll)+5,1) = -del
work(int(indexll)+6,1) = 0.0
fit(jone,jthree) = fit(jone,jthree) - del
fit(jtwo,jfour) = fit(jtwo,jfour) - del
fit(jone,jfour) = fit(jone,jfour) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
indexll = indexll + 6.0
endif
if ((((targ(jone,jfour) + targ(jtwo,jthree)) <=                   (targ(jone,jtwo) +&
      & targ(jthree,jfour))) .and.                   ((targ(jone,jfour) + targ(jtwo,&
      &jthree)) <=                    (targ(jone,jthree) + targ(jtwo,jfour))))) then
if(((fit(jone,jfour) + fit(jtwo,jthree)) <=                   (fit(jone,jthree) + fi&
      &t(jtwo,jfour)))) then
del = 0.0
else
del = (fit(jone,jfour) + fit(jtwo,jthree) -                       (fit(jone,jthree) &
      &+ fit(jtwo,jfour)))/4.0
endif
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = del
work(int(indexll)+3,1) =  -del
work(int(indexll)+4,1) = -del
work(int(indexll)+5,1) = del
work(int(indexll)+6,1) = 0.0
fit(jone,jfour) = fit(jone,jfour) - del
fit(jtwo,jthree) = fit(jtwo,jthree) - del
fit(jone,jthree) = fit(jone,jthree) + del
fit(jtwo,jfour) = fit(jtwo,jfour) + del
indexll = indexll + 6.0
endif
cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree))                   + abs(&
      &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree))                   + abs(p5-fit&
      &(jtwo,jfour)) + abs(p6-fit(jthree,jfour))
enddo
enddo
enddo
enddo
enddo
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
fit(jtwo,jone) = fit(jone,jtwo)
enddo
enddo
aveprox = sum(spread(sum(prox,1),1,1))/(n*(n-1.0))
do i = 1,nint(n)
do j = 1,nint(n)
if(( i /= j)) then
proxave(i,j) = aveprox
else
proxave(i,j) = 0.0
endif
enddo
enddo
diff = sum(spread(sum((prox - fit)**2.0,1),1,1))
denom = sum(spread(sum((prox - proxave)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)
find = fit
aveprox = sum(spread(sum(prox,1),1,1))/(n*(n-1.0))
do i = 1,nint(n)
do j = 1,nint(n)
if(( i /= j)) then
proxave(i,j) = aveprox
else
proxave(i,j) = 0.0
endif
enddo
enddo
diff = sum(spread(sum((prox - find)**2.0,1),1,1))
denom = sum(spread(sum((prox - proxave)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)

deallocate( proxave, targ, fit)
return
end subroutine ms_atreefnd

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 03-May-2007 13:30:39     !
!------------------------------------------------------------------!
