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,dumone_ptr real(8), allocatable :: prox(:,:),inperm(:,:),dumone(:,:) integer(4) :: find_ptr,vaf_ptr real(8), allocatable :: find(:,:),vaf(:,:) ! Any other variables needed integer prox_m,prox_n,inperm_m,inperm_n,dumone_m,dumone_n ! CHECK FOR PROPER NUMBER OF ARGUMENTS if (nrhs .ne. 3) then print *,'ms_atreefnd requires 3 input arguments';return elseif (nlhs .ne. 2) then print *,'ms_atreefnd requires 2 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)) dumone_m=mxGetm(prhs(3));dumone_n=mxGetn(prhs(3)) ! Create matrices for the return argument plhs(1)=mxCreateDoubleMatrix(prox_m,prox_n,0) plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0) find_ptr=mxGetPr(plhs(1)) vaf_ptr=mxGetPr(plhs(2)) ! Copy right hand arguments to local arrays prox_ptr=mxGetPr(prhs(1)) inperm_ptr=mxGetPr(prhs(2)) dumone_ptr=mxGetPr(prhs(3)) ! Allocate and copy data to arrays allocate(prox(prox_m,prox_n),inperm(inperm_m,inperm_n),dumone(dumone_m,dumone_n)) allocate(find(prox_m,prox_n),vaf(dumone_m,dumone_n)) call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n) call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n) call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n) ! Do the actual computations in a subroutine call ms_atreefnd(find,vaf(1,1),prox,inperm,dumone(1,1)) call mxCopyReal8ToPtr(find,find_ptr,prox_m*prox_n) call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n) deallocate(find,vaf,prox,inperm,dumone) return contains subroutine ms_atreefnd(find,vaf,prox,inperm,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) dumone real(8) find(:,:) real(8) vaf ! Matlab function pointers ! All other local variables real(8) denom real(8) diff real(8) aveprox real(8), allocatable :: fit(:,:) real(8), allocatable :: targ(:,:) real(8) del real(8) p6 real(8) p5 real(8) p4 real(8) p3 real(8) p2 real(8) p1 integer jfour integer jthree integer jtwo integer jone real(8) temp integer j integer i integer ifour integer ithree integer itwo integer ione real(8) indexll real(8) indxord(1,4) real(8), allocatable :: proxave(:,:) real(8) iterate real(8) cr real(8) work(657720,1) real(8) n ! Fill in vars going in and out allocate( proxave(size(prox,1),size(prox,2)), targ(size(prox,1),size(prox,2)), fit(s& &ize(prox,1),size(prox,2))) ! --- Main computational routine. ---------------------------------! n = size(prox,1) work = 0. find = prox cr = 1.0 iterate = 0.0 proxave = 0. indxord = 0. do while (cr >= 1.0e-006) cr = 0.0 indexll = 0.0 iterate = iterate + 1.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) indxord(1,1) = inperm(1,ione) indxord(1,2) = inperm(1,itwo) indxord(1,3) = inperm(1,ithree) indxord(1,4) = inperm(1,ifour) do i = 1,3 do j = (i+1),4 if((indxord(1,i) > indxord(1,j))) then temp = indxord(1,i) indxord(1,i) = indxord(1,j) indxord(1,j) = temp endif enddo enddo jone = indxord(1,1) jtwo = indxord(1,2) jthree = indxord(1,3) jfour = indxord(1,4) p1 = find(jone,jtwo) p2 = find(jone,jthree) p3 = find(jone,jfour) p4 = find(jtwo,jthree) p5 = find(jtwo,jfour) p6 = find(jthree,jfour) if((iterate <= 100.0)) then find(jone,jtwo) = find(jone,jtwo) - work(int(indexll)+1,1) find(jone,jthree) = find(jone,jthree) - work(int(indexll)+2,1) find(jone,jfour) = find(jone,jfour) - work(int(indexll)+3,1) find(jtwo,jthree) = find(jtwo,jthree) - work(int(indexll)+4,1) find(jtwo,jfour) = find(jtwo,jfour) - work(int(indexll)+5,1) find(jthree,jfour) = find(jthree,jfour) - work(int(indexll)+6,1) endif if((((find(jone,jtwo) + find(jthree,jfour)) <= (find(jone,jthree) & &+ find(jtwo,jfour))) .and. ((find(jone,jtwo) + find(jthree,& &jfour)) <= (find(jone,jfour) + find(jtwo,jthree))))) then del = (find(jone,jthree) + find(jtwo,jfour) - (find(jone,jfour) + & &find(jtwo,jthree)))/4.0 work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = -del work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = -del work(int(indexll)+6,1) = 0.0 find(jone,jthree) = find(jone,jthree) - del find(jtwo,jfour) = find(jtwo,jfour) - del find(jone,jfour) = find(jone,jfour) + del find(jtwo,jthree) = find(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((find(jone,jthree) + find(jtwo,jfour)) <= (find(jone,jtwo) +& & find(jthree,jfour))) .and. ((find(jone,jthree) + find(jtwo& &,jfour)) <= (find(jone,jfour) + find(jtwo,jthree))))) then del = (find(jone,jtwo) + find(jthree,jfour) - (find(jone,jfour) + & &find(jtwo,jthree)))/4.0 work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = 0.0 work(int(indexll)+6,1) = -del find(jone,jtwo) = find(jone,jtwo) - del find(jthree,jfour) = find(jthree,jfour) - del find(jone,jfour) = find(jone,jfour) + del find(jtwo,jthree) = find(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((find(jone,jfour) + find(jtwo,jthree)) <= (find(jone,jtwo) +& & find(jthree,jfour))) .and. ((find(jone,jfour) + find(jtwo,& &jthree)) <= (find(jone,jthree) + find(jtwo,jfour))))) then del = (find(jone,jtwo) + find(jthree,jfour) - (find(jone,jthree) +& & find(jtwo,jfour)))/4.0 work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del work(int(indexll)+3,1) = 0.0 work(int(indexll)+4,1) = 0.0 work(int(indexll)+5,1) = del work(int(indexll)+6,1) = -del find(jone,jtwo) = find(jone,jtwo) - del find(jthree,jfour) = find(jthree,jfour) - del find(jone,jthree) = find(jone,jthree) + del find(jtwo,jfour) = find(jtwo,jfour) + del indexll = indexll + 6.0 endif cr = cr + abs(p1-find(jone,jtwo)) + abs(p2-find(jone,jthree)) + ab& &s(p3-find(jone,jfour)) + abs(p4-find(jtwo,jthree)) + abs(p5& &-find(jtwo,jfour)) + abs(p6-find(jthree,jfour)) enddo enddo enddo enddo enddo do jone = 1,nint((n-1)) do jtwo = (jone+1),nint(n) find(jtwo,jone) = find(jone,jtwo) enddo enddo targ = find n = size(prox,1) work = 0. fit = prox cr = 1.0 proxave = 0. do while (cr >= 1.0e-006) cr = 0.0 indexll = 0.0 do jone = 1,nint((n-3)) do jtwo = (jone+1),nint((n-2)) do jthree = (jtwo+1),nint((n-1)) do jfour = (jthree+1),nint(n) p1 = fit(jone,jtwo) p2 = fit(jone,jthree) p3 = fit(jone,jfour) p4 = fit(jtwo,jthree) p5 = fit(jtwo,jfour) p6 = fit(jthree,jfour) fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1) fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1) fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1) fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1) fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1) fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1) if ((((targ(jone,jtwo) + targ(jthree,jfour)) <= (targ(jone,jthree)& & + targ(jtwo,jfour))) .and. ((targ(jone,jtwo) + targ(jthree& &,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then del = (fit(jone,jthree) + fit(jtwo,jfour) - (fit(jone,jfour) + fit& &(jtwo,jthree)))/4.0 work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = -del work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = -del work(int(indexll)+6,1) = 0.0 fit(jone,jthree) = fit(jone,jthree) - del fit(jtwo,jfour) = fit(jtwo,jfour) - del fit(jone,jfour) = fit(jone,jfour) + del fit(jtwo,jthree) = fit(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((targ(jone,jthree) + targ(jtwo,jfour)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jthree) + targ(jtwo& &,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then del = (fit(jone,jtwo) + fit(jthree,jfour) - (fit(jone,jfour) + fit& &(jtwo,jthree)))/4.0 work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = 0.0 work(int(indexll)+6,1) = -del fit(jone,jtwo) = fit(jone,jtwo) - del fit(jthree,jfour) = fit(jthree,jfour) - del fit(jone,jfour) = fit(jone,jfour) + del fit(jtwo,jthree) = fit(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((targ(jone,jfour) + targ(jtwo,jthree)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jfour) + targ(jtwo,& &jthree)) <= (targ(jone,jthree) + targ(jtwo,jfour))))) then del = (fit(jone,jtwo) + fit(jthree,jfour) - (fit(jone,jthree) + fi& &t(jtwo,jfour)))/4.0 work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del work(int(indexll)+3,1) = 0.0 work(int(indexll)+4,1) = 0.0 work(int(indexll)+5,1) = del work(int(indexll)+6,1) = -del fit(jone,jtwo) = fit(jone,jtwo) - del fit(jthree,jfour) = fit(jthree,jfour) - del fit(jone,jthree) = fit(jone,jthree) + del fit(jtwo,jfour) = fit(jtwo,jfour) + del indexll = indexll + 6.0 endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) + abs(& &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree)) + abs(p5-fit& &(jtwo,jfour)) + abs(p6-fit(jthree,jfour)) enddo enddo enddo enddo do jone = 1,nint((n-3)) do jtwo = (jone+1),nint((n-2)) do jthree = (jtwo+1),nint((n-1)) do jfour = (jthree+1),nint(n) p1 = fit(jone,jtwo) p2 = fit(jone,jthree) p3 = fit(jone,jfour) p4 = fit(jtwo,jthree) p5 = fit(jtwo,jfour) p6 = fit(jthree,jfour) fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1) fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1) fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1) fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1) fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1) fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1) if ((((targ(jone,jtwo) + targ(jthree,jfour)) <= (targ(jone,jthree)& & + targ(jtwo,jfour))) .and. ((targ(jone,jtwo) + targ(jthree& &,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then if (((fit(jone,jtwo) + fit(jthree,jfour)) <= (fit(jone,jthree) + f& &it(jtwo,jfour)))) then del = 0.0 else del = (fit(jone,jtwo) + fit(jthree,jfour) - (fit(jone,jthree) & &+ fit(jtwo,jfour)))/4.0 endif work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del work(int(indexll)+3,1) = 0.0 work(int(indexll)+4,1) = 0.0 work(int(indexll)+5,1) = del work(int(indexll)+6,1) = -del fit(jone,jtwo) = fit(jone,jtwo) - del fit(jthree,jfour) = fit(jthree,jfour) - del fit(jone,jthree) = fit(jone,jthree) + del fit(jtwo,jfour) = fit(jtwo,jfour) + del indexll = indexll + 6.0 endif if ((((targ(jone,jthree) + targ(jtwo,jfour)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jthree) + targ(jtw& &o,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then if (((fit(jone,jthree) + fit(jtwo,jfour)) <= (fit(jone,jtwo) + fit& &(jthree,jfour)))) then del = 0.0 else del = (fit(jone,jthree) + fit(jtwo,jfour) - (fit(jone,jtwo) + & &fit(jthree,jfour)))/4.0 endif work(int(indexll)+1,1) = del work(int(indexll)+2,1) = -del work(int(indexll)+3,1) = 0.0 work(int(indexll)+4,1) = 0.0 work(int(indexll)+5,1) = -del work(int(indexll)+6,1) = del fit(jone,jthree) = fit(jone,jthree) - del fit(jtwo,jfour) = fit(jtwo,jfour) - del fit(jone,jtwo) = fit(jone,jtwo) + del fit(jthree,jfour) = fit(jthree,jfour) + del indexll = indexll + 6.0 endif if ((((targ(jone,jfour) + targ(jtwo,jthree)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jfour) + targ(jtwo,& &jthree)) <= (targ(jone,jthree) + targ(jtwo,jfour))))) then if (((fit(jone,jfour) + fit(jtwo,jthree)) <= (fit(jone,jtwo) + fit& &(jthree,jfour)))) then del = 0.0 else del = (fit(jone,jfour) + fit(jtwo,jthree) - (fit(jone,jtwo) + & &fit(jthree,jfour)))/4.0 endif work(int(indexll)+1,1) = del work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = -del work(int(indexll)+4,1) = -del work(int(indexll)+5,1) = 0.0 work(int(indexll)+6,1) = del fit(jone,jfour) = fit(jone,jfour) - del fit(jtwo,jthree) = fit(jtwo,jthree) - del fit(jone,jtwo) = fit(jone,jtwo) + del fit(jthree,jfour) = fit(jthree,jfour) + del indexll = indexll + 6.0 endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) + abs(& &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree)) + abs(p5-fit& &(jtwo,jfour)) + abs(p6-fit(jthree,jfour)) enddo enddo enddo enddo do jone = 1,nint((n-3)) do jtwo = (jone+1),nint((n-2)) do jthree = (jtwo+1),nint((n-1)) do jfour = (jthree+1),nint(n) p1 = fit(jone,jtwo) p2 = fit(jone,jthree) p3 = fit(jone,jfour) p4 = fit(jtwo,jthree) p5 = fit(jtwo,jfour) p6 = fit(jthree,jfour) fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1) fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1) fit(jone,jfour) = fit(jone,jfour) - work(int(indexll)+3,1) fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+4,1) fit(jtwo,jfour) = fit(jtwo,jfour) - work(int(indexll)+5,1) fit(jthree,jfour) = fit(jthree,jfour) - work(int(indexll)+6,1) if ((((targ(jone,jtwo) + targ(jthree,jfour)) <= (targ(jone,jthree)& & + targ(jtwo,jfour))) .and. ((targ(jone,jtwo) + targ(jthree& &,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then if(((fit(jone,jtwo) + fit(jthree,jfour)) <= (fit(jone,jfour) + fit& &(jtwo,jthree)))) then del = 0.0 else del = (fit(jone,jtwo) + fit(jthree,jfour) - (fit(jone,jfour) +& & fit(jtwo,jthree)))/4.0 endif work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = 0.0 work(int(indexll)+6,1) = -del fit(jone,jtwo) = fit(jone,jtwo) - del fit(jthree,jfour) = fit(jthree,jfour) - del fit(jone,jfour) = fit(jone,jfour) + del fit(jtwo,jthree) = fit(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((targ(jone,jthree) + targ(jtwo,jfour)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jthree) + targ(jtwo& &,jfour)) <= (targ(jone,jfour) + targ(jtwo,jthree))))) then if(((fit(jone,jthree) + fit(jtwo,jfour)) <= (fit(jone,jfour) + fit& &(jtwo,jthree)))) then del = 0.0 else del = (fit(jone,jthree) + fit(jtwo,jfour) - (fit(jone,jfour) +& & fit(jtwo,jthree)))/4.0 endif work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = -del work(int(indexll)+3,1) = del work(int(indexll)+4,1) = del work(int(indexll)+5,1) = -del work(int(indexll)+6,1) = 0.0 fit(jone,jthree) = fit(jone,jthree) - del fit(jtwo,jfour) = fit(jtwo,jfour) - del fit(jone,jfour) = fit(jone,jfour) + del fit(jtwo,jthree) = fit(jtwo,jthree) + del indexll = indexll + 6.0 endif if ((((targ(jone,jfour) + targ(jtwo,jthree)) <= (targ(jone,jtwo) +& & targ(jthree,jfour))) .and. ((targ(jone,jfour) + targ(jtwo,& &jthree)) <= (targ(jone,jthree) + targ(jtwo,jfour))))) then if(((fit(jone,jfour) + fit(jtwo,jthree)) <= (fit(jone,jthree) + fi& &t(jtwo,jfour)))) then del = 0.0 else del = (fit(jone,jfour) + fit(jtwo,jthree) - (fit(jone,jthree) & &+ fit(jtwo,jfour)))/4.0 endif work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = del work(int(indexll)+3,1) = -del work(int(indexll)+4,1) = -del work(int(indexll)+5,1) = del work(int(indexll)+6,1) = 0.0 fit(jone,jfour) = fit(jone,jfour) - del fit(jtwo,jthree) = fit(jtwo,jthree) - del fit(jone,jthree) = fit(jone,jthree) + del fit(jtwo,jfour) = fit(jtwo,jfour) + del indexll = indexll + 6.0 endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) + abs(& &p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree)) + abs(p5-fit& &(jtwo,jfour)) + abs(p6-fit(jthree,jfour)) enddo enddo enddo enddo enddo do jone = 1,nint((n-1)) do jtwo = (jone+1),nint(n) fit(jtwo,jone) = fit(jone,jtwo) enddo 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 else proxave(i,j) = 0.0 endif enddo enddo diff = sum(spread(sum((prox - fit)**2.0,1),1,1)) denom = sum(spread(sum((prox - proxave)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) find = fit 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 else proxave(i,j) = 0.0 endif enddo enddo diff = sum(spread(sum((prox - find)**2.0,1),1,1)) denom = sum(spread(sum((prox - proxave)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) deallocate( proxave, targ, fit) return end subroutine ms_atreefnd end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 03-May-2007 13:30:39 ! !------------------------------------------------------------------!