       subroutine mexFunction(nlhs,plhs,nrhs,prhs)

       integer plhs(*),prhs(*)
       integer mxCreateDoubleMatrix,mxGetPr
       integer mxGetM,mxGetN,m,n,nsize
       integer fit_pr,addcon_pr,prox_pr
       integer nlhs,nrhs
       real*8 fit(1000*1000),prox(1000*1000),addcon
             
       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)
       prox_pr = mxGetPr(prhs(1))
       fit_pr = mxGetPr(plhs(1))
       addcon_pr = mxGetPr(plhs(2))
       
       call mxCopyPtrToReal8(prox_pr,prox,nsize)
       call linfitac_for(fit,addcon,prox,n)
       call mxCopyReal8toPtr(fit,fit_pr,nsize)
       call mxCopyReal8toPtr(addcon,addcon_pr,1)
       
       return

       end



   