diff --git a/DESCRIPTION b/DESCRIPTION
index 3508988..d804e1b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -7,7 +7,7 @@ Description: Estimation and inference methods for models for conditional quantil
   risk are also now included. See Koenker, R. (2005) Quantile Regression, Cambridge U. Press,
   <doi:10.1017/CBO9780511754098> and Koenker, R. et al. (2017) Handbook of Quantile Regression, 
   CRC Press, <doi:10.1201/9781315120256>. 
-Version: 5.99
+Version: 5.99.1
 Authors@R: c(
     person("Roger", "Koenker",  role = c("cre","aut"), email =  "rkoenker@illinois.edu"),
     person("Stephen", "Portnoy",  role = c("ctb"), 
@@ -29,6 +29,8 @@ Authors@R: c(
 	comment = "contributions to extreme value inference code"),
     person("Ivan", "Fernandez-Val",  role = c("ctb"), 
 	comment = "contributions to extreme value inference code"),
+    person("Martin", "Maechler", role = "ctb", comment = c("tweaks (src/chlfct.f, 'tiny','Large')",
+				ORCID = "0000-0002-8685-9910")),
     person(c("Brian", "D"), "Ripley",  role = c("trl","ctb"), 
 	comment = "Initial (2001) R port from S (to my everlasting shame -- 
 	how could I have been so slow to adopt R!) and for numerous other 
@@ -42,7 +44,7 @@ License: GPL (>= 2)
 URL: https://www.r-project.org
 NeedsCompilation: yes
 VignetteBuilder: R.rsp
-Packaged: 2024-10-22 10:53:40 UTC; roger
+Packaged: 2024-11-22 12:00:57 UTC; ripley
 Author: Roger Koenker [cre, aut],
   Stephen Portnoy [ctb] (Contributions to Censored QR code),
   Pin Tian Ng [ctb] (Contributions to Sparse QR code),
@@ -56,7 +58,9 @@ Author: Roger Koenker [cre, aut],
     code),
   Ivan Fernandez-Val [ctb] (contributions to extreme value inference
     code),
+  Martin Maechler [ctb] (tweaks (src/chlfct.f, 'tiny','Large'),
+    <https://orcid.org/0000-0002-8685-9910>),
   Brian D Ripley [trl, ctb] (Initial (2001) R port from S (to my
     everlasting shame -- how could I have been so slow to adopt R!) and
     for numerous other suggestions and useful advice)
-Date/Publication: 2024-10-22 12:50:02 UTC
+Date/Publication: 2024-11-22 13:27:13 UTC
diff --git a/MD5 b/MD5
index 5fb35e4..e23c8be 100644
--- a/MD5
+++ b/MD5
@@ -1,4 +1,4 @@
-37c2d088aebec596cbbee91d3f86654f *DESCRIPTION
+d4cbfa5a48d1bc69a50a2c38f78e4972 *DESCRIPTION
 59ae171f84eb2be52ee904fe100b1291 *NAMESPACE
 6c1789eaa459e5175106e2eba12bdbb4 *R/ParetoTest.R
 12fc3ee9b6eecda826cff7ed17f25e55 *R/anova.R
@@ -13,11 +13,11 @@ d2f509eb6977017d1b77faae7de97a7c *R/lprq.R
 23abbfdfe6ab274c1ebe89aaa4d825ef *R/qrisk.R
 63c8000c1eb0c126d62a215009245e74 *R/quantreg.R
 a77c6ad37ae6093e30c9d639e27351f0 *R/rqss.R
-f7723bc3ed3a21afb5b60417ddc7baa0 *R/sfn.R
+2dfec3ad6cc67d6024558b3f9a314d55 *R/sfn.R
 8177bc8bd37dd659c3e403f2fa504928 *R/table.R
 11102c1b9076b53ebaac51cc79b90bfb *R/tools.R
 00f6440dcfd495483f57d543f4918de0 *README
-9a494536e7865311fc0c2d15e7e7a410 *build/vignette.rds
+fb1aa25e922a29f7f6e3da91681027e8 *build/vignette.rds
 56ef7db1af717ba44a86062560b034df *data/Bosco.rda
 403622377f41685bab4f044f2dac7785 *data/CobarOre.rda
 3d69ca4dff26067565dddb340c8d12a3 *data/Mammals.rda
@@ -50,12 +50,12 @@ ea9c96de3d670a1fa922585cf85789e3 *demo/predemo.R
 32f68b8e497d02b2cae7cf5f9b83ff52 *demo/rqsslasso.R
 5235d61a73c80e7eb83777d7d923d144 *demo/stack.R
 a582af0bad3d273cb076a346c906cba5 *demo/subset.R
-efd23c63e38a6579b6818e3d99d15768 *inst/ChangeLog
+de7a2746ec618801c4790b614c91acf0 *inst/ChangeLog
 c323f5a63f1fa8a44efe2d32cfab3e20 *inst/FAQ
 f2732c4ab5eb2f4ca3421239a030f0ee *inst/TODO
 3fb70411ec441d151dcae73d63a5e018 *inst/doc/crq.pdf
 48b99baf4f65ac9b9a820769d55a75f0 *inst/doc/crq.pdf.asis
-8d7109abf013b827de1f1cb7a2854a34 *inst/doc/rq.pdf
+0de68334ee7d011f06785bc4926bcdc1 *inst/doc/rq.pdf
 8faba626ac2f6a7ca45473bca52e2e67 *inst/doc/rq.pdf.asis
 64449a66c0d80474c7aefef8907b8f05 *man/Bosco.Rd
 768cf6c848f78a1862b44371c1f3e2a4 *man/CobarOre.Rd
@@ -68,7 +68,7 @@ e0161d8945e1258ff78943b06d7298de *man/MelTemp.Rd
 5e19c6a25ab4ac17d770ee7c24c981e6 *man/ParetoTest.Rd
 1c55b05cfde94862b2503523ccca1998 *man/Peirce.Rd
 114a06055c68efb323cb9dd8a7278050 *man/QTECox.Rd
-9d5c7b53fc04bbbbe5820d2d616e31b5 *man/akj.Rd
+93035639b6851f7246a2cbb98b3abdfc *man/akj.Rd
 a303a4ced70b992220472a10f79fa65c *man/anova.rq.Rd
 742ec68efdfbcf12a2a47e1adcfde3fe *man/bandwidth.rq.Rd
 36991ef697fde62deb19e355974477dd *man/barro.Rd
@@ -129,7 +129,7 @@ c42b5515ccf48527ae43449357b3dfb0 *man/rq.process.object.Rd
 47b5d67bdd5108808d62c0905ab8e58a *man/rqs.fit.Rd
 60adfb48bc8cf2d2728035eb0d19f547 *man/rqss.Rd
 1d5e4c75789a5fac11763df6b7edef03 *man/rqss.object.Rd
-7a198b2a81968dc5da8b7f59ff282714 *man/sfn.control.Rd
+5bee50c05259df5e4da9d55ab9b2fa5b *man/sfn.control.Rd
 2dfb6597245f45ee542f157105db4bc4 *man/srisk.Rd
 cf9fd105c0bd87edbcce9ec4a71780e5 *man/summary.crq.Rd
 3ec599a6f561cd6908bfa5032e59be75 *man/summary.rq.Rd
@@ -142,7 +142,7 @@ bda0ebc1c480b97ce4c9cb3b2ddaf9a9 *src/Makevars
 bb6fb4e8354bba1c16aeaf0cf5f9768e *src/bound.f
 f3a2629172dbce899d97b55d42490e5f *src/boundc.f
 1bcb651bc603316f27c0c6a38f85d70b *src/brute.f
-e897c735aa425937601d47591dbdee1f *src/chlfct.f
+90e2c1e51b9fe243b7d1d03d53694700 *src/chlfct.f
 1662015996d45d980ae598a9e243e760 *src/cholesky.f
 344362611adeaf7c79c8c5a5c344aa21 *src/combos.f
 a812e708aab3cd5e073b4cb0832a07eb *src/crqf.f
@@ -164,7 +164,7 @@ bd02eb1fbf344dd91d0e24a60a1c5269 *src/powell.f
 abef5781d8b3c202a444af9fb291072b *src/pwxy.f
 b8b6da4241edbbf057b398691772e4b0 *src/qfnb.f
 00b5c50c8caa28b22e1434a7436d0bb9 *src/qselect.f
-8e36db97b4d1a23328d5ceba8c882bf7 *src/quantreg_init.c
+11136c6acb68f84acc7813dd5f146ee7 *src/quantreg_init.c
 d2eb1038dddc50e80368118675819272 *src/ratfor/README
 df3e98bf71bf48242ffe8a2939ea98b4 *src/ratfor/boot.r
 abc8e4125098a8ec492c6d56a2e990ed *src/ratfor/brute.r
@@ -195,8 +195,8 @@ cbe8e3c270c3863ecbd3d01ba512b2c4 *src/rqfnc.f
 2e4d138de042e0dc24f5432dea05a7f6 *src/rqs.f
 8420f9a87d703ab5237d1707f25eb3a9 *src/sakj.f
 aef01c414f773d044dc72cb95c8741f4 *src/sparskit2.f
-ec5e92b756a88e077a3de86c8fb7bf86 *src/srqfn.f
-5e11cb4cad5ff02717b26234beae168f *src/srqfnc.f
+50c93379d3a1af3c8d0b5e6ec305e244 *src/srqfn.f
+1ed46340b6a16e19697572e051d160a1 *src/srqfnc.f
 f651f868d227919f2c670255f04fb137 *src/srtpai.f
 a838141113a2049f3ca6b6c5914989bc *tests/panel.R
 62bf15687c35cd5a9fa8d05d94fe28a5 *tests/rq.R
diff --git a/R/sfn.R b/R/sfn.R
index 92c2667..4fb4587 100644
--- a/R/sfn.R
+++ b/R/sfn.R
@@ -11,10 +11,11 @@
 #
 
 
-sfn.control <- function( nsubmax = NULL, tmpmax = NULL, nnzlmax = NULL, 
-    cachsz = 64, small = 1e-6, maxiter=100, warn.mesg=TRUE)
+sfn.control <- function( nsubmax = NULL, tmpmax = NULL, nnzlmax = NULL,
+    cachsz = 64, small = 1e-6, maxiter=100, tiny = 1e-30, Large = 1e128, warn.mesg=TRUE)
     list(nsubmax = nsubmax, tmpmax = tmpmax, nnzlmax = nnzlmax, cachsz = cachsz,
-        small = small, maxiter = maxiter, warn.mesg = warn.mesg)
+         small = small, maxiter = as.integer(maxiter),
+         tiny = tiny, Large = Large, warn.mesg = warn.mesg)
 
 #################################################################
 # Interface for a sparse implementation of LMS interior point method
@@ -29,34 +30,35 @@ rq.fit.sfn <- function(a,y,tau=.5, rhs = (1-tau)*c(t(a) %*% rep(1,length(y))), c
 {
 	y <- -y
 	n <- length(y)
-	m <- a@dimension[2]
+	m <- as.integer(a@dimension[2])
 	if(n != a@dimension[1])
 	     stop("Dimensions of design matrix and the response vector not compatible")
-	u <- rep(1,length=n)
+	u <- rep(1,      length=n)
 	x <- rep((1-tau),length=n)
-	nnzdmax <- nnza <- a@ia[n+1]-1
+	nnzdmax <- nnza <- a@ia[n+1L] -1L
 	iwmax <- 7*m+3
 	ao <- t(a)
 	e <- ao %*% a
-	nnzemax <- e@ia[m+1]-1
+	nnzemax <- as.integer(e@ia[m+1L]) - 1L
         ctrl <- sfn.control()
         if (!missing(control)) {
              control <- as.list(control)
              ctrl[names(control)] <- control
              }
+        if (is.null(ctrl$nsubmax)) ctrl$nsubmax <- nnzemax
+        if (is.null(ctrl$tmpmax )) ctrl$tmpmax  <- 6L * m
+        if (is.null(ctrl$nnzlmax)) ctrl$nnzlmax <- 4L * nnzdmax
 	nsubmax <- ctrl$nsubmax
-	tmpmax <- ctrl$tmpmax
+	tmpmax  <- ctrl$tmpmax
 	nnzlmax <- ctrl$nnzlmax
-        if (is.null(ctrl$nsubmax)) ctrl$nsubmax <- nsubmax <- nnzemax
-        if (is.null(ctrl$tmpmax)) ctrl$tmpmax <- tmpmax <- 6 * m
-        if (is.null(ctrl$nnzlmax)) ctrl$nnzlmax <- nnzlmax <- 4 * nnzdmax
-	wwm <- vector("numeric",3*m)
 	s <- u - x
-	b1 <- solve(e, ao %*% y, tmpmax=tmpmax,nnzlmax=nnzlmax,nsubmax=nsubmax)
+	b1 <- solve(e, ao %*% y, tmpmax=tmpmax, nnzlmax=nnzlmax, nsubmax=nsubmax)
 	r <- y - a %*% b1
-	z <- ifelse(abs(r)<ctrl$small,(r*(r>0)+ctrl$small),r*(r>0))
+	z <- ifelse(abs(r) <  ctrl$small,
+                    r*(r>0) + ctrl$small,
+                    r*(r>0))
 	w <- z - r
-	wwn <- matrix(0,n,14)
+	wwn <- matrix(0, n, 14L)
 	wwn[,1] <- r
 	wwn[,2] <- z
 	wwn[,3] <- w
@@ -64,50 +66,50 @@ rq.fit.sfn <- function(a,y,tau=.5, rhs = (1-tau)*c(t(a) %*% rep(1,length(y))), c
 		n = as.integer(n),
 		m = as.integer(m),
 		nnza = as.integer(nnza),
-		a = as.double(a@ra),
+		a  = as.double (a@ra),
 		ja = as.integer(a@ja),
 		ia = as.integer(a@ia),
-		ao = as.double(ao@ra),
+		ao  = as.double (ao@ra),
 		jao = as.integer(ao@ja),
 		iao = as.integer(ao@ia),
 		nnzdmax = as.integer(nnzdmax),
 		d = double(nnzdmax),
 		jd = integer(nnzdmax),
-		id = integer(m+1),
-		dsub = double(nnzemax+1),
-		jdsub = integer(nnzemax+1),
+		id = integer(m+1L),
+		dsub = double(nnzemax+1L),
+		jdsub = integer(nnzemax+1L),
 		nnzemax = as.integer(nnzemax),
-		e = as.double(e@ra),
+		e  = as.double (e@ra),
 		je = as.integer(e@ja),
 		ie = as.integer(e@ia),
 		nsubmax = as.integer(nsubmax),
 		lindx = integer(nsubmax),
-		xlindx = integer(m+1),
+		xlindx = integer(m+1L),
 		nnzlmax = as.integer(nnzlmax),
 		lnz = double(nnzlmax),
-		xlnz = integer(m+1),
+		xlnz = integer(m+1L),
 		iw = integer(m*5),
 		iwmax = as.integer(iwmax),
 		iwork = integer(iwmax),
-		xsuper = integer(m+1),
+		xsuper = integer(m+1L),
 		tmpmax = as.integer(tmpmax),
 		tmpvec = double(tmpmax),
-		wwm = as.double(wwm),
+		wwm = double(3*m),
 		wwn = as.double(wwn),
 		cachsz = as.integer(ctrl$cachsz),
-		level = as.integer( 8 ),
+		level = 8L,
 		x = as.double(x),
 		s = as.double(s),
 		u = as.double(u),
 		c = as.double(y),
-		sol = as.double(b1),
-		rhs = as.double(rhs),
-		small = as.double(ctrl$small),
+		rhs = as.double(rhs), # 'y'
+		sol = as.double(b1),  # 'b'
+		sm_tn_Lrg = as.double(ctrl[c("small", "tiny", "Large")]),
 		ierr = integer(1),
 		maxiter = as.integer(ctrl$maxiter),
-		time = double(7))[c("sol","ierr","maxiter","time")]
+		time = double(7) )[c("sol","ierr","maxiter","time")]
         ierr <- fit$ierr
-	if(!(ierr==0) && ctrl$warn.mesg)
+	if(ierr != 0 && ctrl$warn.mesg)
             warning(sfnMessage(ierr))
 	coefficients <- -fit$sol
 	residuals <- -y - a %*% coefficients
@@ -116,7 +118,6 @@ rq.fit.sfn <- function(a,y,tau=.5, rhs = (1-tau)*c(t(a) %*% rep(1,length(y))), c
 	     control = ctrl,
              ierr = ierr,
              it = fit$maxiter)
-
 }
 #------------------------------------------------------------------------------
 #################################################################
@@ -137,8 +138,9 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 	y <- -y
 	r <- -r
 	n1 <- length(y)
-	m <- x@dimension[2]
-	if(n1 != x@dimension[1])
+        xd <- as.integer(x@dimension)
+	m <- xd[2]
+	if(n1 != xd[1])
             stop("The design matrix A1' and response vector y are not compatible")
 	n2 <- length(r)
 	if(n2 != R@dimension[1])
@@ -149,8 +151,8 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 	x2 <- rep(1,length=n2)
 	wwm <- vector("numeric",6*m)
 	wwm[1:m] <- rhs
-	nnzx <- x@ia[x@dimension[1]+1]-1
-	nnzR <- R@ia[R@dimension[1]+1]-1
+	nnzx <- x@ia[xd[1]         +1L] -1L
+	nnzR <- R@ia[R@dimension[1]+1L] -1L
 	nnzdmax <- max(nnzx,nnzR)
 	iwmax <- 7*m+3
 	ao1 <- t(x)
@@ -158,25 +160,27 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 	e <- ao1 %*% x
 	g <- ao2 %*% R
 	h <- e + g
-	nnzemax <- e@ia[e@dimension[1]+1]-1
-	nnzgmax <- g@ia[g@dimension[1]+1]-1
-	nnzhmax <- h@ia[h@dimension[1]+1]-1
+	nnzemax <- e@ia[e@dimension[1]+1L] -1L
+	nnzgmax <- g@ia[g@dimension[1]+1L] -1L
+	nnzhmax <- h@ia[h@dimension[1]+1L] -1L
         ctrl <- sfn.control()
         if (!missing(control)) {
              control <- as.list(control)
              ctrl[names(control)] <- control
              }
+        if (is.null(ctrl$nsubmax)) ctrl$nsubmax <- nnzhmax
+        if (is.null(ctrl$tmpmax )) ctrl$tmpmax  <- 6L * m
+        if (is.null(ctrl$nnzlmax)) ctrl$nnzlmax <- 4L * nnzdmax
 	nsubmax <- ctrl$nsubmax
-	tmpmax <- ctrl$tmpmax
+	tmpmax  <- ctrl$tmpmax
 	nnzlmax <- ctrl$nnzlmax
-        if (is.null(ctrl$nsubmax)) nsubmax <- nnzhmax
-        if (is.null(ctrl$tmpmax)) tmpmax <- 6 * m
-        if (is.null(ctrl$nnzlmax)) nnzlmax <- 4 * nnzdmax
 	s <- u - x1
 	chol.o <- chol(e, tmpmax=tmpmax, nsubmax=nsubmax, nnzlmax=nnzlmax)
 	b <- backsolve(chol.o, ao1 %*% y)
 	r1 <- y - x %*% b
-	z1 <- ifelse(abs(r1) < ctrl$small, (r1*(r1>0)+ctrl$small), r1*(r1>0))
+	z1 <- ifelse(abs(r1)  <  ctrl$small,
+                     r1*(r1>0) + ctrl$small,
+                     r1*(r1>0))
 	w <- z1 - r1
 	z2 <- rep(1,n2)
 	wwn1 <- matrix(0,n1,10)
@@ -205,9 +209,9 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 		nnzdmax = as.integer(nnzdmax),
 		d = double(nnzdmax),
 		jd = integer(nnzdmax),
-		id = integer(m+1),
-		dsub = double(nnzhmax+1),
-		jdsub = integer(nnzhmax+1),
+		id = integer(m+1L),
+		dsub = double(nnzhmax+1L),
+		jdsub = integer(nnzhmax+1L),
 		nnzemax = as.integer(nnzemax),
 		e = as.double(e@ra),
 		je = as.integer(e@ja),
@@ -215,21 +219,21 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 		nnzgmax = as.integer(nnzgmax),
 		g = double(nnzgmax),
 		jg = integer(nnzgmax),
-		ig = integer(m+1),
+		ig = integer(m+1L),
 		nnzhmax = as.integer(nnzhmax),
 		h = double(nnzhmax),
 		jh = integer(nnzhmax),
-		ih = integer(m+1),
+		ih = integer(m+1L),
 		nsubmax = as.integer(nsubmax),
 		lindx = integer(nsubmax),
-		xlindx = integer(m+1),
+		xlindx = integer(m+1L),
 		nnzlmax = as.integer(nnzlmax),
 		lnz = double(nnzlmax),
-		xlnz = integer(m+1),
-		iw = integer(m*5),
+		xlnz = integer(m+1L),
+		iw = integer(m*5L),
 		iwmax = as.integer(iwmax),
 		iwork = integer(iwmax),
-		xsuper = integer(m+1),
+		xsuper = integer(m+1L),
 		tmpmax = as.integer(tmpmax),
 		tmpvec = double(tmpmax),
 		maxn1n2 = as.integer(maxn1n2),
@@ -238,22 +242,23 @@ rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
 		wwn1 = as.double(wwn1),
 		wwn2 = as.double(wwn2),
 		cachsz = as.integer(ctrl$cachsz),
-		level = as.integer( 8 ),
+		level = 8L,
 		x1 = as.double(x1),
 		x2 = as.double(x2),
-		s = as.double(s),
-		u = as.double(u),
+		s  = as.double(s),
+		u  = as.double(u),
 		c1 = as.double(y),
 		c2 = as.double(r),
+		sm_tn_Lrg = as.double(ctrl[c("small", "tiny", "Large")]),
+                ## output:
 		sol = as.double(b),
-		small = as.double(ctrl$small),
 		ierr = integer(1),
 		maxiter = as.integer(ctrl$maxiter),
-		time = double(7))[c("sol","ierr","maxiter","time")]
+		time = double(7) )[c("sol","ierr","maxiter","time")]
         ierr <- fit$ierr
 	if(ierr == 13)# stop()
             stop("Increase nnzh.factor")
-	if(!(ierr==0) && ctrl$warn.mesg)
+	if(ierr != 0 && ctrl$warn.mesg)
             warning(sfnMessage(ierr))
 	coefficients <- -fit$sol
 	residuals <- -y - x %*% coefficients
diff --git a/build/vignette.rds b/build/vignette.rds
index 2e62539..a57a79c 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/ChangeLog b/inst/ChangeLog
index 330eceb..2366880 100644
--- a/inst/ChangeLog
+++ b/inst/ChangeLog
@@ -1829,3 +1829,11 @@ version
 
 3.  changed several calls to model.matrix to use contrasts.arg = rather than
 contrasts = 
+
+
+5.100
+
+1. The new `cholesky.f` (from 5.99) needed amendments to correctly pass
+  'tiny, Large' to it e.g., such that package `cobs` tests pass.
+
+2. Most naturally, extend the auxiliary R function  `sfn.control()` accordingly.
diff --git a/inst/doc/rq.pdf b/inst/doc/rq.pdf
index 54e4c77..1dddbea 100644
Binary files a/inst/doc/rq.pdf and b/inst/doc/rq.pdf differ
diff --git a/man/akj.Rd b/man/akj.Rd
index 6aeb864..f4eb0a0 100644
--- a/man/akj.Rd
+++ b/man/akj.Rd
@@ -54,7 +54,7 @@ akj(x, z =, p =, h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0)
       main = expression("score " * hat(psi) * "'" * (x)))
 
  if(require("nor1mix")) {
-  m3 <- norMix(mu= c(-4, 0, 3), sig2 = c(1/3^2, 1, 2^2),
+  m3 <- norMix(mu= c(-4, 0, 3), sigma = c(1/3, 1, 2),
                w = c(.1,.5,.4))
   plot(m3, p.norm = FALSE)
   set.seed(11)
diff --git a/man/sfn.control.Rd b/man/sfn.control.Rd
index 48c5ca7..e3a67e3 100644
--- a/man/sfn.control.Rd
+++ b/man/sfn.control.Rd
@@ -1,36 +1,29 @@
 \name{sfn.control}
 \alias{sfn.control}
-\title{Set Control Parameters for Sparse Fitting }
+\title{Set Control Parameters for Sparse Fitting}
 \description{
-Auxiliary function for setting storage dimensions and other parameters rq.fit.sfn[c]
+  Auxiliary function for setting storage dimensions and other parameters
+  for \code{\link{rq.fit.sfn}()}, \code{\link{rq.fit.sfnc}()} and \code{\link{rqss}()}.
 }
 \usage{
-sfn.control(nsubmax = NULL, tmpmax = NULL, nnzlmax = NULL, cachsz = 64, 
-	small = 1e-06, maxiter = 100, warn.mesg = TRUE)
+sfn.control(nsubmax = NULL, tmpmax = NULL, nnzlmax = NULL, cachsz = 64,
+            small = 1e-06, maxiter = 100,
+            tiny = 1e-30, Large = 1e128,
+            warn.mesg = TRUE)
 }
 \arguments{
-  \item{nsubmax}{
-upper bound for dimension of lindx
-}
-  \item{tmpmax}{
-upper bound for dimension of tmpvec
-}
-  \item{nnzlmax}{
-upper bound for non-zero entries of L stored in lnz, including diagonal
-}
-
-  \item{cachsz}{
-size of cache in kbytes on target machine
-}
-  \item{small}{
-convergence tolerance for interior point algorithm
-}
-  \item{maxiter}{
-maximal number of interior point iterations. 
-}
-  \item{warn.mesg}{
-logical flag controlling printing of warnings. 
-}
+  \item{nsubmax}{upper bound for dimension of lindx}
+  \item{tmpmax}{upper bound for dimension of tmpvec}
+  \item{nnzlmax}{upper bound for non-zero entries of L stored in lnz, including diagonal}
+  \item{cachsz}{size of cache in kbytes on target machine}
+  \item{small}{convergence tolerance for interior point algorithm}
+  \item{maxiter}{maximal number of interior point iterations}
+  \item{tiny}{a tiny positive number; values below \code{tiny * max(diag)} are replaced by \code{Large};
+    originally was \eqn{10^{-30}} hardcoded in Fortran code.}
+  \item{Large}{a large number, practically \dQuote{Infinite} to replace
+    \code{tiny} diagonal entries in Cholesky; was \eqn{10^{128}}, hardcoded in
+    compiled code.}
+  \item{warn.mesg}{logical flag controlling printing of warnings.}
 }
 \details{
 Sparse fitting requires a number of temporary storage arrays whose size depends
@@ -39,7 +32,7 @@ these sizes and some other fitting aspects can be controlled by specifying eleme
 of this control object.
 }
 \value{
-List with components named as the arguments given above.
+  A \code{\link{list}} with components named as the arguments given above.
 }
 
 \author{
@@ -47,6 +40,8 @@ Roger Koenker
 }
 
 \seealso{
-See Also \code{\link{rq.fit.sfn}}
+  \code{\link{rq.fit.sfn}},
+  \code{\link{rq.fit.sfnc}},
+  and \code{\link{rqss}} from which \code{sfn.control()} is called.
 }
 \keyword{ utilities }
diff --git a/src/chlfct.f b/src/chlfct.f
index 3ce2896..59f38c5 100644
--- a/src/chlfct.f
+++ b/src/chlfct.f
@@ -8,12 +8,13 @@
       subroutine chlfct(m,xlindx,lindx,invp,perm,iwork,nnzdsub,jdsub,
      &                  colcnt,nsuper,snode,xsuper,nnzlmax,nsubmax,
      &                  xlnz,lnz,id,jd,d,cachsz,tmpmax,level,
-     &                  tmpvec,split,ierr,it,timewd)
+     &                  tmpvec,split,ierr,it,timewd,
+     &                  tiny,Large)
       integer m,ierr,nsub,iwsiz,nnzdsub,nnzl,nsuper,nnzlmax,nsubmax,
      &        tmpsiz,cachsz,tmpmax,level,it,
      &        xlindx(*),lindx(*),invp(*),perm(*),iwork(*),jdsub(*),
      &        colcnt(*),snode(*),xsuper(*),xlnz(*),id(*),jd(*),split(*)
-      double precision d(*),lnz(*),tmpvec(*),timewd(7)
+      double precision d(*),lnz(*),tmpvec(*),timewd(7),tiny,Large
       real gtimer,timbeg,timend
       external smxpy1,smxpy2,smxpy4,smxpy8
       external mmpy1,mmpy2,mmpy4,mmpy8
diff --git a/src/quantreg_init.c b/src/quantreg_init.c
index f773096..1db8e04 100644
--- a/src/quantreg_init.c
+++ b/src/quantreg_init.c
@@ -58,7 +58,7 @@ static const R_FortranMethodDef FortranEntries[] = {
     {"rqs",       (DL_FUNC) &F77_NAME(rqs),       15},
     {"sakj",      (DL_FUNC) &F77_NAME(sakj),      13},
     {"srqfn",     (DL_FUNC) &F77_NAME(srqfn),     45},
-    {"srqfnc",    (DL_FUNC) &F77_NAME(srqfnc),    65},
+    {"srqfnc",    (DL_FUNC) &F77_NAME(srqfnc),    65}, // 65 = MAX_ARGS -- can NOT be increased
     {"wxy",       (DL_FUNC) &F77_NAME(wxy),       18},
     {"xys",       (DL_FUNC) &F77_NAME(xys),       19},
     {NULL, NULL, 0}
diff --git a/src/srqfn.f b/src/srqfn.f
index a75929b..5da8f6e 100644
--- a/src/srqfn.f
+++ b/src/srqfn.f
@@ -1,17 +1,18 @@
 c        1         2         3         4         5         6         7
 c23456789012345678901234567890123456789012345678901234567890123456789012
-      subroutine srqfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
-     &                  dsub,jdsub,nnzemax,e,je,ie,nsubmax,lindx,xlindx,
-     &                  nnzlmax,lnz,xlnz,iw,iwmax,iwork,xsuper,tmpmax,
-     &                  tmpvec,wwm,wwn,cachsz,level,x,s,u,c,y,b,small,
-     &                  ierr,maxit,timewd)
+      subroutine srqfn(n,m,nnza, a,ja,ia, ao,jao,iao, nnzdmax, d,jd,id,  ! 13
+     &                 dsub,jdsub,nnzemax, e,je,ie,nsubmax,lindx,xlindx, ! 22
+     &                 nnzlmax,lnz,xlnz, iw,iwmax,iwork, xsuper,tmpmax,  ! 30
+     &                 tmpvec,wwm,wwn, cachsz,level, x,s,u,c, y,b,       ! 41
+     &                 sm_tn_Lrg, ierr,maxit, timewd)                    ! 45
+
       integer nnza,m,n,nnzdmax,nnzemax,iwmax,
      &        nnzlmax,nsubmax,cachsz,level,tmpmax,ierr,maxit,
      &        ja(nnza),jao(nnza),jdsub(nnzemax+1),jd(nnzdmax),
      &        ia(n+1),iao(m+1),id(m+1),lindx(nsubmax),xlindx(m+1),
      &        iw(m,5),xlnz(m+1),iwork(iwmax),xsuper(m+1),je(nnzemax),
      &        ie(m+1)
-      double precision small,
+      double precision sm_tn_Lrg(3), ! = c(small, tiny, Large)
      &                 a(nnza),ao(nnza),dsub(nnzemax+1),d(nnzdmax),
      &                 lnz(nnzlmax),c(n),y(m),wwm(m,3),tmpvec(tmpmax),
      &                 wwn(n,14),x(n),s(n),u(n),e(nnzemax),b(m)
@@ -23,36 +24,39 @@ subroutine srqfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
      &          level,x,s,u,c,y,b,wwn(1,1),wwn(1,2),wwn(1,3),
      &          wwn(1,4),nnzemax,e,je,ie,wwm(1,3),wwn(1,5),wwn(1,6),
      &          wwn(1,7),wwn(1,8),wwn(1,9),wwn(1,10),wwn(1,11),
-     &          wwn(1,12),wwn(1,13),wwn(1,14),wwm(1,1),small,ierr,
-     &          maxit,timewd)
+     &          wwn(1,12),wwn(1,13),wwn(1,14),wwm(1,1),  sm_tn_Lrg(1),
+     &                ierr,maxit, timewd, sm_tn_Lrg(2), sm_tn_Lrg(3))
       return
       end
+
       subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
      &                dsub,jdsub,nsubmax,lindx,xlindx,nnzlmax,lnz,
      &                xlnz,invp,perm,iwmax,iwork,colcnt,snode,xsuper,
      &                split,tmpmax,tmpvec,newrhs,cachsz,level,x,s,u,
      &                c,y,b,r,z,w,q,nnzemax,e,je,ie,dy,dx,ds,dz,dw,dxdz,
-     &                dsdw,xi,xinv,sinv,ww1,ww2,small,ierr,maxit,timewd)
+     &                dsdw, xi,xinv,sinv ,ww1,ww2, small,
+     &                ierr,maxit, timewd, tiny, Large)
+
 c        1         2         3         4         5         6         7
 c23456789012345678901234567890123456789012345678901234567890123456789012
-c Sparse implentation of LMS's interior point method via 
-c    Ng-Peyton's sparse Cholesky factorization for sparse 
+c Sparse implentation of LMS's interior point method via
+c    Ng-Peyton's sparse Cholesky factorization for sparse
 c    symmetric positive definite
 c INPUT:
 c     n -- the number of row in the coefficient matrix A'
 c     m -- the number of column in the coefficient matrix A'
 c     nnza -- the number of non-zero elements in A'
-c     a -- an nnza-vector of non-zero values of the design 
-c          matrix (A') stored in csr format 
+c     a -- an nnza-vector of non-zero values of the design
+c          matrix (A') stored in csr format
 c     ja -- an nnza-vector of indices of the non-zero elements of
-c           the coefficient matrix 
+c           the coefficient matrix
 c     ia -- an (n+1)-vector of pointers to the begining of each
 c           row in a and ja
 c     ao -- an nnza-vector of work space for the transpose of
-c           the design matrix stored in csr format or the 
+c           the design matrix stored in csr format or the
 c           design matrix stored in csc format
-c     jao -- an nnza-vector of work space for the indices of the 
-c            transpose of the design matrix 
+c     jao -- an nnza-vector of work space for the indices of the
+c            transpose of the design matrix
 c     iao -- an (n+1)-vector of pointers to the begining of each
 c            column in ao and jao
 c     nnzdmax -- upper bound of the non-zero elements in AA'
@@ -63,52 +67,52 @@ subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
 c     dsub -- the values of e excluding the diagonal elements
 c     jdsub -- the indices to dsub
 c     nsubmax -- upper bound of the dimension of lindx
-c     lindx -- an nsub-vector of interger which contains, in 
+c     lindx -- an nsub-vector of interger which contains, in
 c           column major order, the row subscripts of the nonzero
 c           entries in L in a compressed storage format
 c     xlindx -- an (m+1)-vector of integer of pointers for lindx
-c     nnzlmax -- the upper bound of the non-zero entries in 
+c     nnzlmax -- the upper bound of the non-zero entries in
 c                L stored in lnz, including the diagonal entries
 c     lnz -- First contains the non-zero entries of d; later
 c            contains the entries of the Cholesky factor
 c     xlnz -- column pointer for L stored in lnz
-c     invp -- an n-vector of integer of inverse permutation 
+c     invp -- an n-vector of integer of inverse permutation
 c             vector
 c     perm -- an n-vector of integer of permutation vector
 c     iw -- integer work array of length m
-c     iwmax -- upper bound of the general purpose integer 
+c     iwmax -- upper bound of the general purpose integer
 c              working storage iwork; set at 7*m+3
 c     iwork -- an iwsiz-vector of integer as work space
-c     colcnt -- array of length m, containing the number of 
+c     colcnt -- array of length m, containing the number of
 c               non-zeros in each column of the factor, including
 c               the diagonal entries
-c     snode -- array of length m for recording supernode 
+c     snode -- array of length m for recording supernode
 c              membership
 c     xsuper -- array of length m+1 containing the supernode
 c               partitioning
-c     split -- an m-vector with splitting of supernodes so that 
+c     split -- an m-vector with splitting of supernodes so that
 c              they fit into cache
 c     tmpmax -- upper bound of the dimension of tmpvec
 c     tmpvec -- a tmpmax-vector of temporary vector
-c     newrhs -- extra work vector for right-hand side and 
+c     newrhs -- extra work vector for right-hand side and
 c               solution
-c     cachsz -- size of the cache (in kilobytes) on the target 
+c     cachsz -- size of the cache (in kilobytes) on the target
 c               machine
 c     level -- level of loop unrolling while performing numerical
 c              factorization
 c     x -- an n-vector, the initial feasible solution in the primal
 c              that corresponds to the design matrix A'
-c     s -- an n-vector 
+c     s -- an n-vector
 c     u -- an n-vector of upper bound for x
 c     c -- an n-vector, usually the "negative" of
 c          the pseudo response
-c     y -- an m-vector, the initial dual solution 
+c     y -- an m-vector, the initial dual solution
 c     b -- an n-vector, the rhs of the equality constraint
 c          usually X'a = (1-tau)X'e in the rq setting
 c     r -- an n-vector of residuals
 c     z -- an n-vector of the dual slack variable
-c     w -- an n-vector 
-c     q -- an n-vector of work array containing the diagonal 
+c     w -- an n-vector
+c     q -- an n-vector of work array containing the diagonal
 c          elements of the Q^(-1) matrix
 c     nnzemax -- upper bound of the non-zero elements in AA'
 c     e -- an nnzdmax-vector containing the non-zero entries of
@@ -139,11 +143,11 @@ subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
 c       7 -- insufficient work space in iwork when calling symfct
 c       8 -- inconsistancy in input when calling symfct
 c       9 -- tmpsiz > tmpmax when calling bfinit; increase tmpmax
-c       10 -- nonpositive diagonal encountered when calling 
+c       10 -- nonpositive diagonal encountered when calling
 c            blkfct, the matrix is not positive definite
-c       11 -- insufficient work storage in tmpvec when calling 
+c       11 -- insufficient work storage in tmpvec when calling
 c            blkfct
-c       12 -- insufficient work storage in iwork when calling 
+c       12 -- insufficient work storage in iwork when calling
 c            blkfct
 c     maxit -- upper limit of the iteration; on return holds the
 c              number of iterations
@@ -167,7 +171,7 @@ subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
      &                 e(nnzemax),dy(m),dx(n),ds(n),dz(n),dw(n),
      &                 dxdz(n),dsdw(n),xinv(n),sinv(n),xi(n),
      &                 ww1(n),ww2(m)
-      double precision timewd(7)
+      double precision timewd(7), tiny,Large
       real gtimer,timbeg,timend
       external smxpy1,smxpy2,smxpy4,smxpy8
       external mmpy1,mmpy2,mmpy4,mmpy8
@@ -235,7 +239,7 @@ subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
       call chlfct(m,xlindx,lindx,invp,perm,iwork,nnzdsub,jdsub,
      &            colcnt,nsuper,snode,xsuper,nnzlmax,nsubmax,xlnz,lnz,
      &            ie,je,e,cachsz,tmpmax,level,tmpvec,split,ierr,it,
-     &            timewd)
+     &            timewd, tiny,Large)
       if (ierr .ne. 0) go to 100
 c
 c Call blkslv: Numerical solution for the new rhs stored in b
diff --git a/src/srqfnc.f b/src/srqfnc.f
index 3c7e0a7..dfbf55f 100644
--- a/src/srqfnc.f
+++ b/src/srqfnc.f
@@ -1,12 +1,12 @@
 c        1         2         3         4         5         6         7
 c23456789012345678901234567890123456789012345678901234567890123456789012
-      subroutine srqfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
-     &                  a2,ja2,ia2,ao2,jao2,iao2,nnzdmax,d,jd,id,
-     &                  dsub,jdsub,nnzemax,e,je,ie,nnzgmax,g,jg,ig,
-     &                  nnzhmax,h,jh,ih,nsubmax,lindx,xlindx,nnzlmax,
-     &                  lnz,xlnz,iw,iwmax,iwork,xsuper,tmpmax,tmpvec,
-     &                  maxn1n2,ww1,wwm,wwn1,wwn2,cachsz,level,x1,x2,
-     &                  s,u,c1,c2,y,small,ierr,maxit,timewd)
+      subroutine srqfnc(n1,m,nnza1, a1,ja1,ia1, ao1,jao1,iao1, n2,nnza2, ! 11
+     &                  a2,ja2,ia2, ao2,jao2,iao2, nnzdmax, d,jd,id,     ! 21
+     &                  dsub,jdsub, nnzemax, e,je,ie, nnzgmax, g,jg,ig,  ! 31
+     &                  nnzhmax, h,jh,ih, nsubmax, lindx,xlindx,nnzlmax, ! 39
+     &                  lnz,xlnz, iw,iwmax,iwork, xsuper, tmpmax,tmpvec, ! 47
+     &                  maxn1n2, ww1,wwm,wwn1,wwn2, cachsz, level,x1,x2, ! 56
+     &                  s,u,c1,c2, sm_tn_Lrg,   y, ierr,maxit, timewd)   ! 65 ( = MAX_ARGS !)
       integer nnza1,nnza2,m,n1,n2,nnzdmax,nnzemax,nnzgmax,nnzhmax,iwmax,
      &        nnzlmax,nsubmax,cachsz,level,tmpmax,ierr,maxit,maxn1n2,
      &        ja1(nnza1),jao1(nnza1),ja2(nnza2),jao2(nnza2),
@@ -14,7 +14,7 @@ subroutine srqfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
      &        ia2(n2+1),iao2(m+1),id(m+1),lindx(nsubmax),xlindx(m+1),
      &        iw(m,5),xlnz(m+1),iwork(iwmax),xsuper(m+1),je(nnzemax),
      &        ie(m+1),jg(nnzgmax),ig(m+1),jh(nnzhmax),ih(m+1)
-      double precision small,
+      double precision sm_tn_Lrg(3), !  := c(small, tiny, Large)
      &                 a1(nnza1),ao1(nnza1),a2(nnza2),ao2(nnza2),
      &                 dsub(nnzhmax+1),d(nnzdmax),g(nnzgmax),
      &                 h(nnzhmax),lnz(nnzlmax),c1(n1),c2(n2),y(m),
@@ -33,9 +33,11 @@ subroutine srqfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
      &          ie,nnzgmax,g,jg,ig,nnzhmax,h,jh,ih,wwm(1,4),wwn1(1,4),
      &          wwn2(1,4),wwn1(1,5),wwn1(1,6),wwn2(1,5),wwn1(1,7),
      &          wwn1(1,8),wwn2(1,6),wwn1(1,9),wwn1(1,10),wwn2(1,7),
-     &          maxn1n2,ww1,wwm(1,5),wwm(1,6),small,ierr,maxit,timewd)
+     &     maxn1n2,ww1,wwm(1,5),wwm(1,6), sm_tn_Lrg(1), ierr,maxit,
+     &     timewd, sm_tn_Lrg(2), sm_tn_Lrg(3))
       return
       end
+
 c23456789012345678901234567890123456789012345678901234567890123456789012
       subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
      &                a2,ja2,ia2,ao2,jao2,iao2,nnzdmax,d,jd,id,dsub,
@@ -47,48 +49,49 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
      &                ie,nnzgmax,g,jg,ig,nnzhmax,h,jh,ih,dy,dx1,
      &                dx2,ds,dz1,dz2,dw,
      &                dxdz1,dxdz2,dsdw,xi1,xi2,
-     &                maxn1n2,ww1,ww2,ww3,small,ierr,maxit,timewd)
+     &                maxn1n2, ww1,ww2,ww3, small, ierr,maxit, timewd,
+     &                tiny, Large)
 c        1         2         3         4         5         6         7
 c23456789012345678901234567890123456789012345678901234567890123456789012
-c Sparse implentation of LMS's interior point method via 
-c    Ng-Peyton's sparse Cholesky factorization for sparse 
+c Sparse implentation of LMS's interior point method via
+c    Ng-Peyton's sparse Cholesky factorization for sparse
 c    symmetric positive definite
 c INPUT:
 c     n1 -- the number of row in the coefficient matrix A1'
 c     m -- the number of column in the coefficient matrix A1'
 c     nnza1 -- the number of non-zero elements in A'
-c     a1 -- an nnza1-vector of non-zero values of the design 
-c          matrix (A1') stored in csr format 
+c     a1 -- an nnza1-vector of non-zero values of the design
+c          matrix (A1') stored in csr format
 c     ja1 -- an nnza1-vector of indices of the non-zero elements of
-c           the coefficient matrix 
+c           the coefficient matrix
 c     ia1 -- an (n1+1)-vector of pointers to the begining of each
 c           row in a1 and ja1
 c     ao1 -- an nnza1-vector of work space for the transpose of
-c           the design matrix stored in csr format or the 
+c           the design matrix stored in csr format or the
 c           design matrix stored in csc format
-c     jao1 -- an nnza1-vector of work space for the indices of the 
-c            transpose of the design matrix 
+c     jao1 -- an nnza1-vector of work space for the indices of the
+c            transpose of the design matrix
 c     iao1 -- an (n1+1)-vector of pointers to the begining of each
 c            column in ao1 and jao1
 c     n2 -- the number of row in the constraint matrix A2'
-c     nnza2 -- the number of non-zero elements in A2' 
-c     a2 -- an nnza2-vector of non-zero values of the contraint 
-c          matrix (A2') stored in csr format 
+c     nnza2 -- the number of non-zero elements in A2'
+c     a2 -- an nnza2-vector of non-zero values of the contraint
+c          matrix (A2') stored in csr format
 c     ja2 -- an nnza2-vector of indices of the non-zero elements of
-c           the constraint matrix 
+c           the constraint matrix
 c     ia2 -- an (n2+1)-vector of pointers to the begining of each
 c           row in a2 and ja2
 c     ao2 -- an nnza2-vector of work space for the transpose of
-c           the constraint matrix stored in csr format or the 
+c           the constraint matrix stored in csr format or the
 c           constraint matrix stored in csc format
-c     jao2 -- an nnza2-vector of work space for the indices of the 
-c            transpose of the constraint matrix 
+c     jao2 -- an nnza2-vector of work space for the indices of the
+c            transpose of the constraint matrix
 c     iao2 -- an (n2+1)-vector of pointers to the begining of each
 c            column in ao2 and jao2
 c     nnzdmax -- upper bound of the non-zero elements in A1A1'
 c     d -- an nnzdmax-vector of non-zero values used to store
-c          the transpose of the design matrix multiplied by the design 
-c          matrix (A1A1') stored in csr format; 
+c          the transpose of the design matrix multiplied by the design
+c          matrix (A1A1') stored in csr format;
 c          also used to store A1Q1^(-1) and A2Q2^(-1) later
 c     jd -- an nnzdmax-vector of indices in d
 c     id -- an (m+1)-vector of pointers to the begining of each
@@ -96,59 +99,59 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c     dsub -- the values of d excluding the diagonal elements
 c     jdsub -- the indices to dsub
 c     nsubmax -- upper bound of the dimension of lindx
-c     lindx -- an nsub-vector of interger which contains, in 
+c     lindx -- an nsub-vector of interger which contains, in
 c           column major order, the row subscripts of the nonzero
 c           entries in L in a compressed storage format
 c     xlindx -- an (m+1)-vector of integer of pointers for lindx
-c     nnzlmax -- the upper bound of the non-zero entries in 
+c     nnzlmax -- the upper bound of the non-zero entries in
 c                L stored in lnz, including the diagonal entries
 c     lnz -- First contains the non-zero entries of d; later
 c            contains the entries of the Cholesky factor
 c     xlnz -- column pointer for L stored in lnz
-c     invp -- an n1-vector of integer of inverse permutation 
+c     invp -- an n1-vector of integer of inverse permutation
 c             vector
 c     perm -- an n1-vector of integer of permutation vector
 c     iw -- integer work array of length m
-c     iwmax -- upper bound of the general purpose integer 
+c     iwmax -- upper bound of the general purpose integer
 c              working storage iwork; set at 7*m+3
 c     iwork -- an iwsiz-vector of integer as work space
-c     colcnt -- array of length m, containing the number of 
+c     colcnt -- array of length m, containing the number of
 c               non-zeros in each column of the factor, including
 c               the diagonal entries
-c     snode -- array of length m for recording supernode 
+c     snode -- array of length m for recording supernode
 c              membership
 c     xsuper -- array of length m+1 containing the supernode
 c               partitioning
-c     split -- an m-vector with splitting of supernodes so that 
+c     split -- an m-vector with splitting of supernodes so that
 c              they fit into cache
 c     tmpmax -- upper bound of the dimension of tmpvec
 c     tmpvec -- a tmpmax-vector of temporary vector
 c     rhs -- m-vector to store the rhs
-c     newrhs -- extra work vector for right-hand side and 
+c     newrhs -- extra work vector for right-hand side and
 c               solution
-c     cachsz -- size of the cache (in kilobytes) on the target 
+c     cachsz -- size of the cache (in kilobytes) on the target
 c               machine
 c     level -- level of loop unrolling while performing numerical
 c              factorization
-c     x1 -- an n1-vector, the initial feasible solution for the primal 
+c     x1 -- an n1-vector, the initial feasible solution for the primal
 c           solution that corresponds to the design matrix A1'
-c     x2 -- an n2-vector, the initial feasible solution for the primal 
+c     x2 -- an n2-vector, the initial feasible solution for the primal
 c           solution that corresponds to the constraint matrix A2'
-c     s -- an n1-vector 
+c     s -- an n1-vector
 c     u -- an n1-vector of the upper bound for x1
-c     c1 -- an n1-vector in the primal; negative response in the 
+c     c1 -- an n1-vector in the primal; negative response in the
 c           regression quantile setting
 c     c2 -- an n2-vector, the negative rhs of the inequality constraint
-c     y -- an m-vector, the initial dual solution 
+c     y -- an m-vector, the initial dual solution
 c     b -- an n1-vector, usualy the rhs of the equality constraint
 c          X'a = (1-tau)X'e in the rq setting
 c     r2 -- an n2-vector of residuals
 c     z1 -- an n1-vector of the dual slack variable
-c     z2 -- an n2-vector 
-c     w -- an n-vector 
-c     q1 -- an n1-vector of work array containing the diagonal 
+c     z2 -- an n2-vector
+c     w -- an n-vector
+c     q1 -- an n1-vector of work array containing the diagonal
 c          elements of the Q1^(-1) matrix
-c     q2 -- an n2-vector of work array containing the diagonal 
+c     q2 -- an n2-vector of work array containing the diagonal
 c          elements of the Q2^(-1) matrix
 c     e -- an nnzdmax-vector containing the non-zero entries of
 c          A1Q1^(-1)A1' stored in csr format
@@ -196,11 +199,11 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c       7 -- insufficient work space in iwork when calling symfct
 c       8 -- inconsistancy in input when calling symfct
 c       9 -- tmpsiz > tmpmax when calling symfct; increase tmpmax
-c       10 -- nonpositive diagonal encountered when calling 
+c       10 -- nonpositive diagonal encountered when calling
 c            blkfct
-c       11 -- insufficient work storage in tmpvec when calling 
+c       11 -- insufficient work storage in tmpvec when calling
 c            blkfct
-c       12 -- insufficient work storage in iwork when calling 
+c       12 -- insufficient work storage in iwork when calling
 c            blkfct
 c       13 -- nnzd > nnzdmax in e,je when calling amub
 c       14 -- nnzd > nnzdmax in g,jg when calling amub
@@ -208,6 +211,10 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c     maxit -- upper limit of the iteration; on return holds the
 c              number of iterations
 c     timewd -- amount of time to execute this subroutine
+c     tiny  -- tiny number; values below tiny * max(diag) are replaced by 'Large';
+c               was 10^{-30} hardcoded
+c     Large -- Large number ("Infinite") to replace tiny diagonal entries in Cholesky;
+c               was 10^{128}
 c OUTPUT:
 c     y -- an m-vector of primal solution
 c        1         2         3         4         5         6         7
@@ -232,7 +239,7 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
      &                 ds(n1),dz1(n1),dz2(n2),dw(n1),dxdz1(n1),
      &                 dxdz2(n2),dsdw(n1),
      &                 xi1(n1),xi2(n2),ww1(maxn1n2),ww2(m),ww3(m)
-      double precision timewd(7)
+      double precision timewd(7), tiny,Large
       real gtimer,timbeg,timend
       external smxpy1,smxpy2,smxpy4,smxpy8
       external mmpy1,mmpy2,mmpy4,mmpy8
@@ -257,7 +264,7 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c  and store the residuals in r1 in ds, and r3 in dy temporarily,
 c  and r2 in r2 permanently
 c
-c    Call amux to obtain A1x1 and store the value in ww2 
+c    Call amux to obtain A1x1 and store the value in ww2
 c
          call amux(m,x1,ww2,ao1,jao1,iao1)
 c
@@ -283,13 +290,13 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c  Obtain AQA = A1Q1^(-1)A1' + A2Q2^(-1)A2' in 5 steps
 c
 c  Step1: Obtain A1Q1^(-1) and store the values in d,jd,id in csr format
-c         Also compute A1Q1^(-1)r1 and store the values in ww2 to be used 
+c         Also compute A1Q1^(-1)r1 and store the values in ww2 to be used
 c         to generate r3;
-c  Step2: Compute A1Q1^(-1)A1' and store the values in e,je,ie 
+c  Step2: Compute A1Q1^(-1)A1' and store the values in e,je,ie
 c  Step3: Obtain A2Q2^(-1) and store the values in d,jd,id in csr format
-c         Also compute A2Q2^(-1)r2 and store the values in in ww3 to 
+c         Also compute A2Q2^(-1)r2 and store the values in in ww3 to
 c         be used to generate r3;
-c  Step4: Compute A2Q2^(-1)A2' and store the value in g,jg,ig 
+c  Step4: Compute A2Q2^(-1)A2' and store the value in g,jg,ig
 c  Step5: Compute AQA and store the values in h,jh,ih
 c
 c    Step 1
@@ -349,7 +356,8 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
          call chlfct(m,xlindx,lindx,invp,perm,iwork,nnzdsub,jdsub,
      &            colcnt,nsuper,snode,xsuper,nnzlmax,nsubmax,xlnz,lnz,
      &            ih,jh,h,cachsz,tmpmax,level,tmpvec,split,ierr,it,
-     &            timewd)
+     &            timewd, tiny,Large)
+
          if (ierr .ne. 0) go to 100
 c
 c Call blkslv: Numerical solution for the new rhs stored in rhs
@@ -388,16 +396,16 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
 c
 c Update mu
 c
-            mu = ddot(n1,z1,1,x1,1) + ddot(n2,z2,1,x2,1) 
+            mu = ddot(n1,z1,1,x1,1) + ddot(n2,z2,1,x2,1)
      &           + ddot(n1,w,1,s,1)
             g1 = mu + deltap*ddot(n1,z1,1,dx1,1)
-     &          + deltad*ddot(n1,dz1,1,x1,1)  
+     &          + deltad*ddot(n1,dz1,1,x1,1)
      &          + deltad*deltap*ddot(n1,dz1,1,dx1,1)
      &          + deltap*ddot(n2,z2,1,dx2,1)
-     &          + deltad*ddot(n2,dz2,1,x2,1)  
+     &          + deltad*ddot(n2,dz2,1,x2,1)
      &          + deltad*deltap*ddot(n2,dz2,1,dx2,1)
      &          + deltap*ddot(n1,w,1,ds,1)
-     &          + deltad*ddot(n1,dw,1,s,1) 
+     &          + deltad*ddot(n1,dw,1,s,1)
      &          + deltad*deltap*ddot(n1,dw,1,ds,1)
             mu = mu*((g1/mu)**3)/(2.d0*dble(n1)+dble(n2))
 c
@@ -454,13 +462,13 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
             do i = 1,n1
                dx1(i) = q1(i) * (dx1(i) - xi1(i) - z1(i) + w(i))
                ds(i) = -dx1(i)
-               dz1(i) = -z1(i) + (mu - z1(i)*dx1(i) 
+               dz1(i) = -z1(i) + (mu - z1(i)*dx1(i)
      &                  - dxdz1(i))/x1(i)
                dw(i) = -w(i) + (mu - w(i)*ds(i) - dsdw(i))/s(i)
             enddo
             do i = 1,n2
                dx2(i) = q2(i) * (dx2(i) - xi2(i) - r2(i))
-               dz2(i) = -z2(i) + (mu - z2(i)*dx2(i) - 
+               dz2(i) = -z2(i) + (mu - z2(i)*dx2(i) -
      &                  dxdz2(i))/x2(i)
             enddo
 c
@@ -479,7 +487,7 @@ subroutine slpfnc(n1,m,nnza1,a1,ja1,ia1,ao1,jao1,iao1,n2,nnza2,
          call daxpy(n1,deltad,dz1,1,z1,1)
          call daxpy(n2,deltad,dz2,1,z2,1)
          call daxpy(m,deltad,dy,1,y,1)
-         gap = ddot(n1,z1,1,x1,1) + ddot(n2,z2,1,x2,1) + 
+         gap = ddot(n1,z1,1,x1,1) + ddot(n2,z2,1,x2,1) +
      &         ddot(n1,w,1,s,1)
       goto 20
    30 continue