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) :: proxpermuttm_ptr,fittedtm_ptr,dumone_ptr
real(8), allocatable :: proxpermuttm(:,:),fittedtm(:,:),dumone(:,:)
integer(4) :: monproxpermuttm_ptr,vaf_ptr,diff_ptr
real(8), allocatable :: monproxpermuttm(:,:),vaf(:,:),diff(:,:)
!     Any other variables needed
integer proxpermuttm_m,proxpermuttm_n,fittedtm_m,fittedtm_n,dumone_m,dumone_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 3) then
  print *,'ms_proxmontm requires 3 input arguments';return
elseif (nlhs .ne. 3) then
  print *,'ms_proxmontm requires 3 output arguments';return
endif
!     Get the sizes of all the input variables
proxpermuttm_m=mxGetm(prhs(1));proxpermuttm_n=mxGetn(prhs(1))
fittedtm_m=mxGetm(prhs(2));fittedtm_n=mxGetn(prhs(2))
dumone_m=mxGetm(prhs(3));dumone_n=mxGetn(prhs(3))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(proxpermuttm_m,proxpermuttm_n,0)
plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
plhs(3)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
monproxpermuttm_ptr=mxGetPr(plhs(1))
vaf_ptr=mxGetPr(plhs(2))
diff_ptr=mxGetPr(plhs(3))
!     Copy right hand arguments to local arrays
proxpermuttm_ptr=mxGetPr(prhs(1))
fittedtm_ptr=mxGetPr(prhs(2))
dumone_ptr=mxGetPr(prhs(3))
!     Allocate and copy data to arrays
allocate(proxpermuttm(proxpermuttm_m,proxpermuttm_n),fittedtm(fittedtm_m,fittedtm_n)&
      &,dumone(dumone_m,dumone_n))
allocate(monproxpermuttm(proxpermuttm_m,proxpermuttm_n),vaf(dumone_m,dumone_n),diff(&
      &dumone_m,dumone_n))
call mxCopyPtrToReal8(proxpermuttm_ptr,proxpermuttm,proxpermuttm_m*proxpermuttm_n)
call mxCopyPtrToReal8(fittedtm_ptr,fittedtm,fittedtm_m*fittedtm_n)
call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n)
!     Do the actual computations in a subroutine
call ms_proxmontm(monproxpermuttm,vaf(1,1),diff(1,1),proxpermuttm,fittedtm,dumone(1,&
      &1))
call mxCopyReal8ToPtr(monproxpermuttm,monproxpermuttm_ptr,proxpermuttm_m*proxpermutt&
      &m_n)
call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n)
call mxCopyReal8ToPtr(diff,diff_ptr,dumone_m*dumone_n)
deallocate(monproxpermuttm,vaf,diff,proxpermuttm,fittedtm,dumone)
return

contains


subroutine ms_proxmontm(monproxpermuttm,vaf,diff,proxpermuttm,fittedtm,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) proxpermuttm(:,:)
real(8) fittedtm(:,:)
real(8) dumone
real(8) monproxpermuttm(:,:)
real(8) vaf
real(8) diff
!     Matlab function pointers
!     All other local variables
real(8) denom
real(8) avefit
real(8) ave
real(8) p2
real(8) p1
integer jfour
integer jthree
integer jtwo
integer jone
real(8) indexll
real(8) cr
real(8), allocatable :: fit(:,:)
real(8), allocatable :: targ(:,:)
real(8) work(2433600,1)
real(8) n
real(8) ncol
real(8) nrow
!     Fill in vars going in and out
allocate( targ(size(proxpermuttm,1),size(proxpermuttm,2)), fit(size(proxpermuttm,1),&
      &size(proxpermuttm,2)))
! --- Main computational routine. ---------------------------------!
nrow = size(proxpermuttm,1)
ncol = size(proxpermuttm,2)
n = nrow + ncol
work = 0.
targ = proxpermuttm
fit = fittedtm
cr = 1.0
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
do jone = 1,nint(nrow)
do jtwo = 1,nint(ncol)
do jthree = 1,nint(nrow)
do jfour = 1,nint(ncol)
if (((jone /= jthree) .or. (jtwo /= jfour))) then
p1 = fit(jone,jtwo)
p2 = fit(jthree,jfour)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+2,1)
if (((abs(targ(jone,jtwo) - targ(jthree,jfour))                         > 1.0e-006) &
      &.and. (targ(jone,jtwo) <                         targ(jthree,jfour)))) then
if((fit(jone,jtwo) <= fit(jthree,jfour))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jone,jtwo) > fit(jthree,jfour))) then
ave = (fit(jone,jtwo) + fit(jthree,jfour))/2.0
work(int(indexll)+1,1) = ave - fit(jone,jtwo)
work(int(indexll)+2,1) = ave - fit(jthree,jfour)
fit(jone,jtwo) = ave
fit(jthree,jfour) = ave
endif
endif
if (((abs(targ(jone,jtwo) - targ(jthree,jfour))                         > 1.0e-006) &
      &.and. (targ(jone,jtwo) >                         targ(jthree,jfour)))) then
if((fit(jone,jtwo) >= fit(jthree,jfour))) then
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = 0.0
endif
if ((fit(jone,jtwo) < fit(jthree,jfour))) then
ave = (fit(jone,jtwo) + fit(jthree,jfour))/2.0
work(int(indexll)+1,1) = ave - fit(jone,jtwo)
work(int(indexll)+2,1) = ave - fit(jthree,jfour)
fit(jone,jtwo) = ave
fit(jthree,jfour) = ave
endif
endif
cr = cr + abs(p1-fit(jone,jtwo)) +                   abs(p2-fit(jthree,jfour))
endif
indexll = indexll + 2.0
enddo
enddo
enddo
enddo
enddo
avefit = sum(spread(sum(fit,1),1,1))/(nrow*ncol)
diff = sum(spread(sum((fit - fittedtm)**2.0,1),1,1))
denom = sum(spread(sum((fit - avefit)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)
monproxpermuttm = fit

deallocate( targ, fit)
return
end subroutine ms_proxmontm

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 03-May-2007 13:24:32     !
!------------------------------------------------------------------!
