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,targone_ptr,targtwo_ptr,inpermone_ptr,inpermtwo_ptr,kblock_pt&
      &r,nopt_ptr
real(8), allocatable :: prox(:,:),targone(:,:),targtwo(:,:),inpermone(:,:),inpermtwo(:,&
      &:),kblock(:,:),nopt(:,:)
integer(4) :: outpermone_ptr,outpermtwo_ptr,coordone_ptr,coordtwo_ptr,fitone_ptr,fit&
      &two_ptr,addconone_ptr,addcontwo_ptr,vaf_ptr
real(8), allocatable :: outpermone(:,:),outpermtwo(:,:),coordone(:,:),coordtwo(:,:),fit&
      &one(:,:),fittwo(:,:),addconone(:,:),addcontwo(:,:),vaf(:,:)
!     Any other variables needed
integer prox_m,prox_n,targone_m,targone_n,targtwo_m,targtwo_n,inpermone_m,inpermone_&
      &n,inpermtwo_m,inpermtwo_n,kblock_m,kblock_n,nopt_m,nopt_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 7) then
  print *,'ms_biscalqa requires 7 input arguments';return
elseif (nlhs .ne. 9) then
  print *,'ms_biscalqa requires 9 output arguments';return
endif
!     Get the sizes of all the input variables
prox_m=mxGetm(prhs(1));prox_n=mxGetn(prhs(1))
targone_m=mxGetm(prhs(2));targone_n=mxGetn(prhs(2))
targtwo_m=mxGetm(prhs(3));targtwo_n=mxGetn(prhs(3))
inpermone_m=mxGetm(prhs(4));inpermone_n=mxGetn(prhs(4))
inpermtwo_m=mxGetm(prhs(5));inpermtwo_n=mxGetn(prhs(5))
kblock_m=mxGetm(prhs(6));kblock_n=mxGetn(prhs(6))
nopt_m=mxGetm(prhs(7));nopt_n=mxGetn(prhs(7))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(inpermone_m,inpermone_n,0)
plhs(2)=mxCreateDoubleMatrix(inpermone_m,inpermone_n,0)
plhs(3)=mxCreateDoubleMatrix(inpermtwo_n,nopt_n,0)
plhs(4)=mxCreateDoubleMatrix(inpermtwo_n,nopt_n,0)
plhs(5)=mxCreateDoubleMatrix(prox_m,prox_n,0)
plhs(6)=mxCreateDoubleMatrix(prox_m,prox_n,0)
plhs(7)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
plhs(8)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
plhs(9)=mxCreateDoubleMatrix(kblock_m,kblock_n,0)
outpermone_ptr=mxGetPr(plhs(1))
outpermtwo_ptr=mxGetPr(plhs(2))
coordone_ptr=mxGetPr(plhs(3))
coordtwo_ptr=mxGetPr(plhs(4))
fitone_ptr=mxGetPr(plhs(5))
fittwo_ptr=mxGetPr(plhs(6))
addconone_ptr=mxGetPr(plhs(7))
addcontwo_ptr=mxGetPr(plhs(8))
vaf_ptr=mxGetPr(plhs(9))
!     Copy right hand arguments to local arrays
prox_ptr=mxGetPr(prhs(1))
targone_ptr=mxGetPr(prhs(2))
targtwo_ptr=mxGetPr(prhs(3))
inpermone_ptr=mxGetPr(prhs(4))
inpermtwo_ptr=mxGetPr(prhs(5))
kblock_ptr=mxGetPr(prhs(6))
nopt_ptr=mxGetPr(prhs(7))
!     Allocate and copy data to arrays
allocate(prox(prox_m,prox_n),targone(targone_m,targone_n),targtwo(targtwo_m,targtwo_&
      &n),inpermone(inpermone_m,inpermone_n),inpermtwo(inpermtwo_m,inpermtwo_n),kblo&
      &ck(kblock_m,kblock_n),nopt(nopt_m,nopt_n))
allocate(outpermone(inpermone_m,inpermone_n),outpermtwo(inpermone_m,inpermone_n),coo&
      &rdone(inpermtwo_n,nopt_n),coordtwo(inpermtwo_n,nopt_n),fitone(prox_m,prox_n),&
      &fittwo(prox_m,prox_n),addconone(kblock_m,kblock_n),addcontwo(kblock_m,kblock_&
      &n),vaf(kblock_m,kblock_n))
call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n)
call mxCopyPtrToReal8(targone_ptr,targone,targone_m*targone_n)
call mxCopyPtrToReal8(targtwo_ptr,targtwo,targtwo_m*targtwo_n)
call mxCopyPtrToReal8(inpermone_ptr,inpermone,inpermone_m*inpermone_n)
call mxCopyPtrToReal8(inpermtwo_ptr,inpermtwo,inpermtwo_m*inpermtwo_n)
call mxCopyPtrToReal8(kblock_ptr,kblock,kblock_m*kblock_n)
call mxCopyPtrToReal8(nopt_ptr,nopt,nopt_m*nopt_n)
!     Do the actual computations in a subroutine
call ms_biscalqa(outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo,addconone(1,1&
      &),addcontwo(1,1),vaf(1,1),prox,targone,targtwo,inpermone,inpermtwo,kblock(1,1&
      &),nopt(1,1))
call mxCopyReal8ToPtr(outpermone,outpermone_ptr,inpermone_m*inpermone_n)
call mxCopyReal8ToPtr(outpermtwo,outpermtwo_ptr,inpermone_m*inpermone_n)
call mxCopyReal8ToPtr(coordone,coordone_ptr,inpermtwo_n*nopt_n)
call mxCopyReal8ToPtr(coordtwo,coordtwo_ptr,inpermtwo_n*nopt_n)
call mxCopyReal8ToPtr(fitone,fitone_ptr,prox_m*prox_n)
call mxCopyReal8ToPtr(fittwo,fittwo_ptr,prox_m*prox_n)
call mxCopyReal8ToPtr(addconone,addconone_ptr,kblock_m*kblock_n)
call mxCopyReal8ToPtr(addcontwo,addcontwo_ptr,kblock_m*kblock_n)
call mxCopyReal8ToPtr(vaf,vaf_ptr,kblock_m*kblock_n)
deallocate(outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo,addconone,addcontwo&
      &,vaf,prox,targone,targtwo,inpermone,inpermtwo,kblock,nopt)
return

contains


subroutine ms_biscalqa(outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo,addcono&
      &ne,addcontwo,vaf,prox,targone,targtwo,inpermone,inpermtwo,kblock,nopt)
!     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) targone(:,:)
real(8) targtwo(:,:)
real(8) inpermone(:,:)
real(8) inpermtwo(:,:)
real(8) kblock
real(8) nopt
real(8) outpermone(:,:)
real(8) outpermtwo(:,:)
real(8) coordone(:,:)
real(8) coordtwo(:,:)
real(8) fitone(:,:)
real(8) fittwo(:,:)
real(8) addconone
real(8) addcontwo
real(8) vaf
!     Matlab function pointers
!     All other local variables
real(8) vaftwo
real(8) begindextwo
real(8) vafone
real(8), allocatable :: coord(:,:)
real(8) denom
real(8) diff
real(8) aveprox
real(8) del
real(8) p3
real(8) p2
real(8) p1
integer jthree
real(8) index
real(8) cr
real(8) work(24360,1)
real(8) addconpv
real(8), allocatable :: fmadd(:,:)
real(8), allocatable :: fit(:,:)
real(8) addcon
real(8) acondiff
integer jone
integer jtwo
integer nlimlow
integer insertpt
real(8) tryindex
real(8), allocatable :: intrperm(:,:)
integer k
real(8) nchange
real(8) begindexone
real(8) vafprev
real(8) vafdiff
integer j
integer i
real(8) aveproxorig
real(8), allocatable :: proxthree(:,:)
real(8), allocatable :: proxave(:,:)
real(8), allocatable :: proxaveorig(:,:)
real(8), allocatable :: proxtwo(:,:)
real(8), allocatable :: proxone(:,:)
real(8), allocatable :: fittwodim(:,:)
real(8), allocatable :: fitonedim(:,:)
real(8) n
!     Fill in vars going in and out
allocate( fitonedim(size(prox,1),size(prox,2)), fittwodim(size(prox,1),size(prox,2))&
      &, proxone(size(prox,1),size(prox,2)), proxtwo(size(prox,1),size(prox,2)), pro&
      &xaveorig(size(prox,1),size(prox,2)), proxave(size(prox,1),size(prox,2)), prox&
      &three(size(prox,1),size(prox,2)), intrperm(size(inpermone,1),size(inpermone,2&
      &)), fit(size(prox,1),size(prox,2)), fmadd(size(prox,1),size(prox,2)), coord(s&
      &ize(inpermtwo,2),1))
! --- Main computational routine. ---------------------------------!
n = size(prox,1)
outpermone = inpermone
outpermtwo = inpermtwo
coordone = 0.
coordtwo = 0.
fitone = targone
fittwo = targtwo
addconone = 0.0
addcontwo = 0.0
fitonedim = 0.
fittwodim = 0.
proxone = prox
proxtwo = 0.
proxaveorig = 0.
proxave = 0.
proxthree = 0.
aveproxorig = 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
proxaveorig(i,j) = aveproxorig
else
proxaveorig(i,j) = 0.0
endif
enddo
enddo
vafdiff = 1.0
vaf = 0.0
do while (vafdiff >= 1.0e-005)
vafprev = vaf
if ((nopt == 1.0)) then
begindexone = sum(spread(sum(proxone([int(outpermone)],[int(outpermone)])*fitone,1),&
      &1,1))
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint((n-1))
do j = (k+1),nint(n)
intrperm = outpermone
intrperm(1,k) = outpermone(1,j)
intrperm(1,j) = outpermone(1,k)
tryindex = sum(spread(sum(proxone([int(intrperm)],[int(intrperm)])*fitone,1),1,1))
if((tryindex > (begindexone + 1.0e-008))) then
nchange = 1.0
begindexone = tryindex
outpermone = intrperm
endif
enddo
enddo
do k = 1,nint(kblock)
do insertpt = 1,nint((n+1))
do nlimlow = 1,nint((n+1-k))
intrperm = outpermone
if ((nlimlow > insertpt)) then
jtwo = 0.0
do j = insertpt,(insertpt+k-1)
intrperm(1,j) =outpermone(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = (insertpt+k),(nlimlow+k-1)
intrperm(1,j) = outpermone(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) = outpermone(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = nlimlow,(insertpt-k-1)
intrperm(1,j) = outpermone(1,nlimlow+k+jone)
jone = jone + 1.0
enddo
endif
tryindex = sum(spread(sum(proxone([int(intrperm)],[int(intrperm)])*fitone,1),1,1))
if ((tryindex > (begindexone + 1.0e-008))) then
nchange = 1.0
begindexone = tryindex
outpermone = intrperm
endif
enddo
enddo
enddo
if ((kblock > 1.0)) then
do k = 2,nint(kblock)
do nlimlow = 1,nint((n+1-k))
intrperm = outpermone
do j = 1,k
intrperm(1,nlimlow+j-1) = outpermone(1,nlimlow+k-j)
enddo
tryindex = sum(spread(sum(proxone([int(intrperm)],[int(intrperm)])*fitone,1),1,1))
if ((tryindex > (begindexone + 1.0e-008))) then
nchange = 1.0
begindexone = tryindex
outpermone = intrperm
endif
enddo
enddo
endif
enddo
endif
n=size(prox,1)
acondiff = 1.0
addcon = 0.0
proxave = 0.
fit = 0.
fmadd = 0.
do while (acondiff >= 1.0e-005)
do i =1,nint(n)
do j=1,nint(n)
if ((i /= j)) then
fit(i,j) = proxone(int(outpermone(1,i)),int(outpermone(1,j))) + addcon
else
fit(i,j) = 0.0
endif
enddo
enddo
addconpv = addcon
work = 0.
cr = 1.0
do while (cr >= 1.0e-006)
cr = 0.0
index = 0.0
do jone = 1,nint((n-2))
do jtwo = (jone+1),nint((n-1))
do jthree = (jtwo+1),nint(n)
p1 = fit(jone,jtwo)
p2 = fit(jone,jthree)
p3 = fit(jtwo,jthree)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(index)+1,1)
fit(jone,jthree) = fit(jone,jthree) - work(int(index)+2,1)
fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(index)+3,1)
del = (fit(jone,jthree) - fit(jone,jtwo) -                fit(jtwo,jthree))/3.0
fit(jone,jthree) = fit(jone,jthree) - del
fit(jone,jtwo) = fit(jone,jtwo) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
work(int(index)+1,1) = del
work(int(index)+2,1) = -del
work(int(index)+3,1) = del
index = index + 3.0
cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree))                + abs(p3 &
      &- fit(jtwo,jthree))
enddo
enddo
enddo
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
p1 = fit(jone,jtwo)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(index)+1,1)
if ((fit(jone,jtwo) < 0.0)) then
work(int(index)+1,1) = -fit(jone,jtwo)
fit(jone,jtwo) = 0.0
else
work(int(index)+1,1) = 0.0
endif
index = index + 1.0
cr= cr + abs(p1-fit(jone,jtwo))
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(proxone([int(outpermone)],[int(outpermone)]) - fit,1),1,1))&
      &/(n*(n-1.0))
acondiff = abs(addcon - addconpv)
enddo
aveprox = sum(spread(sum(proxone,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((proxone([int(outpermone)],[int(outpermone)]) - (fmadd))**2.0,&
      &1),1,1))
denom = sum(spread(sum((proxone([int(outpermone)],[int(outpermone)]) - proxave)**2.0&
      &,1),1,1))
vaf = 1.0 - (diff/denom)
coord = 0.
do i = 1,nint(n)
do j = 1,nint(n)
if ((i>j)) then
coord(i,1) = coord(i,1) + fit(i,j)
endif
if ((i<j)) then
coord(i,1) = coord(i,1) - fit(i,j)
endif
enddo
coord(i,1) = coord(i,1)/n
enddo
fitone = fit
coordone = coord
vafone = vaf
addconone = addcon
fitonedim = fitone
do i = 1,nint(n)
do j = 1,nint(n)
if((i /= j)) then
proxtwo(int(outpermone(1,i)),int(outpermone(1,j))) = proxone(int(outpermone(1,i)),in&
      &t(outpermone(1,j))) -                 fitonedim(i,j) + addconone
else
proxtwo(int(outpermone(1,i)),int(outpermone(1,j))) = 0.0
endif
enddo
enddo
do i = 1,nint(n)
do j = 1,nint(n)
if ((i /= j)) then
proxtwo(int(outpermtwo(1,i)),int(outpermtwo(1,j))) =                 proxtwo(int(out&
      &permtwo(1,i)),int(outpermtwo(1,j))) + fittwodim(i,j) - addcontwo
else
endif
enddo
enddo
if ((nopt == 1.0)) then
begindextwo = sum(spread(sum(proxtwo([int(outpermtwo)],[int(outpermtwo)])*fittwo,1),&
      &1,1))
nchange = 1.0
do while (nchange == 1)
nchange=0.0
do k = 1,nint((n-1))
do j = (k+1),nint(n)
intrperm = outpermtwo
intrperm(1,k) = outpermtwo(1,j)
intrperm(1,j) = outpermtwo(1,k)
tryindex = sum(spread(sum(proxtwo([int(intrperm)],[int(intrperm)])*fittwo,1),1,1))
if ((tryindex > (begindextwo + 1.0e-008))) then
nchange = 1.0
begindextwo = tryindex
outpermtwo = intrperm
endif
enddo
enddo
do k = 1,nint(kblock)
do insertpt = 1,nint((n+1))
do nlimlow = 1,nint((n+1-k))
intrperm = outpermtwo
if ((nlimlow > insertpt)) then
jtwo = 0.0
do j = insertpt,(insertpt+k-1)
intrperm(1,j) =outpermtwo(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = (insertpt+k),(nlimlow+k-1)
intrperm(1,j) = outpermtwo(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) = outpermtwo(1,nlimlow+jtwo)
jtwo = jtwo + 1.0
enddo
jone = 0.0
do j = nlimlow,(insertpt-k-1)
intrperm(1,j) = outpermtwo(1,nlimlow+k+jone)
jone = jone + 1.0
enddo
endif
tryindex = sum(spread(sum(proxtwo([int(intrperm)],[int(intrperm)])*fittwo,1),1,1))
if ((tryindex > (begindextwo + 1.0e-008))) then
nchange = 1.0
begindextwo = tryindex
outpermtwo = intrperm
endif
enddo
enddo
enddo
if ((kblock > 1.0)) then
do k = 2,nint(kblock)
do nlimlow = 1,nint((n+1-k))
intrperm = outpermtwo
do j = 1,k
intrperm(1,nlimlow+j-1) = outpermtwo(1,nlimlow+k-j)
enddo
tryindex = sum(spread(sum(proxtwo([int(intrperm)],[int(intrperm)])*fittwo,1),1,1))
if ((tryindex > (begindextwo + 1.0e-008))) then
nchange = 1.0
begindextwo = tryindex
outpermtwo = intrperm
endif
enddo
enddo
endif
enddo
endif
n=size(prox,1)
acondiff = 1.0
addcon = 0.0
proxave = 0.
fit = 0.
fmadd = 0.
do while (acondiff >= 1.0e-005)
do i =1,nint(n)
do j=1,nint(n)
if ((i /= j)) then
fit(i,j) = proxtwo(int(outpermtwo(1,i)),int(outpermtwo(1,j))) + addcon
else
fit(i,j) = 0.0
endif
enddo
enddo
addconpv = addcon
work = 0.
cr = 1.0
do while (cr >= 1.0e-006)
cr = 0.0
index = 0.0
do jone = 1,nint((n-2))
do jtwo = (jone+1),nint((n-1))
do jthree = (jtwo+1),nint(n)
p1 = fit(jone,jtwo)
p2 = fit(jone,jthree)
p3 = fit(jtwo,jthree)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(index)+1,1)
fit(jone,jthree) = fit(jone,jthree) - work(int(index)+2,1)
fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(index)+3,1)
del = (fit(jone,jthree) - fit(jone,jtwo) -                fit(jtwo,jthree))/3.0
fit(jone,jthree) = fit(jone,jthree) - del
fit(jone,jtwo) = fit(jone,jtwo) + del
fit(jtwo,jthree) = fit(jtwo,jthree) + del
work(int(index)+1,1) = del
work(int(index)+2,1) = -del
work(int(index)+3,1) = del
index = index + 3.0
cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree))                + abs(p3 &
      &- fit(jtwo,jthree))
enddo
enddo
enddo
do jone = 1,nint((n-1))
do jtwo = (jone+1),nint(n)
p1 = fit(jone,jtwo)
fit(jone,jtwo) = fit(jone,jtwo) - work(int(index)+1,1)
if ((fit(jone,jtwo) < 0.0)) then
work(int(index)+1,1) = -fit(jone,jtwo)
fit(jone,jtwo) = 0.0
else
work(int(index)+1,1) = 0.0
endif
index = index + 1.0
cr= cr + abs(p1-fit(jone,jtwo))
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(proxtwo([int(outpermtwo)],[int(outpermtwo)]) - fit,1),1,1))&
      &/(n*(n-1.0))
acondiff = abs(addcon - addconpv)
enddo
aveprox = sum(spread(sum(proxtwo,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((proxtwo([int(outpermtwo)],[int(outpermtwo)]) - (fmadd))**2.0,&
      &1),1,1))
denom = sum(spread(sum((proxtwo([int(outpermtwo)],[int(outpermtwo)]) - proxave)**2.0&
      &,1),1,1))
vaf = 1.0 - (diff/denom)
coord = 0.
do i = 1,nint(n)
do j = 1,nint(n)
if ((i>j)) then
coord(i,1) = coord(i,1) + fit(i,j)
endif
if ((i<j)) then
coord(i,1) = coord(i,1) - fit(i,j)
endif
enddo
coord(i,1) = coord(i,1)/n
enddo
fittwo = fit
coordtwo = coord
vaftwo = vaf
addcontwo = addcon
fittwodim = fittwo
do i = 1,nint(n)
do j = 1,nint(n)
if ((i /= j)) then
proxthree(int(outpermtwo(1,i)),int(outpermtwo(1,j))) = proxtwo(int(outpermtwo(1,i)),&
      &int(outpermtwo(1,j))) -                 fittwodim(i,j) + addcontwo
else
proxthree(int(outpermtwo(1,i)),int(outpermtwo(1,j))) = 0.0
endif
enddo
enddo
diff = sum(spread(sum(proxthree**2.0,1),1,1))
denom = sum(spread(sum((prox-proxaveorig)**2.0,1),1,1))
vaf = 1.0 - (diff/denom)
vafdiff = abs(vaf-vafprev)
do i = 1,nint(n)
do j = 1,nint(n)
if ((i /= j)) then
proxone(int(outpermone(1,i)),int(outpermone(1,j))) = fitone(i,j) - addconone +      &
      &       proxthree(int(outpermone(1,i)),int(outpermone(1,j)))
endif
if ((i == j)) then
proxone(int(outpermone(1,i)),int(outpermone(1,j))) = 0.0
endif
enddo
enddo
enddo

deallocate( fitonedim, fittwodim, proxone, proxtwo, proxaveorig, proxave, proxthree,&
      & intrperm, fit, fmadd, coord)
return
end subroutine ms_biscalqa

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 03-May-2007 09:16:27     !
!------------------------------------------------------------------!
