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) :: proxpermut_ptr,fitted_ptr,dumone_ptr real(8), allocatable :: proxpermut(:,:),fitted(:,:),dumone(:,:) integer(4) :: monproxpermut_ptr,vaf_ptr,diff_ptr real(8), allocatable :: monproxpermut(:,:),vaf(:,:),diff(:,:) ! Any other variables needed integer proxpermut_m,proxpermut_n,fitted_m,fitted_n,dumone_m,dumone_n ! CHECK FOR PROPER NUMBER OF ARGUMENTS if (nrhs .ne. 3) then print *,'ms_proxmon requires 3 input arguments';return elseif (nlhs .ne. 3) then print *,'ms_proxmon requires 3 output arguments';return endif ! Get the sizes of all the input variables proxpermut_m=mxGetm(prhs(1));proxpermut_n=mxGetn(prhs(1)) fitted_m=mxGetm(prhs(2));fitted_n=mxGetn(prhs(2)) dumone_m=mxGetm(prhs(3));dumone_n=mxGetn(prhs(3)) ! Create matrices for the return argument plhs(1)=mxCreateDoubleMatrix(proxpermut_m,proxpermut_n,0) plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0) plhs(3)=mxCreateDoubleMatrix(dumone_m,dumone_n,0) monproxpermut_ptr=mxGetPr(plhs(1)) vaf_ptr=mxGetPr(plhs(2)) diff_ptr=mxGetPr(plhs(3)) ! Copy right hand arguments to local arrays proxpermut_ptr=mxGetPr(prhs(1)) fitted_ptr=mxGetPr(prhs(2)) dumone_ptr=mxGetPr(prhs(3)) ! Allocate and copy data to arrays allocate(proxpermut(proxpermut_m,proxpermut_n),fitted(fitted_m,fitted_n),dumone(dumo& &ne_m,dumone_n)) allocate(monproxpermut(proxpermut_m,proxpermut_n),vaf(dumone_m,dumone_n),diff(dumone& &_m,dumone_n)) call mxCopyPtrToReal8(proxpermut_ptr,proxpermut,proxpermut_m*proxpermut_n) call mxCopyPtrToReal8(fitted_ptr,fitted,fitted_m*fitted_n) call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n) ! Do the actual computations in a subroutine call ms_proxmon(monproxpermut,vaf(1,1),diff(1,1),proxpermut,fitted,dumone(1,1)) call mxCopyReal8ToPtr(monproxpermut,monproxpermut_ptr,proxpermut_m*proxpermut_n) call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n) call mxCopyReal8ToPtr(diff,diff_ptr,dumone_m*dumone_n) deallocate(monproxpermut,vaf,diff,proxpermut,fitted,dumone) return contains subroutine ms_proxmon(monproxpermut,vaf,diff,proxpermut,fitted,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) proxpermut(:,:) real(8) fitted(:,:) real(8) dumone real(8) monproxpermut(:,:) real(8) vaf real(8) diff ! Matlab function pointers ! All other local variables real(8) denom integer j integer i real(8) avefit real(8) ave real(8) p2 real(8) p1 integer jfour integer jthree integer jtwo integer jone real(8) indexll real(8), allocatable :: proxave(:,:) real(8) cr real(8), allocatable :: fit(:,:) real(8), allocatable :: targ(:,:) real(8) work(756900,1) real(8) n ! Fill in vars going in and out allocate( targ(size(proxpermut,1),size(proxpermut,2)), fit(size(proxpermut,1),size(p& &roxpermut,2)), proxave(size(proxpermut,1),size(proxpermut,2))) ! --- Main computational routine. ---------------------------------! n = size(proxpermut,1) work = 0. targ = proxpermut fit = fitted cr = 1.0 proxave = 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) 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) & &.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 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 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 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 endif endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jthree,jfour)) endif indexll = indexll + 2.0 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 avefit = sum(spread(sum(fit,1),1,1))/(n*(n-1.0)) do i = 1,nint(n) do j = 1,nint(n) if(( i /= j)) then proxave(i,j) = avefit else proxave(i,j) = 0.0 endif enddo enddo diff = sum(spread(sum((fit - fitted)**2.0,1),1,1)) denom = sum(spread(sum((fit - proxave)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) monproxpermut = fit diff = (.5)*diff deallocate( targ, fit, proxave) return end subroutine ms_proxmon end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 03-May-2007 13:20:58 ! !------------------------------------------------------------------!