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) :: proxpermuttm_ptr,fittedtm_ptr,dumone_ptr real(8), allocatable :: proxpermuttm(:,:),fittedtm(:,:),dumone(:,:) integer(4) :: monproxpermuttm_ptr,vaf_ptr,diff_ptr real(8), allocatable :: monproxpermuttm(:,:),vaf(:,:),diff(:,:) ! Any other variables needed integer proxpermuttm_m,proxpermuttm_n,fittedtm_m,fittedtm_n,dumone_m,dumone_n ! CHECK FOR PROPER NUMBER OF ARGUMENTS if (nrhs .ne. 3) then print *,'ms_proxmontm requires 3 input arguments';return elseif (nlhs .ne. 3) then print *,'ms_proxmontm requires 3 output arguments';return endif ! Get the sizes of all the input variables proxpermuttm_m=mxGetm(prhs(1));proxpermuttm_n=mxGetn(prhs(1)) fittedtm_m=mxGetm(prhs(2));fittedtm_n=mxGetn(prhs(2)) dumone_m=mxGetm(prhs(3));dumone_n=mxGetn(prhs(3)) ! Create matrices for the return argument plhs(1)=mxCreateDoubleMatrix(proxpermuttm_m,proxpermuttm_n,0) plhs(2)=mxCreateDoubleMatrix(dumone_m,dumone_n,0) plhs(3)=mxCreateDoubleMatrix(dumone_m,dumone_n,0) monproxpermuttm_ptr=mxGetPr(plhs(1)) vaf_ptr=mxGetPr(plhs(2)) diff_ptr=mxGetPr(plhs(3)) ! Copy right hand arguments to local arrays proxpermuttm_ptr=mxGetPr(prhs(1)) fittedtm_ptr=mxGetPr(prhs(2)) dumone_ptr=mxGetPr(prhs(3)) ! Allocate and copy data to arrays allocate(proxpermuttm(proxpermuttm_m,proxpermuttm_n),fittedtm(fittedtm_m,fittedtm_n)& &,dumone(dumone_m,dumone_n)) allocate(monproxpermuttm(proxpermuttm_m,proxpermuttm_n),vaf(dumone_m,dumone_n),diff(& &dumone_m,dumone_n)) call mxCopyPtrToReal8(proxpermuttm_ptr,proxpermuttm,proxpermuttm_m*proxpermuttm_n) call mxCopyPtrToReal8(fittedtm_ptr,fittedtm,fittedtm_m*fittedtm_n) call mxCopyPtrToReal8(dumone_ptr,dumone,dumone_m*dumone_n) ! Do the actual computations in a subroutine call ms_proxmontm(monproxpermuttm,vaf(1,1),diff(1,1),proxpermuttm,fittedtm,dumone(1,& &1)) call mxCopyReal8ToPtr(monproxpermuttm,monproxpermuttm_ptr,proxpermuttm_m*proxpermutt& &m_n) call mxCopyReal8ToPtr(vaf,vaf_ptr,dumone_m*dumone_n) call mxCopyReal8ToPtr(diff,diff_ptr,dumone_m*dumone_n) deallocate(monproxpermuttm,vaf,diff,proxpermuttm,fittedtm,dumone) return contains subroutine ms_proxmontm(monproxpermuttm,vaf,diff,proxpermuttm,fittedtm,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) proxpermuttm(:,:) real(8) fittedtm(:,:) real(8) dumone real(8) monproxpermuttm(:,:) real(8) vaf real(8) diff ! Matlab function pointers ! All other local variables real(8) denom 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) cr real(8), allocatable :: fit(:,:) real(8), allocatable :: targ(:,:) real(8) work(2433600,1) real(8) n real(8) ncol real(8) nrow ! Fill in vars going in and out allocate( targ(size(proxpermuttm,1),size(proxpermuttm,2)), fit(size(proxpermuttm,1),& &size(proxpermuttm,2))) ! --- Main computational routine. ---------------------------------! nrow = size(proxpermuttm,1) ncol = size(proxpermuttm,2) n = nrow + ncol work = 0. targ = proxpermuttm fit = fittedtm cr = 1.0 do while (cr >= 1.0e-006) cr = 0.0 indexll = 0.0 do jone = 1,nint(nrow) do jtwo = 1,nint(ncol) do jthree = 1,nint(nrow) do jfour = 1,nint(ncol) 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 avefit = sum(spread(sum(fit,1),1,1))/(nrow*ncol) diff = sum(spread(sum((fit - fittedtm)**2.0,1),1,1)) denom = sum(spread(sum((fit - avefit)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) monproxpermuttm = fit deallocate( targ, fit) return end subroutine ms_proxmontm end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 03-May-2007 13:24:32 ! !------------------------------------------------------------------!