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,inperm_ptr,invar_ptr,addcontol_ptr,maxiteration_ptr
real(8), allocatable :: proxtm(:,:),inperm(:,:),invar(:,:),addcontol(:,:),maxiteration(&
      &:,:)
integer(4) :: fit_ptr,coord_ptr,outvar_ptr
real(8), allocatable :: fit(:,:),coord(:,:),outvar(:,:)
!     Any other variables needed
integer proxtm_m,proxtm_n,inperm_m,inperm_n,invar_m,invar_n,addcontol_m,addcontol_n,&
      &maxiteration_m,maxiteration_n
!     CHECK FOR PROPER NUMBER OF ARGUMENTS
if (nrhs .ne. 5) then
  print *,'ms_linfittmac requires 5 input arguments';return
elseif (nlhs .ne. 3) then
  print *,'ms_linfittmac requires 3 output arguments';return
endif
!     Get the sizes of all the input variables
proxtm_m=mxGetm(prhs(1));proxtm_n=mxGetn(prhs(1))
inperm_m=mxGetm(prhs(2));inperm_n=mxGetn(prhs(2))
invar_m=mxGetm(prhs(3));invar_n=mxGetn(prhs(3))
addcontol_m=mxGetm(prhs(4));addcontol_n=mxGetn(prhs(4))
maxiteration_m=mxGetm(prhs(5));maxiteration_n=mxGetn(prhs(5))
!     Create matrices for the return argument
plhs(1)=mxCreateDoubleMatrix(proxtm_m,proxtm_n,0)
plhs(2)=mxCreateDoubleMatrix(inperm_n,maxiteration_n,0)
plhs(3)=mxCreateDoubleMatrix(invar_m,invar_n,0)
fit_ptr=mxGetPr(plhs(1))
coord_ptr=mxGetPr(plhs(2))
outvar_ptr=mxGetPr(plhs(3))
!     Copy right hand arguments to local arrays
proxtm_ptr=mxGetPr(prhs(1))
inperm_ptr=mxGetPr(prhs(2))
invar_ptr=mxGetPr(prhs(3))
addcontol_ptr=mxGetPr(prhs(4))
maxiteration_ptr=mxGetPr(prhs(5))
!     Allocate and copy data to arrays
allocate(proxtm(proxtm_m,proxtm_n),inperm(inperm_m,inperm_n),invar(invar_m,invar_n),&
      &addcontol(addcontol_m,addcontol_n),maxiteration(maxiteration_m,maxiteration_n&
      &))
allocate(fit(proxtm_m,proxtm_n),coord(inperm_n,maxiteration_n),outvar(invar_m,invar_&
      &n))
call mxCopyPtrToReal8(proxtm_ptr,proxtm,proxtm_m*proxtm_n)
call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n)
call mxCopyPtrToReal8(invar_ptr,invar,invar_m*invar_n)
call mxCopyPtrToReal8(addcontol_ptr,addcontol,addcontol_m*addcontol_n)
call mxCopyPtrToReal8(maxiteration_ptr,maxiteration,maxiteration_m*maxiteration_n)
!     Do the actual computations in a subroutine
call ms_linfittmac(fit,coord,outvar,proxtm,inperm,invar,addcontol(1,1),maxiteration(&
      &1,1))
call mxCopyReal8ToPtr(fit,fit_ptr,proxtm_m*proxtm_n)
call mxCopyReal8ToPtr(coord,coord_ptr,inperm_n*maxiteration_n)
call mxCopyReal8ToPtr(outvar,outvar_ptr,invar_m*invar_n)
deallocate(fit,coord,outvar,proxtm,inperm,invar,addcontol,maxiteration)
return

contains


subroutine ms_linfittmac(fit,coord,outvar,proxtm,inperm,invar,addcontol,maxiteration)
!     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) inperm(:,:)
real(8) invar(:,:)
real(8) addcontol
real(8) maxiteration
real(8) fit(:,:)
real(8) coord(:,:)
real(8) outvar(:,:)
!     Matlab function pointers
!     All other local variables
real(8) vaf
real(8) denom
real(8) diff
real(8), allocatable :: fittedmatrix(:,:)
real(8) aveprox
real(8), allocatable :: aveproxmatrix(:,:)
real(8), allocatable :: addconmatrix(:,:)
real(8), allocatable :: colperm(:,:)
real(8), allocatable :: rowperm(:,:)
real(8) colidx
real(8) rowidx
real(8) del
real(8) p4
real(8) p3
real(8) p2
real(8) p1
integer ifour
integer ithree
integer itwo
integer ione
real(8) indexll
real(8) iteration
real(8) addconpv
real(8) work(657720,1)
real(8) cr
integer j
integer i
real(8) addcon
real(8) acondiff
integer icol
integer irow
real(8), allocatable :: mfit(:,:)
real(8), allocatable :: tempfit(:,:)
real(8) n
real(8) ncol
real(8) nrow
!     Fill in vars going in and out
allocate( tempfit(size(inperm,2),size(inperm,2)), mfit(size(inperm,2),size(inperm,2)&
      &), rowperm(size(proxtm,1),1), colperm(size(proxtm,2),1), addconmatrix(size(pr&
      &oxtm,1),size(proxtm,2)), aveproxmatrix(size(proxtm,1),size(proxtm,2)), fitted&
      &matrix(size(proxtm,1),size(proxtm,2)))
! --- Main computational routine. ---------------------------------!
nrow  = size(proxtm,1)
ncol = size(proxtm,2)
n = nrow + ncol
tempfit = 0.
fit = 0.
mfit = 0.
outvar = 0.
do irow = 1,nint(nrow)
do icol = 1,nint(ncol)
tempfit(irow,icol+int(nrow)) = proxtm(irow,icol)
tempfit(icol+int(nrow),irow) = proxtm(irow,icol)
enddo
enddo
acondiff = 1.0
addcon = 0.0
do while (acondiff >= addcontol)
do i =1,nint(n)
do j=1,nint(n)
if ((i /= j)) then
mfit(i,j) = tempfit(int(inperm(1,i)),int(inperm(1,j))) + addcon
else
mfit(i,j) = 0.0
endif
enddo
enddo
cr = 1.0
work = 0.
addconpv = addcon
iteration = 0.0
do while ((cr >= 1.0e-005) .and. (iteration < maxiteration))
iteration = iteration +1.0
cr = 0.0
indexll=0.0
do ione = 1,nint((n-3))
do itwo = (ione+1),nint((n-2))
do ithree = (itwo+1),nint((n-1))
do ifour = (ithree+1),nint(n)
if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) <= nrow))                     &
      & .and. ((inperm(1,ithree) > nrow) .and. (inperm(1,ifour)                     &
      & > nrow))) .or.                      (((inperm(1,ione) > nrow) .and. (inperm(&
      &1,itwo) > nrow))                      .and. ((inperm(1,ithree) <= nrow) .and.&
      & (inperm(1,ifour) <= nrow))))) then
p1 = mfit(ione,ithree)
p2 = mfit(ione,ifour)
p3 = mfit(itwo,ithree)
p4 = mfit(itwo,ifour)
mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+1,1)
mfit(ione,ifour) = mfit(ione,ifour) - work(int(indexll)+2,1)
mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+3,1)
mfit(itwo,ifour) = mfit(itwo,ifour) - work(int(indexll)+4,1)
del = (mfit(itwo,ithree) + mfit(ione,ifour) -                      mfit(ione,ithree)&
      & - mfit(itwo,ifour))/4.0
mfit(itwo,ithree) = mfit(itwo,ithree) - del
mfit(ione,ifour) = mfit(ione,ifour) - del
mfit(ione,ithree) = mfit(ione,ithree) + del
mfit(itwo,ifour) = mfit(itwo,ifour) + del
work(int(indexll)+1,1) = del
work(int(indexll)+2,1) = -del
work(int(indexll)+3,1) = -del
work(int(indexll)+4,1) = del
cr = cr + abs(p1 - mfit(ione,ithree)) +                      abs(p2 - mfit(ione,ifou&
      &r)) +                      abs(p3 - mfit(itwo,ithree)) +                     &
      & abs(p4 - mfit(itwo,ifour))
indexll = indexll + 4.0
endif
if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow))                   .an&
      &d. ((inperm(1,ithree) <= nrow) .and. (inperm(1,ifour)                   > nro&
      &w))) .or.                   (((inperm(1,ione) > nrow) .and. (inperm(1,itwo) <&
      &= nrow))                   .and. ((inperm(1,ithree) > nrow) .and. (inperm(1,i&
      &four) <= nrow))))) then
p1 = mfit(ione,itwo)
p2 = mfit(ione,ifour)
p3 = mfit(itwo,ithree)
p4 = mfit(ithree,ifour)
mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1)
mfit(ione,ifour) = mfit(ione,ifour) - work(int(indexll)+2,1)
mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+3,1)
mfit(ithree,ifour) = mfit(ithree,ifour) - work(int(indexll)+4,1)
del = (mfit(ione,itwo) + mfit(itwo,ithree) +                      mfit(ithree,ifour)&
      & - mfit(ione,ifour))/4.0
mfit(itwo,ithree) = mfit(itwo,ithree) - del
mfit(ione,ifour) = mfit(ione,ifour) + del
mfit(ione,itwo) = mfit(ione,itwo) - del
mfit(ithree,ifour) = mfit(ithree,ifour) - del
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = del
work(int(indexll)+3,1) = -del
work(int(indexll)+4,1) = -del
cr = cr + abs(p1 - mfit(ione,itwo)) +                      abs(p2 - mfit(ione,ifour)&
      &) +                      abs(p3 - mfit(itwo,ithree)) +                      a&
      &bs(p4 - mfit(ithree,ifour))
indexll = indexll + 4.0
endif
if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow))                   .an&
      &d. ((inperm(1,ithree) > nrow) .and. (inperm(1,ifour)                   <= nro&
      &w))) .or.                   (((inperm(1,ione) > nrow) .and. (inperm(1,itwo) <&
      &= nrow))                   .and. ((inperm(1,ithree) <= nrow) .and. (inperm(1,&
      &ifour) > nrow))))) then
p1 = mfit(ione,itwo)
p2 = mfit(ione,ithree)
p3 = mfit(itwo,ifour)
p4 = mfit(ithree,ifour)
mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1)
mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+2,1)
mfit(itwo,ifour) = mfit(itwo,ifour) - work(int(indexll)+3,1)
mfit(ithree,ifour) = mfit(ithree,ifour) - work(int(indexll)+4,1)
del = (mfit(ione,itwo) + mfit(itwo,ifour) -                      mfit(ione,ithree) -&
      & mfit(ithree,ifour))/4.0
mfit(ione,itwo) = mfit(ione,itwo) - del
mfit(itwo,ifour) = mfit(itwo,ifour) - del
mfit(ione,ithree) = mfit(ione,ithree) + del
mfit(ithree,ifour) = mfit(ithree,ifour) + del
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = del
work(int(indexll)+3,1) = -del
work(int(indexll)+4,1) = del
cr = cr + abs(p1 - mfit(ione,itwo)) +                      abs(p2 - mfit(ione,ithree&
      &)) +                      abs(p3 - mfit(itwo,ifour)) +                      a&
      &bs(p4 - mfit(ithree,ifour))
indexll = indexll + 4.0
endif
enddo
enddo
enddo
enddo
do ione = 1,nint((n-2))
do itwo = (ione+1),nint((n-1))
do ithree = (itwo+1),nint(n)
if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) <= nrow))                     &
      & .and. ((inperm(1,ithree) > nrow) )) .or.                      (((inperm(1,io&
      &ne) > nrow) .and. (inperm(1,itwo) > nrow))                      .and. ((inper&
      &m(1,ithree) <= nrow) )))) then
p1 = mfit(ione,ithree)
p2 = mfit(itwo,ithree)
mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+1,1)
mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+2,1)
if ((mfit(ione,ithree) >= mfit(itwo,ithree))) then
del = 0.0
else
del = (mfit(itwo,ithree) - mfit(ione,ithree))/2.0
endif
mfit(ione,ithree) = mfit(ione,ithree) + del
mfit(itwo,ithree) = mfit(itwo,ithree) - del
work(int(indexll)+1,1) = del
work(int(indexll)+2,1) = -del
cr = cr + abs(p1 - mfit(ione,ithree)) +                      abs(p2 - mfit(itwo,ithr&
      &ee))
indexll = indexll + 2.0
endif
if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow))                   .an&
      &d. ((inperm(1,ithree) > nrow) )) .or.                   (((inperm(1,ione) > n&
      &row) .and. (inperm(1,itwo) <= nrow))                   .and. ((inperm(1,ithre&
      &e) <= nrow) )))) then
p1 = mfit(ione,itwo)
p2 = mfit(ione,ithree)
mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1)
mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+2,1)
if ((mfit(ione,ithree) >= mfit(ione,itwo))) then
del = 0.0
else
del = (mfit(ione,itwo) - mfit(ione,ithree))/2.0
endif
mfit(ione,ithree) = mfit(ione,ithree) + del
mfit(ione,itwo) = mfit(ione,itwo) - del
work(int(indexll)+1,1) = -del
work(int(indexll)+2,1) = del
cr = cr + abs(p1 - mfit(ione,itwo)) +                      abs(p2 - mfit(ione,ithree&
      &))
indexll = indexll + 2.0
endif
enddo
enddo
enddo
do ione = 1,nint((n-1))
do itwo = (ione+1),nint(n)
if ((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow)) .or.                ((&
      &inperm(1,ione) > nrow) .and. (inperm(1,itwo) <= nrow)))) then
p1 = mfit(ione,itwo)
mfit(ione,itwo) = mfit(ione,itwo) - work( int(indexll)+1,1)
if((mfit(ione,itwo) < 0.0)) then
work(int(indexll)+1,1) = -mfit(ione,itwo)
mfit(ione,itwo) = 0.0
endif
indexll = indexll + 1.0
cr = cr + abs(p1 - mfit(ione,itwo))
endif
enddo
enddo
enddo
do i = 1,nint((n-1))
do j = (i+1),nint(n)
mfit(j,i) = mfit(i,j)
enddo
enddo
rowidx = 0.0
do i = 1,nint(n)
if ((inperm(1,i) <= nrow)) then
rowidx = rowidx + 1.0
colidx = 0.0
do j = 1,nint(n)
if ((inperm(1,j) > nrow)) then
colidx = colidx + 1.0
fit(int(rowidx),int(colidx)) = mfit(i,j)
endif
enddo
endif
enddo
rowperm = 0.
colperm = 0.
colidx = 0.0
rowidx = 0.0
do i = 1,nint(n)
if ((inperm(1,i) <= nrow)) then
rowidx = rowidx + 1.0
rowperm(int(rowidx),1) = inperm(1,i)
endif
if ((inperm(1,i) > nrow)) then
colidx = colidx + 1.0
colperm(int(colidx),1) = inperm(1,i) - nrow
endif
enddo
addcon = (-1.0)*(sum(spread(sum(proxtm([int(rowperm)],[int(colperm)]) - fit,1),1,1))&
      &/(nrow*ncol))
acondiff = abs(addcon - addconpv)
addconmatrix = 0.
aveproxmatrix = 0.
aveprox = sum(spread(sum(proxtm,1),1,1))/(nrow*ncol)
do i = 1,nint(nrow)
do j = 1,nint(ncol)
addconmatrix(i,j) = addcon
aveproxmatrix(i,j) = aveprox
enddo
enddo
fittedmatrix = fit - addconmatrix
diff = sum(spread(sum((proxtm([int(rowperm)],[int(colperm)]) - fittedmatrix)**2.0,1)&
      &,1,1))
denom = sum(spread(sum((proxtm([int(rowperm)],[int(colperm)]) - aveproxmatrix)**2.0,&
      &1),1,1))
vaf = 1.0 - (diff/denom)
enddo
coord = 0.
if ((inperm(1,1) <= nrow)) then
do i = 2,nint(n)
if((inperm(1,i) > nrow)) then
coord(i,1) = mfit(1,i)
endif
enddo
do i = 2,nint((n-1))
do j = (i+1),nint(n)
if (( (inperm(1,i) <= nrow) .and. (inperm(1,j) > nrow ) )) then
coord(i,1) = coord(j,1) - mfit(i,j)
endif
if (( (inperm(1,i) > nrow) .and. (inperm(1,j) <= nrow) )) then
coord(j,1) = coord(i,1) + mfit(i,j)
endif
enddo
enddo
endif
if ((inperm(1,1) >  nrow)) then
do i = 2,nint(n)
if((inperm(1,i) <= nrow)) then
coord(i,1) = mfit(1,i)
endif
enddo
do i = 2,nint((n-1))
do j = (i+1),nint(n)
if (( (inperm(1,i) <= nrow) .and. (inperm(1,j) > nrow ) )) then
coord(j,1) = coord(i,1) + mfit(i,j)
endif
if (( (inperm(1,i) > nrow) .and. (inperm(1,j) <= nrow) )) then
coord(i,1) = coord(j,1) - mfit(i,j)
endif
enddo
enddo
endif
outvar(1,1) = vaf
outvar(2,1) = addcon

deallocate( tempfit, mfit, rowperm, colperm, addconmatrix, aveproxmatrix, fittedmatr&
      &ix)
return
end subroutine ms_linfittmac

end subroutine mexfunction




!------------------------------------------------------------------!
!     This file generated by matlab2fmex: 18-Feb-2007 14:23:14     !
!------------------------------------------------------------------!
