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_ultrafnd requires 3 input arguments';return elseif (nlhs .ne. 2) then print *,'ms_ultrafnd 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_ultrafnd(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_ultrafnd(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 integer jfour real(8), allocatable :: fit(:,:) real(8), allocatable :: targ(:,:) real(8) ave real(8) p3 real(8) p2 real(8) p1 integer jthree integer jtwo integer jone real(8) temp integer j integer i integer ithree integer itwo integer ione real(8) indexll real(8) indxord(1,3) real(8), allocatable :: proxave(:,:) real(8) iterate real(8) cr real(8) work(756900,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-2)) do itwo = (ione+1),nint((n-1)) do ithree = (itwo+1),nint(n) indxord(1,1) = inperm(1,ione) indxord(1,2) = inperm(1,itwo) indxord(1,3) = inperm(1,ithree) do i = 1,2 do j = (i+1),3 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) p1 = find(jone,jtwo) p2 = find(jone,jthree) p3 = find(jtwo,jthree) 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(jtwo,jthree) = find(jtwo,jthree) - work(int(indexll)+3,1) endif if (((find(jone,jtwo) <= find(jone,jthree)) .and. (find(jone,jtwo)& & <= find(jtwo,jthree)))) then ave = (find(jone,jthree) + find(jtwo,jthree))/2.0 work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = ave - find(jone,jthree) work(int(indexll)+3,1) = ave - find(jtwo,jthree) find(jone,jthree) = ave find(jtwo,jthree) = ave indexll = indexll + 3.0 endif if (((find(jone,jthree) <= find(jone,jtwo)) .and. (find(jone,jthre& &e) <= find(jtwo,jthree)))) then ave = (find(jone,jtwo) + find(jtwo,jthree))/2.0 work(int(indexll)+1,1) = ave - find(jone,jtwo) work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = ave - find(jtwo,jthree) find(jone,jtwo) = ave find(jtwo,jthree) = ave indexll = indexll + 3.0 endif if (((find(jtwo,jthree) <= find(jone,jthree)) .and. (find(jtwo,jth& &ree) <= find(jone,jtwo)))) then ave = (find(jone,jthree) + find(jone,jtwo))/2.0 work(int(indexll)+1,1) = ave - find(jone,jtwo) work(int(indexll)+2,1) = ave - find(jone,jthree) work(int(indexll)+3,1) = 0.0 find(jone,jtwo) = ave find(jone,jthree) = ave indexll = indexll + 3.0 endif cr = cr + abs(p1-find(jone,jtwo)) + abs(p2-find(jone,jthree)) + abs(p& &3-find(jtwo,jthree)) 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-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(indexll)+1,1) fit(jone,jthree) = fit(jone,jthree) - work(int(indexll)+2,1) fit(jtwo,jthree) = fit(jtwo,jthree) - work(int(indexll)+3,1) if (((targ(jone,jtwo) <= targ(jone,jthree)) .and. (targ(jone,jtwo)& & <= targ(jtwo,jthree)))) then ave = (fit(jone,jthree) + fit(jtwo,jthree))/2.0 work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = ave - fit(jone,jthree) work(int(indexll)+3,1) = ave - fit(jtwo,jthree) fit(jone,jthree) = ave fit(jtwo,jthree) = ave indexll = indexll + 3.0 endif if (((targ(jone,jthree) <= targ(jone,jtwo)) .and. (targ(jone,jthre& &e) <= targ(jtwo,jthree)))) then ave = (fit(jone,jtwo) + fit(jtwo,jthree))/2.0 work(int(indexll)+1,1) = ave - fit(jone,jtwo) work(int(indexll)+2,1) = 0.0 work(int(indexll)+3,1) = ave - fit(jtwo,jthree) fit(jone,jtwo) = ave fit(jtwo,jthree) = ave indexll = indexll + 3.0 endif if (((targ(jtwo,jthree) <= targ(jone,jthree)) .and. (targ(jtwo,jth& &ree) <= targ(jone,jtwo)))) then ave = (fit(jone,jthree) + fit(jone,jtwo))/2.0 work(int(indexll)+1,1) = ave - fit(jone,jtwo) work(int(indexll)+2,1) = ave - fit(jone,jthree) work(int(indexll)+3,1) = 0.0 fit(jone,jtwo) = ave fit(jone,jthree) = ave indexll = indexll + 3.0 endif 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) do jthree = 1,nint((n-1)) do jfour = (jthree+1),nint(n) 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))& & 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 indexll = indexll + 2.0 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 indexll = indexll + 2.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 indexll = indexll + 2.0 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 indexll = indexll + 2.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 indexll = indexll + 2.0 endif endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jthree,jfour)) endif 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_ultrafnd end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 09-Feb-2007 14:52:50 ! !------------------------------------------------------------------!