monPrimal<-function(y,x=1:length(y),d=diff(diag(length(y))),ups=1e-10,verbose=FALSE) {
fx<-colSums(x*t(d)); i<-is.active(fx,ups=ups)
if (verbose) 
	cat("     start","\n",
		"i: ",formatC(i,digits=3,width=3),"\n",
		"x: ",formatC(x,digits=6,width=10,format="f"),
		"\n")
repeat {
	li<-length(i)
	if (li == 0) z<-y
		else {
			di<-d[i,,drop=FALSE]; lb<--lsfit(t(di),y,intercept=FALSE)$coef
			mlb<-min(lb); z<-y+colSums(lb*di)
			}
	if (verbose) 
		cat("       try","\n",
			"z: ",formatC(z,digits=6,width=10,format="f"),
			"\n")
	fz<-colSums(z*t(d)); fmz<-min(fz)
	if (is.pos(fmz) && (li == 0)) return(z)
	if (is.pos(fmz) && is.pos(mlb)) return(z)
	if (is.pos(fmz) && is.neg(mlb)) {
			i<-i[-which.min(lb)]
			x<-z; fx<-fz
			if (verbose) 
				cat("  feasible","\n",
					"i: ",formatC(i,digits=3,width=3),"\n",
					"x: ",formatC(x,digits=6,width=10,format="f"),
					"\n")
			}
	else {
		case<-"infeasible"
		k<-which((fx>0)&(fz<0))
		alw<-max(-fz[k]/(fx[k]-fz[k]))
		x<-z+alw*(x-z)
		fx<-colSums(x*t(d))
		i<-is.active(fx,ups=ups)
		if (verbose) 
			cat("infeasible","\n",
				"i: ",formatC(i,digits=3,width=3),"\n",
				"x: ",formatC(x,digits=6,width=10,format="f"),
				"\n")
		}
	}
}

is.active<-function(f,ups=1e-10) which(abs(f) < ups)

is.pos<-function(x,ups=1e-10) x > -ups

is.neg<-function(x,ups=1e-10) x < ups