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,targ_ptr,dumone_ptr
real(8), allocatable :: prox(:,:),inperm(:,:),targ(:,:),dumone(:,:)
integer(4) :: fit_ptr,vaf_ptr,addcon_ptr
real(8), allocatable :: fit(:,:),vaf(:,:),addcon(:,:)
!     Any other variables needed
integer prox_m,prox_n,inperm_m,inperm_n,targ_m,targ_n,dumone_m,dumone_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 4) then
  print *,'ms_cirfitac_ftarg requires 4 input arguments';return
elseif (nlhs .ne. 3) then
  print *,'ms_cirfitac_ftarg 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))
targ_m=mxGetm(prhs(3));targ_n=mxGetn(prhs(3))
dumone_m=mxGetm(prhs(4));dumone_n=mxGetn(prhs(4))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(prox_m,prox_n,0)
plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
plhs(3)=mxCreateDoubleMatrix(dumone_m,dumone_n,0)
fit_ptr=mxGetPr(plhs(1))
vaf_ptr=mxGetPr(plhs(2))
addcon_ptr=mxGetPr(plhs(3))
!     Copy right hand arguments to local arrays
prox_ptr=mxGetPr(prhs(1))
inperm_ptr=mxGetPr(prhs(2))
targ_ptr=mxGetPr(prhs(3))
dumone_ptr=mxGetPr(prhs(4))
!     Allocate and copy data to arrays
allocate(prox(prox_m,prox_n),inperm(inperm_m,inperm_n),targ(targ_m,targ_n),dumone(du&
      &mone_m,dumone_n))
allocate(fit(prox_m,prox_n),vaf(dumone_m,dumone_n),addcon(dumone_m,dumone_n))
call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n)
call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n)
call mxCopyPtrToReal8(targ_ptr,targ,targ_m*targ_n)
call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n)
!     Do the actual computations in a subroutine
call ms_cirfitac_ftarg(fit,vaf(1,1),addcon(1,1),prox,inperm,targ,dumone(1,1))
call mxCopyReal8ToPtr(fit,fit_ptr,prox_m*prox_n)
call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n)
call mxCopyReal8ToPtr(addcon,addcon_ptr,dumone_m*dumone_n)
deallocate(fit,vaf,addcon,prox,inperm,targ,dumone)
return

contains


subroutine ms_cirfitac_ftarg(fit,vaf,addcon,prox,inperm,targ,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) targ(:,:)
real(8) dumone
real(8) fit(:,:)
real(8) vaf
real(8) addcon
!     Matlab function pointers
!     All other local variables
real(8) denom
real(8) diff
real(8) aveprox
real(8) del
real(8) bbc
real(8) aac
real(8) bb
real(8) aa
real(8) nhaf
real(8) nhto
real(8) p1
integer jtwo
integer jone
real(8) indexll
real(8), allocatable :: coor(:,:)
real(8) cr
real(8) work(24360,1)
real(8) addconpv
integer j
integer i
real(8), allocatable :: fmadd(:,:)
real(8), allocatable :: proxave(:,:)
real(8) acondiff
real(8) n
!     Fill in vars going in and out
allocate( proxave(size(prox,1),size(prox,2)), fmadd(size(prox,1),size(prox,2)), coor&
      &(size(targ,2),1))
! --- Main computational routine. ---------------------------------!
n = size(prox,1)
acondiff = 1.0
addcon = 0.0
fit = prox([int(inperm)],[int(inperm)])
proxave = 0.
fmadd = 0.
do while (acondiff >= 1.0e-004)
do i = 1,nint(n)
do j = 1,nint(n)
if((i /= j)) then
fit(i,j) = prox(int(inperm(1,i)),int(inperm(1,j))) + addcon
else
fit(i,j) = 0.0
endif
enddo
enddo
addconpv = addcon
work = 0.
cr = 1.0
coor = 0.
do while (cr >= 1.0e-006)
cr = 0.0
indexll = 0.0
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
p1 = fit(jone,jtwo)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1)
if ((fit(jone,jtwo) < 0.0)) then
work(int(indexll)+1,1) = -fit(jone,jtwo)
fit(jone,jtwo) = 0.0
else
work(int(indexll)+1,1) = 0.0
endif
indexll = indexll + 1.0
cr = cr + abs(p1-fit(jone,jtwo))
enddo
enddo
do jone = 1,nint((n-2))
do jtwo = (jone+2),nint(n)
if (((jone /= 1.0) .or. (jtwo /= n))) then
nhto = jtwo - jone
nhaf = n - nhto
do i = 1,nint((n-1))
coor(i,1) = fit(i,i+1)
enddo
p1 = fit(jone,jtwo)
coor(int(n),1) = fit(1,int(n))
fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+3,1)
do i = 1,nint((n-1))
if ((i < jone)) then
fit(i,i+1) = fit(i,i+1) - work(int(indexll)+2,1)
endif
if (((i >= jone) .and. (i < jtwo))) then
fit(i,i+1) = fit(i,i+1) - work(int(indexll)+1,1)
endif
if ((i >= jtwo)) then
fit(i,i+1) = fit(i,i+1) - work(int(indexll)+2,1)
endif
enddo
fit(1,int(n)) = fit(1,int(n)) - work(int(indexll)+2,1)
aa = 0.0
bb = 0.0
aac = 0.0
bbc = 0.0
do i = 1,nint((n-1))
if ((i < jone)) then
bbc = bbc + targ(i,i+1)
endif
if (((i >= jone) .and. (i < jtwo))) then
aac = aac + targ(i,i+1)
endif
if ((i >= jtwo)) then
bbc = bbc + targ(i,i+1)
endif
enddo
bbc = bbc + targ(1,int(n))
do i = 1,nint((n-1))
if ((i < jone)) then
bb = bb + fit(i,i+1)
endif
if (((i >= jone) .and. (i < jtwo))) then
aa = aa + fit(i,i+1)
endif
if ((i >= jtwo)) then
bb = bb + fit(i,i+1)
endif
enddo
bb = bb + fit(1,int(n))
if ((aac <= bbc)) then
del = (aa - fit(jone,jtwo))/(nhto+1.0)
fit(jone,jtwo) = fit(jone,jtwo) + del
do i = 1,nint(nhto)
fit(jone+i-1,jone+i) =                         fit(jone+i-1,jone+i) - del
enddo
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = 0.0
work(int(indexll)+3,1) = del
cr = cr + abs(p1-fit(jone,jtwo))
indexll = indexll + 3.0
endif
if ((aac > bbc)) then
del = (bb - fit(jone,jtwo))/(nhaf+1.0)
fit(jone,jtwo) = fit(jone,jtwo) + del
if ((jone /= 1.0)) then
do i = 1,(jone-1)
fit(i,i+1) = fit(i,i+1) - del
enddo
endif
if ((jtwo /= n)) then
do i = jtwo,nint((n-1))
fit(i,i+1) = fit(i,i+1) - del
enddo
endif
fit(1,int(n)) = fit(1,int(n)) - del
work(int(indexll)+1,1) = 0.0
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) = del
cr = cr + abs(p1 - fit(jone,jtwo))
indexll = indexll +3.0
endif
endif
enddo
enddo
enddo
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
fit(jtwo,jone) = fit(jone,jtwo)
enddo
enddo
addcon = -sum(spread(sum(prox([int(inperm)],[int(inperm)]) - fit,1),1,1))/(n*(n-1.0))
acondiff = abs(addcon - addconpv)
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
fmadd(i,j) = fit(i,j) - addcon
else
proxave(i,j) = 0.0
fmadd(i,j) = 0.0
endif
enddo
enddo
diff = sum(spread(sum((prox([int(inperm)],[int(inperm)]) - (fmadd))**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)

deallocate( proxave, fmadd, coor)
return
end subroutine ms_cirfitac_ftarg

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 06-Dec-2007 10:32:29     !
!------------------------------------------------------------------!
