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) :: proxtm_ptr,targ_ptr,inperm_ptr,kblock_ptr
real(8), allocatable :: proxtm(:,:),targ(:,:),inperm(:,:),kblock(:,:)
integer(4) :: outperm_ptr,rawindex_ptr,index_ptr,squareprox_ptr
real(8), allocatable :: outperm(:,:),rawindex(:,:),index(:,:),squareprox(:,:)
!     Any other variables needed
integer proxtm_m,proxtm_n,targ_m,targ_n,inperm_m,inperm_n,kblock_m,kblock_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 4) then
  print *,'ms_ordertm requires 4 input arguments';return
elseif (nlhs .ne. 4) then
  print *,'ms_ordertm requires 4 output arguments';return
endif
!     Get the sizes of all the input variables
proxtm_m=mxGetm(prhs(1));proxtm_n=mxGetn(prhs(1))
targ_m=mxGetm(prhs(2));targ_n=mxGetn(prhs(2))
inperm_m=mxGetm(prhs(3));inperm_n=mxGetn(prhs(3))
kblock_m=mxGetm(prhs(4));kblock_n=mxGetn(prhs(4))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(inperm_m,inperm_n,0)
plhs(2)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
plhs(3)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
plhs(4)=mxCreateDoubleMatrix(targ_m,targ_n,0)
outperm_ptr=mxGetPr(plhs(1))
rawindex_ptr=mxGetPr(plhs(2))
index_ptr=mxGetPr(plhs(3))
squareprox_ptr=mxGetPr(plhs(4))
!     Copy right hand arguments to local arrays
proxtm_ptr=mxGetPr(prhs(1))
targ_ptr=mxGetPr(prhs(2))
inperm_ptr=mxGetPr(prhs(3))
kblock_ptr=mxGetPr(prhs(4))
!     Allocate and copy data to arrays
allocate(proxtm(proxtm_m,proxtm_n),targ(targ_m,targ_n),inperm(inperm_m,inperm_n),kbl&
      &ock(kblock_m,kblock_n))
allocate(outperm(inperm_m,inperm_n),rawindex(kblock_m,kblock_n),index(kblock_m,kbloc&
      &k_n),squareprox(targ_m,targ_n))
call mxCopyPtrToReal8(proxtm_ptr,proxtm,proxtm_m*proxtm_n)
call mxCopyPtrToReal8(targ_ptr,targ,targ_m*targ_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_ordertm(outperm,rawindex(1,1),index(1,1),squareprox,proxtm,targ,inperm,kbloc&
      &k(1,1))
call mxCopyReal8ToPtr(outperm,outperm_ptr,inperm_m*inperm_n)
call mxCopyReal8ToPtr(rawindex,rawindex_ptr,kblock_m*kblock_n)
call mxCopyReal8ToPtr(index,index_ptr,kblock_m*kblock_n)
call mxCopyReal8ToPtr(squareprox,squareprox_ptr,targ_m*targ_n)
deallocate(outperm,rawindex,index,squareprox,proxtm,targ,inperm,kblock)
return

contains


subroutine ms_ordertm(outperm,rawindex,index,squareprox,proxtm,targ,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) proxtm(:,:)
real(8) targ(:,:)
real(8) inperm(:,:)
real(8) kblock
real(8) outperm(:,:)
real(8) rawindex
real(8) index
real(8) squareprox(:,:)
!     Matlab function pointers
!     All other local variables
real(8) jone
real(8) jtwo
integer nlimlow
integer insertpt
real(8) tryindex
real(8), allocatable :: intrperm(:,:)
integer j
integer k
real(8) nchange
integer iterate
real(8) begindex
integer icol
integer irow
real(8) aveprox
real(8), allocatable :: prox(:,:)
real(8) n
real(8) ncol
real(8) nrow
!     Fill in vars going in and out
allocate( prox(size(targ,1),size(targ,2)), intrperm(size(inperm,1),size(inperm,2)))
! --- Main computational routine. ---------------------------------!
nrow = size(proxtm,1)
ncol = size(proxtm,2)
n = nrow + ncol
prox = 0.
aveprox = sum(spread(sum(proxtm,1),1,1))/(nrow*ncol)
do irow = 1,nint(nrow)
do icol = 1,nint(ncol)
prox(irow,icol+int(nrow)) = proxtm(irow,icol) - aveprox
prox(icol+int(nrow),irow) = proxtm(irow,icol) - aveprox
enddo
enddo
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+int(jtwo))
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = (insertpt+k),(nlimlow+k-1)
intrperm(1,j) = outperm(1,insertpt+int(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+int(jtwo))
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = nlimlow,(insertpt-k-1)
intrperm(1,j) = outperm(1,nlimlow+k+int(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
squareprox = prox

deallocate( prox, intrperm)
return
end subroutine ms_ordertm

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 16-Feb-2007 11:25:11     !
!------------------------------------------------------------------!
