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,kblock_ptr real(8), allocatable :: prox(:,:),inperm(:,:),kblock(:,:) integer(4) :: find_ptr,vaf_ptr,outperm_ptr real(8), allocatable :: find(:,:),vaf(:,:),outperm(:,:) ! Any other variables needed integer prox_m,prox_n,inperm_m,inperm_n,kblock_m,kblock_n ! CHECK FOR PROPER NUMBER OF ARGUMENTS if (nrhs .ne. 3) then print *,'ms_arobfnd requires 3 input arguments';return elseif (nlhs .ne. 3) then print *,'ms_arobfnd requires 3 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)) kblock_m=mxGetm(prhs(3));kblock_n=mxGetn(prhs(3)) ! Create matrices for the return argument plhs(1)=mxCreateDoubleMatrix(prox_m,prox_n,0) plhs(2)=mxCreateDoubleMatrix(kblock_m,kblock_n,0) plhs(3)=mxCreateDoubleMatrix(inperm_m,inperm_n,0) find_ptr=mxGetPr(plhs(1)) vaf_ptr=mxGetPr(plhs(2)) outperm_ptr=mxGetPr(plhs(3)) ! Copy right hand arguments to local arrays prox_ptr=mxGetPr(prhs(1)) inperm_ptr=mxGetPr(prhs(2)) kblock_ptr=mxGetPr(prhs(3)) ! Allocate and copy data to arrays allocate(prox(prox_m,prox_n),inperm(inperm_m,inperm_n),kblock(kblock_m,kblock_n)) allocate(find(prox_m,prox_n),vaf(kblock_m,kblock_n),outperm(inperm_m,inperm_n)) call mxCopyPtrToReal8(prox_ptr,prox,prox_m*prox_n) call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n) call mxCopyPtrToReal8(kblock_ptr,kblock,kblock_m*kblock_n) ! Do the actual computations in a subroutine call ms_arobfnd(find,vaf(1,1),outperm,prox,inperm,kblock(1,1)) call mxCopyReal8ToPtr(find,find_ptr,prox_m*prox_n) call mxCopyReal8ToPtr(vaf,vaf_ptr,kblock_m*kblock_n) call mxCopyReal8ToPtr(outperm,outperm_ptr,inperm_m*inperm_n) deallocate(find,vaf,outperm,prox,inperm,kblock) return contains subroutine ms_arobfnd(find,vaf,outperm,prox,inperm,kblock) ! 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) kblock real(8) find(:,:) real(8) vaf real(8) outperm(:,:) ! Matlab function pointers ! All other local variables real(8), allocatable :: prevperm(:,:) real(8) nprevperm real(8) denom real(8) diff real(8) aveprox real(8) ave real(8) p2 real(8) p1 real(8) indexll real(8), allocatable :: proxave(:,:) real(8) cr real(8), allocatable :: fit(:,:) real(8) work(24360,1) integer jone integer jtwo integer nlimlow integer insertpt real(8) rawindex real(8) tryindex real(8), allocatable :: intrperm(:,:) integer k real(8) nchange integer iterate real(8) begindex real(8) index integer j integer i real(8), allocatable :: targ(:,:) real(8) n ! Fill in vars going in and out allocate( targ(size(prox,1),size(prox,2)), intrperm(size(inperm,1),size(inperm,2)), & &fit(size(prox,1),size(prox,2)), proxave(size(prox,1),size(prox,2)), prevperm(& &size(inperm,1),size(inperm,2))) ! --- Main computational routine. ---------------------------------! n = size(prox,1) targ = 0. do i = 1,nint(n-1) do j = (i+1),nint(n) targ(i,j) = abs(i-j) targ(j,i) = targ(i,j) enddo enddo n = size(prox,1) outperm = inperm index = 1.0 begindex = sum(spread(sum(prox([int(inperm)],[int(inperm)])*targ,1),1,1)) do iterate = 1,2 nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 1,nint((n-1)) do j = (k+1),nint(n) intrperm = outperm intrperm(1,k) = outperm(1,j) intrperm(1,j) = outperm(1,k) tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if ((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index + 1.0 endif enddo enddo enddo rawindex = begindex nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 1,nint(kblock) do insertpt = 1,nint((n+1)) do nlimlow = 1,nint((n+1-k)) intrperm = outperm if ((nlimlow > insertpt)) then jtwo = 0.0 do j = insertpt,(insertpt+k-1) intrperm(1,j) =outperm(1,nlimlow+jtwo) jtwo = jtwo + 1.0 enddo jone = 0.0 do j = (insertpt+k),(nlimlow+k-1) intrperm(1,j) = outperm(1,insertpt+jone) jone = jone + 1.0 enddo endif if (((nlimlow+k) < insertpt)) then jtwo = 0.0 do j = (insertpt-k),(insertpt-1) intrperm(1,j) = outperm(1,nlimlow+jtwo) jtwo = jtwo + 1.0 enddo jone = 0.0 do j = nlimlow,(insertpt-k-1) intrperm(1,j) = outperm(1,nlimlow+k+jone) jone = jone + 1.0 enddo endif tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index +1.0 endif enddo enddo enddo enddo rawindex = begindex nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 2,nint(kblock) do nlimlow = 1,nint((n+1-k)) intrperm = outperm do j = 1,k intrperm(1,nlimlow+j-1) = outperm(1,nlimlow+k-j) enddo tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index + 1.0 endif enddo enddo enddo rawindex = begindex nchange = 1.0 enddo inperm = outperm n = size(prox,1) work = 0. fit = prox([int(inperm)],[int(inperm)]) 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)) p1 = fit(jone,jtwo) p2 = fit(jone,jtwo+1) fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1) fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(int(indexll)+2,1) if((fit(jone,jtwo) <= fit(jone,jtwo+1))) then work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = 0.0 endif if ((fit(jone,jtwo) > fit(jone,jtwo+1))) then ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0 work(int(indexll)+1,1) = ave - fit(jone,jtwo) work(int(indexll)+2,1) = ave - fit(jone,jtwo+1) fit(jone,jtwo) = ave fit(jone,jtwo+1) = ave endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jtwo+1)) indexll = indexll + 2.0 enddo enddo do jone = 3,nint(n) do jtwo = 1,(jone-2) p1 = fit(jtwo,jone) p2 = fit(jtwo+1,jone) fit(jtwo,jone) = fit(jtwo,jone) - work(int(indexll)+1,1) fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(int(indexll)+2,1) if((fit(jtwo+1,jone) <= fit(jtwo,jone))) then work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = 0.0 endif if ((fit(jtwo+1,jone) > fit(jtwo,jone))) then ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0 work(int(indexll)+1,1) = ave - fit(jtwo,jone) work(int(indexll)+2,1) = ave - fit(jtwo+1,jone) fit(jtwo,jone) = ave fit(jtwo+1,jone) = ave endif cr = cr + abs(p1-fit(jtwo,jone)) + abs(p2-fit(jtwo+1,jone)) indexll = indexll + 2.0 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([int(inperm)],[int(inperm)]) - fit)**2.0,1),1,1)) denom = sum(spread(sum((prox([int(inperm)],[int(inperm)]) - proxave)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) nprevperm = 1.0 do while (nprevperm == 1) nprevperm = 0.0 prevperm = outperm inperm = outperm n = size(prox,1) outperm = inperm index = 1.0 targ = fit begindex = sum(spread(sum(prox([int(inperm)],[int(inperm)])*targ,1),1,1)) do iterate = 1,2 nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 1,nint((n-1)) do j = (k+1),nint(n) intrperm = outperm intrperm(1,k) = outperm(1,j) intrperm(1,j) = outperm(1,k) tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if ((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index + 1.0 endif enddo enddo enddo rawindex = begindex nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 1,nint(kblock) do insertpt = 1,nint((n+1)) do nlimlow = 1,nint((n+1-k)) intrperm = outperm if ((nlimlow > insertpt)) then jtwo = 0.0 do j = insertpt,(insertpt+k-1) intrperm(1,j) =outperm(1,nlimlow+jtwo) jtwo = jtwo + 1.0 enddo jone = 0.0 do j = (insertpt+k),(nlimlow+k-1) intrperm(1,j) = outperm(1,insertpt+jone) jone = jone + 1.0 enddo endif if (((nlimlow+k) < insertpt)) then jtwo = 0.0 do j = (insertpt-k),(insertpt-1) intrperm(1,j) = outperm(1,nlimlow+jtwo) jtwo = jtwo + 1.0 enddo jone = 0.0 do j = nlimlow,(insertpt-k-1) intrperm(1,j) = outperm(1,nlimlow+k+jone) jone = jone + 1.0 enddo endif tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index +1.0 endif enddo enddo enddo enddo rawindex = begindex nchange = 1.0 do while (nchange == 1) nchange=0.0 do k = 2,nint(kblock) do nlimlow = 1,nint((n+1-k)) intrperm = outperm do j = 1,k intrperm(1,nlimlow+j-1) = outperm(1,nlimlow+k-j) enddo tryindex = sum(spread(sum(prox([int(intrperm)],[int(intrperm)])*targ,1),1,1)) if((tryindex > (begindex + 1.0e-008))) then nchange = 1.0 begindex = tryindex outperm = intrperm index = index + 1.0 endif enddo enddo enddo rawindex = begindex nchange = 1.0 enddo inperm = outperm n = size(prox,1) work = 0. fit = prox([int(inperm)],[int(inperm)]) 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)) p1 = fit(jone,jtwo) p2 = fit(jone,jtwo+1) fit(jone,jtwo) = fit(jone,jtwo) - work(int(indexll)+1,1) fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(int(indexll)+2,1) if((fit(jone,jtwo) <= fit(jone,jtwo+1))) then work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = 0.0 endif if ((fit(jone,jtwo) > fit(jone,jtwo+1))) then ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0 work(int(indexll)+1,1) = ave - fit(jone,jtwo) work(int(indexll)+2,1) = ave - fit(jone,jtwo+1) fit(jone,jtwo) = ave fit(jone,jtwo+1) = ave endif cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jtwo+1)) indexll = indexll + 2.0 enddo enddo do jone = 3,nint(n) do jtwo = 1,(jone-2) p1 = fit(jtwo,jone) p2 = fit(jtwo+1,jone) fit(jtwo,jone) = fit(jtwo,jone) - work(int(indexll)+1,1) fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(int(indexll)+2,1) if((fit(jtwo+1,jone) <= fit(jtwo,jone))) then work(int(indexll)+1,1) = 0.0 work(int(indexll)+2,1) = 0.0 endif if ((fit(jtwo+1,jone) > fit(jtwo,jone))) then ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0 work(int(indexll)+1,1) = ave - fit(jtwo,jone) work(int(indexll)+2,1) = ave - fit(jtwo+1,jone) fit(jtwo,jone) = ave fit(jtwo+1,jone) = ave endif cr = cr + abs(p1-fit(jtwo,jone)) + abs(p2-fit(jtwo+1,jone)) indexll = indexll + 2.0 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([int(inperm)],[int(inperm)]) - fit)**2.0,1),1,1)) denom = sum(spread(sum((prox([int(inperm)],[int(inperm)]) - proxave)**2.0,1),1,1)) vaf = 1.0 - (diff/denom) if ((sum(prevperm - outperm) /= 0.0)) then nprevperm = 1.0 endif enddo find = fit deallocate( targ, intrperm, fit, proxave, prevperm) return end subroutine ms_arobfnd end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 09-Feb-2007 15:01:06 ! !------------------------------------------------------------------!