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,kblock_ptr
real(8), allocatable :: prox(:,:),inperm(:,:),kblock(:,:)
integer(4) :: find_ptr,vaf_ptr,outperm_ptr
real(8), allocatable :: find(:,:),vaf(:,:),outperm(:,:)
!     Any other variables needed
integer prox_m,prox_n,inperm_m,inperm_n,kblock_m,kblock_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 3) then
  print *,'ms_arobfnd requires 3 input arguments';return
elseif (nlhs .ne. 3) then
  print *,'ms_arobfnd requires 3 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))
kblock_m=mxGetm(prhs(3));kblock_n=mxGetn(prhs(3))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(prox_m,prox_n,0)
plhs(2)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
plhs(3)=mxCreateDoubleMatrix(inperm_m,inperm_n,0)
find_ptr=mxGetPr(plhs(1))
vaf_ptr=mxGetPr(plhs(2))
outperm_ptr=mxGetPr(plhs(3))
!     Copy right hand arguments to local arrays
prox_ptr=mxGetPr(prhs(1))
inperm_ptr=mxGetPr(prhs(2))
kblock_ptr=mxGetPr(prhs(3))
!     Allocate and copy data to arrays
allocate(prox(prox_m,prox_n),inperm(inperm_m,inperm_n),kblock(kblock_m,kblock_n))
allocate(find(prox_m,prox_n),vaf(kblock_m,kblock_n),outperm(inperm_m,inperm_n))
call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n)
call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n)
call mxCopyPtrToReal8(kblock_ptr,kblock,kblock_m*kblock_n)
!     Do the actual computations in a subroutine
call ms_arobfnd(find,vaf(1,1),outperm,prox,inperm,kblock(1,1))
call mxCopyReal8ToPtr(find,find_ptr,prox_m*prox_n)
call mxCopyReal8ToPtr(vaf,vaf_ptr,kblock_m*kblock_n)
call mxCopyReal8ToPtr(outperm,outperm_ptr,inperm_m*inperm_n)
deallocate(find,vaf,outperm,prox,inperm,kblock)
return

contains


subroutine ms_arobfnd(find,vaf,outperm,prox,inperm,kblock)
!     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) kblock
real(8) find(:,:)
real(8) vaf
real(8) outperm(:,:)
!     Matlab function pointers
!     All other local variables
real(8), allocatable :: prevperm(:,:)
real(8) nprevperm
real(8) denom
real(8) diff
real(8) aveprox
real(8) ave
real(8) p2
real(8) p1
real(8) indexll
real(8), allocatable :: proxave(:,:)
real(8) cr
real(8), allocatable :: fit(:,:)
real(8) work(24360,1)
integer jone
integer jtwo
integer nlimlow
integer insertpt
real(8) rawindex
real(8) tryindex
real(8), allocatable :: intrperm(:,:)
integer k
real(8) nchange
integer iterate
real(8) begindex
real(8) index
integer j
integer i
real(8), allocatable :: targ(:,:)
real(8) n
!     Fill in vars going in and out
allocate( targ(size(prox,1),size(prox,2)), intrperm(size(inperm,1),size(inperm,2)), &
      &fit(size(prox,1),size(prox,2)), proxave(size(prox,1),size(prox,2)), prevperm(&
      &size(inperm,1),size(inperm,2)))
! --- Main computational routine. ---------------------------------!
n = size(prox,1)
targ = 0.
do i = 1,nint(n-1)
do j = (i+1),nint(n)
targ(i,j) = abs(i-j)
targ(j,i) = targ(i,j)
enddo
enddo
n = size(prox,1)
outperm = inperm
index = 1.0
begindex = sum(spread(sum(prox([int(inperm)],[int(inperm)])*targ,1),1,1))
do iterate = 1,2
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint((n-1))
do j = (k+1),nint(n)
intrperm = outperm
intrperm(1,k) = outperm(1,j)
intrperm(1,j) = outperm(1,k)
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if ((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index + 1.0
endif
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint(kblock)
do insertpt = 1,nint((n+1))
do nlimlow = 1,nint((n+1-k))
intrperm = outperm
if ((nlimlow > insertpt)) then
jtwo = 0.0
do j = insertpt,(insertpt+k-1)
intrperm(1,j) =outperm(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = (insertpt+k),(nlimlow+k-1)
intrperm(1,j) = outperm(1,insertpt+jone)
jone = jone + 1.0
enddo
endif
if (((nlimlow+k) < insertpt)) then
jtwo = 0.0
do j = (insertpt-k),(insertpt-1)
intrperm(1,j) = outperm(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = nlimlow,(insertpt-k-1)
intrperm(1,j) = outperm(1,nlimlow+k+jone)
jone = jone + 1.0
enddo
endif
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index +1.0
endif
enddo
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 2,nint(kblock)
do nlimlow = 1,nint((n+1-k))
intrperm = outperm
do j = 1,k
intrperm(1,nlimlow+j-1) = outperm(1,nlimlow+k-j)
enddo
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index + 1.0
endif
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
enddo
inperm = outperm
n = size(prox,1)
work = 0.
fit = prox([int(inperm)],[int(inperm)])
cr = 1.0
proxave = 0.
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
do jone = 1,nint((n-2))
do jtwo = (jone+1),nint((n-1))
p1 = fit(jone,jtwo)
p2 = fit(jone,jtwo+1)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(int(indexll)+2,1)
if((fit(jone,jtwo) <= fit(jone,jtwo+1))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jone,jtwo) > fit(jone,jtwo+1))) then
ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0
work(int(indexll)+1,1) = ave - fit(jone,jtwo)
work(int(indexll)+2,1) = ave - fit(jone,jtwo+1)
fit(jone,jtwo) = ave
fit(jone,jtwo+1) = ave
endif
cr = cr + abs(p1-fit(jone,jtwo)) +             abs(p2-fit(jone,jtwo+1))
indexll = indexll + 2.0
enddo
enddo
do jone = 3,nint(n)
do jtwo = 1,(jone-2)
p1 = fit(jtwo,jone)
p2 = fit(jtwo+1,jone)
fit(jtwo,jone) = fit(jtwo,jone) - work(int(indexll)+1,1)
fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(int(indexll)+2,1)
if((fit(jtwo+1,jone) <= fit(jtwo,jone))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jtwo+1,jone) > fit(jtwo,jone))) then
ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0
work(int(indexll)+1,1) = ave - fit(jtwo,jone)
work(int(indexll)+2,1) = ave - fit(jtwo+1,jone)
fit(jtwo,jone) = ave
fit(jtwo+1,jone) = ave
endif
cr = cr + abs(p1-fit(jtwo,jone)) +             abs(p2-fit(jtwo+1,jone))
indexll = indexll + 2.0
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([int(inperm)],[int(inperm)]) - fit)**2.0,1),1,1))
denom = sum(spread(sum((prox([int(inperm)],[int(inperm)]) - proxave)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)
nprevperm = 1.0
do while (nprevperm == 1)
nprevperm = 0.0
prevperm = outperm
inperm = outperm
n = size(prox,1)
outperm = inperm
index = 1.0
targ = fit
begindex = sum(spread(sum(prox([int(inperm)],[int(inperm)])*targ,1),1,1))
do iterate = 1,2
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint((n-1))
do j = (k+1),nint(n)
intrperm = outperm
intrperm(1,k) = outperm(1,j)
intrperm(1,j) = outperm(1,k)
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if ((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index + 1.0
endif
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint(kblock)
do insertpt = 1,nint((n+1))
do nlimlow = 1,nint((n+1-k))
intrperm = outperm
if ((nlimlow > insertpt)) then
jtwo = 0.0
do j = insertpt,(insertpt+k-1)
intrperm(1,j) =outperm(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = (insertpt+k),(nlimlow+k-1)
intrperm(1,j) = outperm(1,insertpt+jone)
jone = jone + 1.0
enddo
endif
if (((nlimlow+k) < insertpt)) then
jtwo = 0.0
do j = (insertpt-k),(insertpt-1)
intrperm(1,j) = outperm(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = nlimlow,(insertpt-k-1)
intrperm(1,j) = outperm(1,nlimlow+k+jone)
jone = jone + 1.0
enddo
endif
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index +1.0
endif
enddo
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 2,nint(kblock)
do nlimlow = 1,nint((n+1-k))
intrperm = outperm
do j = 1,k
intrperm(1,nlimlow+j-1) = outperm(1,nlimlow+k-j)
enddo
tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1))
if((tryindex > (begindex + 1.0e-008))) then
nchange = 1.0
begindex = tryindex
outperm = intrperm
index = index + 1.0
endif
enddo
enddo
enddo
rawindex = begindex
nchange = 1.0
enddo
inperm = outperm
n = size(prox,1)
work = 0.
fit = prox([int(inperm)],[int(inperm)])
cr = 1.0
proxave = 0.
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
do jone = 1,nint((n-2))
do jtwo = (jone+1),nint((n-1))
p1 = fit(jone,jtwo)
p2 = fit(jone,jtwo+1)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(int(indexll)+2,1)
if((fit(jone,jtwo) <= fit(jone,jtwo+1))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jone,jtwo) > fit(jone,jtwo+1))) then
ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0
work(int(indexll)+1,1) = ave - fit(jone,jtwo)
work(int(indexll)+2,1) = ave - fit(jone,jtwo+1)
fit(jone,jtwo) = ave
fit(jone,jtwo+1) = ave
endif
cr = cr + abs(p1-fit(jone,jtwo)) +             abs(p2-fit(jone,jtwo+1))
indexll = indexll + 2.0
enddo
enddo
do jone = 3,nint(n)
do jtwo = 1,(jone-2)
p1 = fit(jtwo,jone)
p2 = fit(jtwo+1,jone)
fit(jtwo,jone) = fit(jtwo,jone) - work(int(indexll)+1,1)
fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(int(indexll)+2,1)
if((fit(jtwo+1,jone) <= fit(jtwo,jone))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jtwo+1,jone) > fit(jtwo,jone))) then
ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0
work(int(indexll)+1,1) = ave - fit(jtwo,jone)
work(int(indexll)+2,1) = ave - fit(jtwo+1,jone)
fit(jtwo,jone) = ave
fit(jtwo+1,jone) = ave
endif
cr = cr + abs(p1-fit(jtwo,jone)) +             abs(p2-fit(jtwo+1,jone))
indexll = indexll + 2.0
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([int(inperm)],[int(inperm)]) - fit)**2.0,1),1,1))
denom = sum(spread(sum((prox([int(inperm)],[int(inperm)]) - proxave)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)
if ((sum(prevperm - outperm) /= 0.0)) then
nprevperm = 1.0
endif
enddo
find = fit

deallocate( targ, intrperm, fit, proxave, prevperm)
return
end subroutine ms_arobfnd

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 03-May-2007 13:26:16     !
!------------------------------------------------------------------!
