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
  module procedure callback0_
 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

 subroutine callback0_(funstr)
 character (len=*) :: funstr
 integer lhs(50), rhs(50)
  call mexCallMATLAB(0,NULL,0,NULL,funstr)
 end subroutine callback0_

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,dumone_ptr,dumtwo_ptr
real(8), allocatable :: prox(:,:),dumone(:,:),dumtwo(:,:)
integer(4) :: coord_ptr,permut_ptr,cumobfun_ptr,diff_ptr
real(8), allocatable :: coord(:,:),permut(:,:),cumobfun(:,:),diff(:,:)
!     Any other variables needed
integer prox_m,prox_n,dumone_m,dumone_n,dumtwo_m,dumtwo_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 3) then
  print *,'ms_uniscaldp requires 3 input arguments';return
elseif (nlhs .ne. 4) then
  print *,'ms_uniscaldp requires 4 output arguments';return
endif
!     Get the sizes of all the input variables
prox_m=mxGetm(prhs(1));prox_n=mxGetn(prhs(1))
dumone_m=mxGetm(prhs(2));dumone_n=mxGetn(prhs(2))
dumtwo_m=mxGetm(prhs(3));dumtwo_n=mxGetn(prhs(3))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
plhs(3)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
plhs(4)=mxCreateDoubleMatrix(dumtwo_m,dumtwo_n,0)
coord_ptr=mxGetPr(plhs(1))
permut_ptr=mxGetPr(plhs(2))
cumobfun_ptr=mxGetPr(plhs(3))
diff_ptr=mxGetPr(plhs(4))
!     Copy right hand arguments to local arrays
prox_ptr=mxGetPr(prhs(1))
dumone_ptr=mxGetPr(prhs(2))
dumtwo_ptr=mxGetPr(prhs(3))
!     Allocate and copy data to arrays
allocate(prox(prox_m,prox_n),dumone(dumone_m,dumone_n),dumtwo(dumtwo_m,dumtwo_n))
allocate(coord(dumone_m,dumone_n),permut(dumone_m,dumone_n),cumobfun(dumone_m,dumone&
      &_n),diff(dumtwo_m,dumtwo_n))
call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n)
call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n)
call mxCopyPtrToReal8(dumtwo_ptr,dumtwo,dumtwo_m*dumtwo_n)
!     Do the actual computations in a subroutine
call ms_uniscaldp(coord,permut,cumobfun,diff(1,1),prox,dumone,dumtwo(1,1))
call mxCopyReal8ToPtr(coord,coord_ptr,dumone_m*dumone_n)
call mxCopyReal8ToPtr(permut,permut_ptr,dumone_m*dumone_n)
call mxCopyReal8ToPtr(cumobfun,cumobfun_ptr,dumone_m*dumone_n)
call mxCopyReal8ToPtr(diff,diff_ptr,dumtwo_m*dumtwo_n)
deallocate(coord,permut,cumobfun,diff,prox,dumone,dumtwo)
return

contains


subroutine ms_uniscaldp(coord,permut,cumobfun,diff,prox,dumone,dumtwo)
!     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) dumone(:,:)
real(8) dumtwo
real(8) coord(:,:)
real(8) permut(:,:)
real(8) cumobfun(:,:)
real(8) diff
!     Matlab function pointers
!     All other local variables
real(8) lastint
real(8) comp
real(8) idxtmp
real(8) temp
real(8) incre
real(8) ione
real(8) sum
real(8) jone
real(8) nk
real(8) idxcomp
integer jj
real(8) nh
real(8) m2
real(8) nfirst
integer k
real(8) mtc
real(8) index
integer j
real(8) rowsum
integer i
real(8), allocatable :: subcomp(:,:)
real(8), allocatable :: sub(:,:)
real(8) idxstore(16777215,1)
real(8) valstore(16777215,1)
real(8) nsum
real(8) bigneg
real(8) n
!     Fill in vars going in and out
allocate( sub(size(dumone,1),size(dumone,2)), subcomp(size(dumone,1),size(dumone,2)))
! --- Main computational routine. ---------------------------------!
call mlcall0('tic')
n = size(prox,1)
bigneg = -1.0e+20
nsum = (2.0**n) - 1.0
valstore = 1.
idxstore = 0.
sub = 0.
subcomp = 0.
coord = 0.
permut = 0.
cumobfun = 0.
do i = 1,nint(n)
rowsum = 0.0
do j = 1,nint(n)
rowsum = rowsum + prox(i,j)
enddo
index = 2.0**(i-1.0)
idxstore(int(index),1) = i
valstore(int(index),1) = rowsum**2.0
enddo
mtc = 0.0
do k = 1,nint((n-1))
nfirst = 0.0
if((mtc == 0.0)) then
m2 = 0.0
nh = k
do j = 1,nint(nh)
sub(k+j-int(nh),1) = m2 + j
enddo
if((sub(1,1) /= (n-k+1.0))) then
mtc = 1.0
endif
endif
do while (mtc == 1)
if((nfirst == 1.0)) then
if((m2 < (n-nh))) then
nh = 0.0
endif
nh = nh + 1.0
m2 = sub(k+1-int(nh),1)
do j = 1,nint(nh)
sub(k+j-int(nh),1) = m2 + j
enddo
if((sub(1,1) /= (n-k+1.0))) then
mtc = 1.0
else
mtc = 0.0
endif
endif
if((nfirst == 0.0)) then
nfirst = 1.0
endif
index = 0.0
do i = 1,k
index = index + 2.0**(sub(i,1)-1.0)
enddo
jj = 1.0
subcomp = 0.
do i = 1,nint(n)
idxcomp = 1.0
do j = 1,k
if((sub(j,1) == i)) then
idxcomp = 0.0
endif
enddo
if((idxcomp == 1.0)) then
subcomp(jj,1) = i
jj = jj + 1.0
endif
enddo
nk = n - k
do jj = 1,nint(nk)
jone = subcomp(jj,1)
sum = 0.0
do i = 1,k
ione = sub(i,1)
sum = sum + prox(int(jone),int(ione))
enddo
do i = 1,nint(nk)
ione = subcomp(i,1)
if((ione /= jone)) then
sum = sum - prox(int(jone),int(ione))
endif
enddo
incre = sum**2.0
temp = valstore(int(index),1) + incre
idxtmp = index + 2.0**(jone-1.0)
comp = valstore(int(idxtmp),1)
if((temp > comp)) then
valstore(int(idxtmp),1) = temp
idxstore(int(idxtmp),1) = jone
endif
enddo
enddo
enddo
permut(int(n),1) = idxstore(int(nsum),1)
cumobfun(int(n),1) = valstore(int(nsum),1)
index = nsum
lastint = permut(int(n),1)
do i = 1,nint((n-1))
index = index - (2.0**(lastint-1.0))
lastint = idxstore(int(index),1)
permut(int(n)-i,1) = lastint
cumobfun(int(n)-i,1) = valstore(int(index),1)
enddo
do i = 1,nint(n)
do j = 1,nint(n)
if((i > j)) then
coord(i,1) = coord(i,1) + prox(int(permut(i,1)),int(permut(j,1)))
endif
if((i < j)) then
coord(i,1) = coord(i,1) - prox(int(permut(i,1)),int(permut(j,1)))
endif
enddo
coord(i,1) = coord(i,1)/n
enddo
diff = 0.0
do i = 1,nint((n-1))
do j = (i+1),nint(n)
diff = diff + (prox(int(permut(i,1)),int(permut(j,1))) -          abs(coord(i,1) - c&
      &oord(j,1)))**2.0
enddo
enddo
call mlcall0('toc')

deallocate( sub, subcomp)
return
end subroutine ms_uniscaldp

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 27-Jul-2007 18:05:41     !
!------------------------------------------------------------------!
