       subroutine mexFunction(nlhs,plhs,nrhs,prhs)

       integer plhs(*),prhs(*)
       integer mxCreateDoubleMatrix,mxGetPr
       integer mxGetM,mxGetN,m,n,size
       integer fit_pr,vaf_pr,coord_pr,addcon_pr,prox_pr,inperm_pr
       integer nlhs,nrhs
       real*8 fit(1000*1000),prox(1000*1000),coord(1000)
       real*8 vaf,addcon
       integer inperm(1000)

      
       m = mxGetM(prhs(1))
       n = mxGetN(prhs(1))
       nsize = m*n
       if(n.gt.1000) then
          call mexErrMsgTxt('n must be <= 1000')
       end if

       plhs(1) = mxCreateDoubleMatrix(m,n,0)
       plhs(2) = mxCreateDoubleMatrix(1,1,0)
       plhs(3) = mxCreateDoubleMatrix(m,1,0)
       plhs(4) = mxCreateDoubleMatrix(1,1,0)
       prox_pr = mxGetPr(prhs(1))
       inperm_pr = mxGetPr(prhs(2))
       fit_pr = mxGetPr(plhs(1))
       vaf_pr = mxGetPr(plhs(2))
       coord_pr = mxGetPr(plhs(3))
       addcon_pr = mxGetPr(plhs(4))

     
       call mxCopyPtrToReal8(prox_pr,prox,nsize)
       call mxCopyPtrToReal8(inperm_pr,inperm,m)
       call linfitac_for(fit,vaf,coord,addcon,prox,inperm,n)
       call mxCopyReal8toPtr(fit,fit_pr,1)
       call mxCopyReal8toPtr(vaf,vaf_pr,1)
       call mxCopyReal8ToPtr(coord,coord_pr,m)
       call mxCopyReal8ToPtr(addcon,addcon_pr,1)
    

       return

       end



   