wilcox_R: What the package does (short line)

tsreg <-
function(x,y,xout=FALSE,outfun=outpro,iter=10,varfun=pbvar,
corfun=pbcor,plotit=FALSE,WARN=TRUE,HD=FALSE,OPT=FALSE,xlab='X',ylab='Y',...){
#
#  Compute Theil-Sen regression estimator
#
#  Use Gauss-Seidel algorithm
#  when there is more than one predictor
#
#  OPT=TRUE, compute the intercept using median(y)-b_1median(X)
#  OPT=FALSE compute the intercept using median of y-b_1X
#
#  Starting with version Rallfun-v29, OPT=F is the default, which is consistent with 
#  other R functions that have been supplied for computing the Theil--Sen estimator. 
#
#
x<-as.matrix(x)
xx<-cbind(x,y)
xx<-elimna(xx)
x<-xx[,1:ncol(x)]
x<-as.matrix(x)
y<-xx[,ncol(x)+1]
temp<-NA
x<-as.matrix(x)
n=nrow(x)
n.keep=n
if(xout){
x<-as.matrix(x)
flag<-outfun(x,plotit=FALSE,...)$keep
x<-x[flag,]
y<-y[flag]
x<-as.matrix(x)
n.keep=nrow(x)
}
if(ncol(x)==1){
temp1<-tsp1reg(x,y,HD=HD,OPT=OPT)
coef<-temp1$coef
res<-temp1$res
}
if(ncol(x)>1){
for(p in 1:ncol(x)){
temp[p]<-tsp1reg(x[,p],y,OPT=OPT)$coef[2]
}
res<-y-x%*%temp
if(!HD)alpha<-median(res)
if(HD)alpha<-hd(res)
r<-matrix(NA,ncol=ncol(x),nrow=nrow(x))
tempold<-temp
for(it in 1:iter){
for(p in 1:ncol(x)){
r[,p]<-y-x%*%temp-alpha+temp[p]*x[,p]
temp[p]<-tsp1reg(x[,p],r[,p],plotit=FALSE,OPT=OPT)$coef[2]
}
if(!HD)alpha<-median(y-x%*%temp)
if(HD)alpha<-hd(y-x%*%temp)
tempold<-temp
}
coef<-c(alpha,temp)
res<-y-x%*%temp-alpha
}
yhat<-y-res
stre=NULL
temp=varfun(y)
if(temp==0){
if(WARN)print("Warning: When computing strength of association, measure of variation=0")
}
e.pow=NULL
if(temp>0){
e.pow<-varfun(yhat)/varfun(y)
if(!is.na(e.pow)){
if(e.pow>=1)e.pow<-corfun(yhat,y)$cor^2
e.pow=as.numeric(e.pow)
stre=sqrt(e.pow)
}}
if(plotit){
if(ncol(x)==1){
plot(x,y,xlab=xlab,ylab=ylab)
abline(coef)
}}
list(n=n,n.keep=n.keep,
coef=coef,residuals=res,Strength.Assoc=stre,Explanatory.Power=e.pow)
}