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) :: proxtm_ptr,inperm_ptr,invar_ptr,addcontol_ptr,maxiteration_ptr real(8), allocatable :: proxtm(:,:),inperm(:,:),invar(:,:),addcontol(:,:),maxiteration(& &:,:) integer(4) :: fit_ptr,coord_ptr,outvar_ptr real(8), allocatable :: fit(:,:),coord(:,:),outvar(:,:) ! Any other variables needed integer proxtm_m,proxtm_n,inperm_m,inperm_n,invar_m,invar_n,addcontol_m,addcontol_n,& &maxiteration_m,maxiteration_n ! CHECK FOR PROPER NUMBER OF ARGUMENTS if (nrhs .ne. 5) then print *,'ms_linfittmac requires 5 input arguments';return elseif (nlhs .ne. 3) then print *,'ms_linfittmac requires 3 output arguments';return endif ! Get the sizes of all the input variables proxtm_m=mxGetm(prhs(1));proxtm_n=mxGetn(prhs(1)) inperm_m=mxGetm(prhs(2));inperm_n=mxGetn(prhs(2)) invar_m=mxGetm(prhs(3));invar_n=mxGetn(prhs(3)) addcontol_m=mxGetm(prhs(4));addcontol_n=mxGetn(prhs(4)) maxiteration_m=mxGetm(prhs(5));maxiteration_n=mxGetn(prhs(5)) ! Create matrices for the return argument plhs(1)=mxCreateDoubleMatrix(proxtm_m,proxtm_n,0) plhs(2)=mxCreateDoubleMatrix(inperm_n,maxiteration_n,0) plhs(3)=mxCreateDoubleMatrix(invar_m,invar_n,0) fit_ptr=mxGetPr(plhs(1)) coord_ptr=mxGetPr(plhs(2)) outvar_ptr=mxGetPr(plhs(3)) ! Copy right hand arguments to local arrays proxtm_ptr=mxGetPr(prhs(1)) inperm_ptr=mxGetPr(prhs(2)) invar_ptr=mxGetPr(prhs(3)) addcontol_ptr=mxGetPr(prhs(4)) maxiteration_ptr=mxGetPr(prhs(5)) ! Allocate and copy data to arrays allocate(proxtm(proxtm_m,proxtm_n),inperm(inperm_m,inperm_n),invar(invar_m,invar_n),& &addcontol(addcontol_m,addcontol_n),maxiteration(maxiteration_m,maxiteration_n& &)) allocate(fit(proxtm_m,proxtm_n),coord(inperm_n,maxiteration_n),outvar(invar_m,invar_& &n)) call mxCopyPtrToReal8(proxtm_ptr,proxtm,proxtm_m*proxtm_n) call mxCopyPtrToReal8(inperm_ptr,inperm,inperm_m*inperm_n) call mxCopyPtrToReal8(invar_ptr,invar,invar_m*invar_n) call mxCopyPtrToReal8(addcontol_ptr,addcontol,addcontol_m*addcontol_n) call mxCopyPtrToReal8(maxiteration_ptr,maxiteration,maxiteration_m*maxiteration_n) ! Do the actual computations in a subroutine call ms_linfittmac(fit,coord,outvar,proxtm,inperm,invar,addcontol(1,1),maxiteration(& &1,1)) call mxCopyReal8ToPtr(fit,fit_ptr,proxtm_m*proxtm_n) call mxCopyReal8ToPtr(coord,coord_ptr,inperm_n*maxiteration_n) call mxCopyReal8ToPtr(outvar,outvar_ptr,invar_m*invar_n) deallocate(fit,coord,outvar,proxtm,inperm,invar,addcontol,maxiteration) return contains subroutine ms_linfittmac(fit,coord,outvar,proxtm,inperm,invar,addcontol,maxiteration) ! COMPUTATIONAL SUBROUTINE ! Before anything, list what modules we will use. use mexfunctions; use mexoperators; use mexcallback ! size variables ! Input/Output local mirrors real(8) proxtm(:,:) real(8) inperm(:,:) real(8) invar(:,:) real(8) addcontol real(8) maxiteration real(8) fit(:,:) real(8) coord(:,:) real(8) outvar(:,:) ! Matlab function pointers ! All other local variables real(8) vaf real(8) denom real(8) diff real(8), allocatable :: fittedmatrix(:,:) real(8) aveprox real(8), allocatable :: aveproxmatrix(:,:) real(8), allocatable :: addconmatrix(:,:) real(8), allocatable :: colperm(:,:) real(8), allocatable :: rowperm(:,:) real(8) colidx real(8) rowidx real(8) del real(8) p4 real(8) p3 real(8) p2 real(8) p1 integer ifour integer ithree integer itwo integer ione real(8) indexll real(8) iteration real(8) addconpv real(8) work(657720,1) real(8) cr integer j integer i real(8) addcon real(8) acondiff integer icol integer irow real(8), allocatable :: mfit(:,:) real(8), allocatable :: tempfit(:,:) real(8) n real(8) ncol real(8) nrow ! Fill in vars going in and out allocate( tempfit(size(inperm,2),size(inperm,2)), mfit(size(inperm,2),size(inperm,2)& &), rowperm(size(proxtm,1),1), colperm(size(proxtm,2),1), addconmatrix(size(pr& &oxtm,1),size(proxtm,2)), aveproxmatrix(size(proxtm,1),size(proxtm,2)), fitted& &matrix(size(proxtm,1),size(proxtm,2))) ! --- Main computational routine. ---------------------------------! nrow = size(proxtm,1) ncol = size(proxtm,2) n = nrow + ncol tempfit = 0. fit = 0. mfit = 0. outvar = 0. do irow = 1,nint(nrow) do icol = 1,nint(ncol) tempfit(irow,icol+int(nrow)) = proxtm(irow,icol) tempfit(icol+int(nrow),irow) = proxtm(irow,icol) enddo enddo acondiff = 1.0 addcon = 0.0 do while (acondiff >= addcontol) do i =1,nint(n) do j=1,nint(n) if ((i /= j)) then mfit(i,j) = tempfit(int(inperm(1,i)),int(inperm(1,j))) + addcon else mfit(i,j) = 0.0 endif enddo enddo cr = 1.0 work = 0. addconpv = addcon iteration = 0.0 do while ((cr >= 1.0e-005) .and. (iteration < maxiteration)) iteration = iteration +1.0 cr = 0.0 indexll=0.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) if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) <= nrow)) & & .and. ((inperm(1,ithree) > nrow) .and. (inperm(1,ifour) & & > nrow))) .or. (((inperm(1,ione) > nrow) .and. (inperm(& &1,itwo) > nrow)) .and. ((inperm(1,ithree) <= nrow) .and.& & (inperm(1,ifour) <= nrow))))) then p1 = mfit(ione,ithree) p2 = mfit(ione,ifour) p3 = mfit(itwo,ithree) p4 = mfit(itwo,ifour) mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+1,1) mfit(ione,ifour) = mfit(ione,ifour) - work(int(indexll)+2,1) mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+3,1) mfit(itwo,ifour) = mfit(itwo,ifour) - work(int(indexll)+4,1) del = (mfit(itwo,ithree) + mfit(ione,ifour) - mfit(ione,ithree)& & - mfit(itwo,ifour))/4.0 mfit(itwo,ithree) = mfit(itwo,ithree) - del mfit(ione,ifour) = mfit(ione,ifour) - del mfit(ione,ithree) = mfit(ione,ithree) + del mfit(itwo,ifour) = mfit(itwo,ifour) + del work(int(indexll)+1,1) = del work(int(indexll)+2,1) = -del work(int(indexll)+3,1) = -del work(int(indexll)+4,1) = del cr = cr + abs(p1 - mfit(ione,ithree)) + abs(p2 - mfit(ione,ifou& &r)) + abs(p3 - mfit(itwo,ithree)) + & & abs(p4 - mfit(itwo,ifour)) indexll = indexll + 4.0 endif if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow)) .an& &d. ((inperm(1,ithree) <= nrow) .and. (inperm(1,ifour) > nro& &w))) .or. (((inperm(1,ione) > nrow) .and. (inperm(1,itwo) <& &= nrow)) .and. ((inperm(1,ithree) > nrow) .and. (inperm(1,i& &four) <= nrow))))) then p1 = mfit(ione,itwo) p2 = mfit(ione,ifour) p3 = mfit(itwo,ithree) p4 = mfit(ithree,ifour) mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1) mfit(ione,ifour) = mfit(ione,ifour) - work(int(indexll)+2,1) mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+3,1) mfit(ithree,ifour) = mfit(ithree,ifour) - work(int(indexll)+4,1) del = (mfit(ione,itwo) + mfit(itwo,ithree) + mfit(ithree,ifour)& & - mfit(ione,ifour))/4.0 mfit(itwo,ithree) = mfit(itwo,ithree) - del mfit(ione,ifour) = mfit(ione,ifour) + del mfit(ione,itwo) = mfit(ione,itwo) - del mfit(ithree,ifour) = mfit(ithree,ifour) - del work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del work(int(indexll)+3,1) = -del work(int(indexll)+4,1) = -del cr = cr + abs(p1 - mfit(ione,itwo)) + abs(p2 - mfit(ione,ifour)& &) + abs(p3 - mfit(itwo,ithree)) + a& &bs(p4 - mfit(ithree,ifour)) indexll = indexll + 4.0 endif if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow)) .an& &d. ((inperm(1,ithree) > nrow) .and. (inperm(1,ifour) <= nro& &w))) .or. (((inperm(1,ione) > nrow) .and. (inperm(1,itwo) <& &= nrow)) .and. ((inperm(1,ithree) <= nrow) .and. (inperm(1,& &ifour) > nrow))))) then p1 = mfit(ione,itwo) p2 = mfit(ione,ithree) p3 = mfit(itwo,ifour) p4 = mfit(ithree,ifour) mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1) mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+2,1) mfit(itwo,ifour) = mfit(itwo,ifour) - work(int(indexll)+3,1) mfit(ithree,ifour) = mfit(ithree,ifour) - work(int(indexll)+4,1) del = (mfit(ione,itwo) + mfit(itwo,ifour) - mfit(ione,ithree) -& & mfit(ithree,ifour))/4.0 mfit(ione,itwo) = mfit(ione,itwo) - del mfit(itwo,ifour) = mfit(itwo,ifour) - del mfit(ione,ithree) = mfit(ione,ithree) + del mfit(ithree,ifour) = mfit(ithree,ifour) + del work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del work(int(indexll)+3,1) = -del work(int(indexll)+4,1) = del cr = cr + abs(p1 - mfit(ione,itwo)) + abs(p2 - mfit(ione,ithree& &)) + abs(p3 - mfit(itwo,ifour)) + a& &bs(p4 - mfit(ithree,ifour)) indexll = indexll + 4.0 endif enddo enddo enddo enddo do ione = 1,nint((n-2)) do itwo = (ione+1),nint((n-1)) do ithree = (itwo+1),nint(n) if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) <= nrow)) & & .and. ((inperm(1,ithree) > nrow) )) .or. (((inperm(1,io& &ne) > nrow) .and. (inperm(1,itwo) > nrow)) .and. ((inper& &m(1,ithree) <= nrow) )))) then p1 = mfit(ione,ithree) p2 = mfit(itwo,ithree) mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+1,1) mfit(itwo,ithree) = mfit(itwo,ithree) - work(int(indexll)+2,1) if ((mfit(ione,ithree) >= mfit(itwo,ithree))) then del = 0.0 else del = (mfit(itwo,ithree) - mfit(ione,ithree))/2.0 endif mfit(ione,ithree) = mfit(ione,ithree) + del mfit(itwo,ithree) = mfit(itwo,ithree) - del work(int(indexll)+1,1) = del work(int(indexll)+2,1) = -del cr = cr + abs(p1 - mfit(ione,ithree)) + abs(p2 - mfit(itwo,ithr& &ee)) indexll = indexll + 2.0 endif if (((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow)) .an& &d. ((inperm(1,ithree) > nrow) )) .or. (((inperm(1,ione) > n& &row) .and. (inperm(1,itwo) <= nrow)) .and. ((inperm(1,ithre& &e) <= nrow) )))) then p1 = mfit(ione,itwo) p2 = mfit(ione,ithree) mfit(ione,itwo) = mfit(ione,itwo) - work(int(indexll)+1,1) mfit(ione,ithree) = mfit(ione,ithree) - work(int(indexll)+2,1) if ((mfit(ione,ithree) >= mfit(ione,itwo))) then del = 0.0 else del = (mfit(ione,itwo) - mfit(ione,ithree))/2.0 endif mfit(ione,ithree) = mfit(ione,ithree) + del mfit(ione,itwo) = mfit(ione,itwo) - del work(int(indexll)+1,1) = -del work(int(indexll)+2,1) = del cr = cr + abs(p1 - mfit(ione,itwo)) + abs(p2 - mfit(ione,ithree& &)) indexll = indexll + 2.0 endif enddo enddo enddo do ione = 1,nint((n-1)) do itwo = (ione+1),nint(n) if ((((inperm(1,ione) <= nrow) .and. (inperm(1,itwo) > nrow)) .or. ((& &inperm(1,ione) > nrow) .and. (inperm(1,itwo) <= nrow)))) then p1 = mfit(ione,itwo) mfit(ione,itwo) = mfit(ione,itwo) - work( int(indexll)+1,1) if((mfit(ione,itwo) < 0.0)) then work(int(indexll)+1,1) = -mfit(ione,itwo) mfit(ione,itwo) = 0.0 endif indexll = indexll + 1.0 cr = cr + abs(p1 - mfit(ione,itwo)) endif enddo enddo enddo do i = 1,nint((n-1)) do j = (i+1),nint(n) mfit(j,i) = mfit(i,j) enddo enddo rowidx = 0.0 do i = 1,nint(n) if ((inperm(1,i) <= nrow)) then rowidx = rowidx + 1.0 colidx = 0.0 do j = 1,nint(n) if ((inperm(1,j) > nrow)) then colidx = colidx + 1.0 fit(int(rowidx),int(colidx)) = mfit(i,j) endif enddo endif enddo rowperm = 0. colperm = 0. colidx = 0.0 rowidx = 0.0 do i = 1,nint(n) if ((inperm(1,i) <= nrow)) then rowidx = rowidx + 1.0 rowperm(int(rowidx),1) = inperm(1,i) endif if ((inperm(1,i) > nrow)) then colidx = colidx + 1.0 colperm(int(colidx),1) = inperm(1,i) - nrow endif enddo addcon = (-1.0)*(sum(spread(sum(proxtm([int(rowperm)],[int(colperm)]) - fit,1),1,1))& &/(nrow*ncol)) acondiff = abs(addcon - addconpv) addconmatrix = 0. aveproxmatrix = 0. aveprox = sum(spread(sum(proxtm,1),1,1))/(nrow*ncol) do i = 1,nint(nrow) do j = 1,nint(ncol) addconmatrix(i,j) = addcon aveproxmatrix(i,j) = aveprox enddo enddo fittedmatrix = fit - addconmatrix diff = sum(spread(sum((proxtm([int(rowperm)],[int(colperm)]) - fittedmatrix)**2.0,1)& &,1,1)) denom = sum(spread(sum((proxtm([int(rowperm)],[int(colperm)]) - aveproxmatrix)**2.0,& &1),1,1)) vaf = 1.0 - (diff/denom) enddo coord = 0. if ((inperm(1,1) <= nrow)) then do i = 2,nint(n) if((inperm(1,i) > nrow)) then coord(i,1) = mfit(1,i) endif enddo do i = 2,nint((n-1)) do j = (i+1),nint(n) if (( (inperm(1,i) <= nrow) .and. (inperm(1,j) > nrow ) )) then coord(i,1) = coord(j,1) - mfit(i,j) endif if (( (inperm(1,i) > nrow) .and. (inperm(1,j) <= nrow) )) then coord(j,1) = coord(i,1) + mfit(i,j) endif enddo enddo endif if ((inperm(1,1) > nrow)) then do i = 2,nint(n) if((inperm(1,i) <= nrow)) then coord(i,1) = mfit(1,i) endif enddo do i = 2,nint((n-1)) do j = (i+1),nint(n) if (( (inperm(1,i) <= nrow) .and. (inperm(1,j) > nrow ) )) then coord(j,1) = coord(i,1) + mfit(i,j) endif if (( (inperm(1,i) > nrow) .and. (inperm(1,j) <= nrow) )) then coord(i,1) = coord(j,1) - mfit(i,j) endif enddo enddo endif outvar(1,1) = vaf outvar(2,1) = addcon deallocate( tempfit, mfit, rowperm, colperm, addconmatrix, aveproxmatrix, fittedmatr& &ix) return end subroutine ms_linfittmac end subroutine mexfunction !------------------------------------------------------------------! ! This file generated by matlab2fmex: 03-May-2007 13:39:44 ! !------------------------------------------------------------------!