diff -Nru r-cran-gss-2.2-0/ChangeLog r-cran-gss-2.2-2/ChangeLog
--- r-cran-gss-2.2-0/ChangeLog	2020-05-03 18:34:45.000000000 +0000
+++ r-cran-gss-2.2-2/ChangeLog	2020-05-26 04:21:51.000000000 +0000
@@ -1,3 +1,17 @@
+Tue May 26 00:21:21 EDT 2020, Chong Gu  <chong@purdue.edu>
+
+	* DESCRIPTION: Version 2.2-2.
+
+	* R: Added sscomp and sscomp2.
+
+	* man: Added documentation for sscomp and sscomp2.
+
+Sat May 16 15:49:33 EDT 2020, Chong Gu  <chong@purdue.edu>
+
+	* DESCRIPTION: Version 2.2-1.
+
+	* R: Bug fixes, code cleanings in numerous places.
+
 Sun May  3 14:34:28 EDT 2020, Chong Gu  <chong@purdue.edu>
 
 	* DESCRIPTION: Version 2.2-0.
diff -Nru r-cran-gss-2.2-0/debian/changelog r-cran-gss-2.2-2/debian/changelog
--- r-cran-gss-2.2-0/debian/changelog	2020-05-05 01:36:33.000000000 +0000
+++ r-cran-gss-2.2-2/debian/changelog	2020-06-01 20:13:45.000000000 +0000
@@ -1,3 +1,17 @@
+r-cran-gss (2.2-2-1) unstable; urgency=medium
+
+  * New upstream release
+
+ -- Dirk Eddelbuettel <edd@debian.org>  Mon, 01 Jun 2020 15:13:45 -0500
+
+r-cran-gss (2.2-1-1) unstable; urgency=medium
+
+  * New upstream release
+
+  * debian/control: Set Build-Depends: to current R version
+
+ -- Dirk Eddelbuettel <edd@debian.org>  Thu, 21 May 2020 12:17:20 -0500
+
 r-cran-gss (2.2-0-1) unstable; urgency=medium
 
   * New upstream release
diff -Nru r-cran-gss-2.2-0/debian/control r-cran-gss-2.2-2/debian/control
--- r-cran-gss-2.2-0/debian/control	2020-03-05 19:45:13.000000000 +0000
+++ r-cran-gss-2.2-2/debian/control	2020-05-21 17:17:01.000000000 +0000
@@ -2,7 +2,7 @@
 Section: gnu-r
 Priority: optional
 Maintainer: Dirk Eddelbuettel <edd@debian.org>
-Build-Depends: debhelper (>= 7.0.0), r-base-dev (>= 3.6.3), dh-r
+Build-Depends: debhelper (>= 7.0.0), r-base-dev (>= 4.0.0), dh-r
 Standards-Version: 4.5.0
 Vcs-Browser: https://salsa.debian.org/edd/r-cran-gss
 Vcs-Git: https://salsa.debian.org/edd/r-cran-gss.git
diff -Nru r-cran-gss-2.2-0/DESCRIPTION r-cran-gss-2.2-2/DESCRIPTION
--- r-cran-gss-2.2-0/DESCRIPTION	2020-05-04 05:20:23.000000000 +0000
+++ r-cran-gss-2.2-2/DESCRIPTION	2020-05-26 05:20:03.000000000 +0000
@@ -1,6 +1,6 @@
 Package: gss
-Version: 2.2-0
-Date: 2020-05-03
+Version: 2.2-2
+Date: 2020-05-26
 Title: General Smoothing Splines
 Author: Chong Gu <chong@purdue.edu>
 Maintainer: Chong Gu <chong@purdue.edu>
@@ -8,7 +8,7 @@
 Description: A comprehensive package for structural multivariate
         function estimation using smoothing splines.
 License: GPL (>= 2)
-Packaged: 2020-05-03 22:48:29 UTC; chong
+Packaged: 2020-05-26 04:26:48 UTC; chong
 NeedsCompilation: yes
 Repository: CRAN
-Date/Publication: 2020-05-04 05:20:23 UTC
+Date/Publication: 2020-05-26 05:20:03 UTC
diff -Nru r-cran-gss-2.2-0/INDEX r-cran-gss-2.2-2/INDEX
--- r-cran-gss-2.2-0/INDEX	2020-04-30 16:56:11.000000000 +0000
+++ r-cran-gss-2.2-2/INDEX	2020-05-26 03:45:28.000000000 +0000
@@ -42,6 +42,8 @@
 cpssden                     Evaluating 1-D conditional cdf of ssden estimates
 cqssden                     Evaluating 1-D conditional quantiles of ssden estimates
 print.ssden                 Print function for ssden objects
+sscomp                      Estimating composition with one sample
+sscomp2                     Estimating composition with multiple samples
 
 sscden                      Estimating conditional density using smoothing splines
 d.sscden                    Evaluating pdf of sscden estimates
diff -Nru r-cran-gss-2.2-0/man/sscomp.Rd r-cran-gss-2.2-2/man/sscomp.Rd
--- r-cran-gss-2.2-0/man/sscomp.Rd	1970-01-01 00:00:00.000000000 +0000
+++ r-cran-gss-2.2-2/man/sscomp.Rd	2020-05-26 04:26:40.000000000 +0000
@@ -0,0 +1,33 @@
+\name{sscomp}
+\alias{sscomp}
+\alias{sscomp2}
+\title{Composition Estimation}
+\description{
+    Estimate composition using multinomial counts.
+}
+\usage{
+sscomp(x,wt=rep(1,length(x)),alpha=1.4)
+
+sscomp2(x,alpha=1.4)
+}
+\arguments{
+    \item{x}{Numerical vector or matrix of multinomial counts.}
+    \item{wt}{Numerical vector of integration weights.}
+    \item{alpha}{Parameter defining cross-validation score for smoothing
+      parameter selection.}
+}
+\details{
+    \code{sscomp} takes a vector \code{x} to estimate composition using
+    density estimation on a nominal discrete domain; zero counts must be
+    included in \code{x} to specify the domain.  \code{wt} mimicking the
+    shape of the unknown density could improve performance.
+
+    \code{sscomp2} takes a matrix \code{x}, collapses columns to
+    estimate a density using \code{sscomp}, then using that as \code{wt}
+    in further \code{sscomp} calls to estimate composition for each
+    column.
+}
+\references{
+    Gu, C. (2020), Composition estimation via shrinkage.
+    \emph{manuscript}.
+}
diff -Nru r-cran-gss-2.2-0/MD5 r-cran-gss-2.2-2/MD5
--- r-cran-gss-2.2-0/MD5	2020-05-04 05:20:23.000000000 +0000
+++ r-cran-gss-2.2-2/MD5	2020-05-26 05:20:03.000000000 +0000
@@ -1,15 +1,15 @@
-4beab7f87614ba21ae92f4879fe0e0df *ChangeLog
-5f847a7ff2964a8322acd73461e492bc *DESCRIPTION
-0f6b858d42fc1e6af7774e5d4a9c3300 *INDEX
-53174290c2c5d3bddfe81254c7abfa53 *NAMESPACE
+22fc0985d2cc9446dd1a3c39ab6526d7 *ChangeLog
+aa20d4dabee3fd8acde750d7fc5d3d43 *DESCRIPTION
+24ff2b8f15521a49b59f6d19146ae6df *INDEX
+20c494cb2d9aaa1a65f87cf5afd95df6 *NAMESPACE
 1d00015b10d8256aaede525722621ffc *R/cdsscden.R
 03873ecc380382307d04da914fb5a866 *R/cdssden.R
 a4cf561453300bc881519be811e6e36c *R/dsscden.R
 28afb3bc3469b09d1db335b5bfe4c88e *R/dsscopu.R
 c887032085de3b8aa5728103301014ab *R/dssden.R
-66c8a55c9254540eb7cae2930446ae77 *R/family.R
-53e1b86f0b4f1981ec2da0fc5e7fb1d0 *R/family.cv.R
-bed7aa9ccb46dce38cd135aefe71349d *R/family.proj.R
+559874c71d989cb315af06eca7797c56 *R/family.R
+8cf03c37c38512d84831d3ddad9c33b7 *R/family.cv.R
+e845e29cbedeefa26e421410d897ebf5 *R/family.proj.R
 630dbcbb2a0ced555d40e693b7b3f4d4 *R/family.surv.R
 d67212f47d1a4287b478740f0fb89992 *R/fitted.R
 7b557a66a85c9a87034d408e71e30a95 *R/gauss.quad.R
@@ -36,19 +36,20 @@
 500416a79d3297d43b5f35c78b640fda *R/project.ssanova.R
 023ed96bb4991e538c3f38a7a8bb3224 *R/project.ssanova9.R
 5a78c5b0ac10960f45efb305efa74b00 *R/project.sscden.R
-523f48d90ebe5792f08cf63e387de5f1 *R/project.sscden1.R
+efcfe6ca48d01c38d4db04c26e138805 *R/project.sscden1.R
 50ed5e000696a4ec366834a09673e448 *R/project.sscox.R
 89a6175fe772e65c76ab3eb644acd316 *R/project.ssden.R
-c95e17013f7411b89d0f8e14b6d2a3b1 *R/project.ssden1.R
+8133bee08081d95af25fd6d2a1f905de *R/project.ssden1.R
 fdc26fdef9a730f792d3c3bd146c0940 *R/project.sshzd.R
-5d77baee61428e3359defc8a7599c68b *R/project.sshzd1.R
+794230b5f6c1e47dfca0fef94ebbf311 *R/project.sshzd1.R
 bf5f875bbef9fb78e75bafa5aa4be563 *R/project.ssllrm.R
 d44fab9e998870ebdb9a3c393d0ef756 *R/smolyak.R
-c1e008d561780b72c7b07a39e7562b3f *R/ssanova.R
+912a7cb4edd645076b72eaae69d220e6 *R/ssanova.R
 dd4aa54b10ec5e2a8b0cdff18780d562 *R/ssanova0.R
 9be84d7ff1b4e8a5922bd678ace8b40f *R/ssanova9.R
 edb14e97f6de2cb36d610cc02725ba3a *R/sscden.R
 40ddb4f3d1511288c96182020673fa15 *R/sscden1.R
+1cbca9d45e7d735cd750a1f29c8b786a *R/sscomp.R
 58bfc87ec8debc40bdf3ad95f4df085e *R/sscopu.R
 111d09029fc69db43413326a418fb286 *R/sscopu2.R
 b9fad02c0244197802cc8fd0fcdde239 *R/sscox.R
@@ -136,6 +137,7 @@
 e2f739c757a04a968ec837b15ed91985 *man/ssanova0.Rd
 e54aa57ed8c94851cb790cd9de84b871 *man/ssanova9.Rd
 a56a56fb459ddc58c95c27e8b9e37ef8 *man/sscden.Rd
+a89272e0fa273b6524da0e23574f96bb *man/sscomp.Rd
 142a1be80d200caf0cbd80218c2db7bc *man/sscopu.Rd
 4df29e6fff4a4655285f24904cb212c3 *man/sscox.Rd
 97f4d0c57cfa8fac287275a9588eca88 *man/ssden.Rd
@@ -153,7 +155,6 @@
 0423a47eb8c93a22c1ec15136a2419f5 *src/cdennewton.f
 ac411cbbffa81ff956cb59e79f67403c *src/cdennewton10.f
 08b22168d8e387134de819afae34c70d *src/copu2newton.f
-1d0848b9a9c7ec70927d1385723b00fd *src/dchdc0.f
 7660189020bc629e33531dfccc58721f *src/dcoef.f
 b95abc103e26f67dc5c543096d0a9967 *src/dcore.f
 b05ac68613d1536384c72933d4d089be *src/dcrdr.f
@@ -178,7 +179,7 @@
 53df31681eb9b0911ec8bb3b5e9cc0c5 *src/hzdaux.f
 44cf970d0d0271b391ccdfc5475d3561 *src/hzdnewton.f
 fa5c780f4bb48e1aa23fa7cd9be53e57 *src/hzdnewton10.f
-16be7665de52ce87120138bcda100182 *src/init.c
+f3314da76ddac6f573346ba745289e56 *src/init.c
 946e7257e372eac1629e256df6afde4a *src/llrmnewton.f
 88109ec6b4f9ff4bd8b525f0a866c332 *src/ratfor/cdennewton.r
 025d7e6ce3a32679a292ebb44bb3e4ec *src/ratfor/cdennewton10.r
diff -Nru r-cran-gss-2.2-0/NAMESPACE r-cran-gss-2.2-2/NAMESPACE
--- r-cran-gss-2.2-0/NAMESPACE	2020-05-02 19:09:23.000000000 +0000
+++ r-cran-gss-2.2-2/NAMESPACE	2020-05-26 03:42:50.000000000 +0000
@@ -10,7 +10,7 @@
        predict1, predict9, project, psscden, pssden, qsscden, qssden,
        smolyak.quad, smolyak.size, ssanova, ssanova0, ssanova9,
        sscden, sscden1, sscox, ssden, ssden1, sshzd, sshzd1, ssllrm,
-       survexp.sshzd)
+       survexp.sshzd,sscomp,sscomp2)
 export(sscopu, sscopu2, dsscopu, cdsscopu, cpsscopu, cqsscopu,
        sshzd2d, sshzd2d1, hzdrate.sshzd2d, survexp.sshzd2d)
 # export internal functions used in examples
diff -Nru r-cran-gss-2.2-0/R/family.cv.R r-cran-gss-2.2-2/R/family.cv.R
--- r-cran-gss-2.2-0/R/family.cv.R	2019-06-26 05:03:48.000000000 +0000
+++ r-cran-gss-2.2-2/R/family.cv.R	2020-05-16 13:32:59.000000000 +0000
@@ -46,19 +46,15 @@
     v <- t(sr)%*%(wt*w*sr)/sum(wt*w)-outer(mu,mu)
     v[(nnull+1):nn,(nnull+1):nn] <- v[(nnull+1):nn,(nnull+1):nn]+q/sum(wt*y)
     ## Cholesky decomposition of H
-    z <- .Fortran("dchdc0",
-                  v=as.double(v), as.integer(nn), as.integer(nn),
-                  double(nn), jpvt=as.integer(rep(0,nn)),
-                  as.integer(1), rkv=integer(1),
-                  PACKAGE="gss")[c("v","jpvt","rkv")]
-    v <- matrix(z$v,nn,nn)
-    rkv <- z$rkv
+    suppressWarnings(z <- chol(v,pivot=TRUE))
+    v <- z
+    rkv <- attr(z,"rank")
     while (v[rkv,rkv]<v[1,1]*sqrt(.Machine$double.eps)) rkv <- rkv-1
     if (rkv<nn) v[(rkv+1):nn,(rkv+1):nn] <- diag(v[1,1],nn-rkv)
     ## trace
     mu <- apply(wt*y*sr,2,sum)/sum(wt*y)
     sr <- sqrt(wt*y)*t(t(sr)-mu)
-    sr <- backsolve(v,t(sr[,z$jpvt]),transpose=TRUE)
+    sr <- backsolve(v,t(sr[,attr(z,"pivot")]),transpose=TRUE)
     aux1 <- sum(sr^2)
     aux2 <- 1/sum(wt*y)/(sum(wt*y)-1)
     list(score=lkhd+abs(alpha)*aux1*aux2,varht=1,w=as.vector(wt*w))
diff -Nru r-cran-gss-2.2-0/R/family.proj.R r-cran-gss-2.2-2/R/family.proj.R
--- r-cran-gss-2.2-0/R/family.proj.R	2019-07-11 16:52:42.000000000 +0000
+++ r-cran-gss-2.2-2/R/family.proj.R	2020-05-16 14:15:08.000000000 +0000
@@ -224,7 +224,7 @@
     for (i in 2:(nnu+1)) y <- cbind(y,y0[,i]-y0[,i-1])
     y <- cbind(y,1-y0[,nnu+1])
     lkhd <- function(log.nu) {
-        nu <- exp(nu)
+        nu <- exp(log.nu)
         G <- c(0,cumsum(nu))
         P <- exp(outer(eta,G,"+"))
         lkhd <- 0
@@ -256,7 +256,8 @@
     P <- P/(1+P)
     for (i in 1:length(eta)) {
         tmp <- diff(c(0,P[i,],1))
-        kl <- kl+wt[i]*sum(y[i,]*log(y[i,]/tmp))
+        if (min(tmp)<=0) kl <- Inf
+        else kl <- kl+wt[i]*sum(y[i,]*log(y[i,]/tmp))
     }
     kl <- kl/sum(wt)
     wt <- w*wt
diff -Nru r-cran-gss-2.2-0/R/family.R r-cran-gss-2.2-2/R/family.R
--- r-cran-gss-2.2-0/R/family.R	2020-05-03 22:45:37.000000000 +0000
+++ r-cran-gss-2.2-2/R/family.R	2020-05-16 19:40:01.000000000 +0000
@@ -301,6 +301,9 @@
     dd <- log(nu)
     repeat {
         ## gradient and hessian
+        nu <- exp(dd)
+        G <- c(0,cumsum(nu))
+        P <- exp(outer(eta,G,"+"))
         grad <- hess.wk <- NULL
         for (i in 1:nnu) {
             g.wk <- h.wk <- 0
@@ -335,29 +338,36 @@
         else hess <- abs(hess)
         ## update nu
         mumax <- max(abs(grad))
+        dd.diff <- solve(hess,grad)
         repeat {
-            ddnew <- dd-solve(hess,grad)
-            nu <- exp(ddnew)
-            G <- c(0,cumsum(nu))
-            P <- exp(outer(eta,G,"+"))
-            lkhdnew <- 0
-            for (i in 1:(nnu+1))
-                lkhdnew <- lkhdnew+sum(wt*(y[,i]+y[,i+1])*log(1+P[,i]))/sum(wt)
-            for (i in 1:nnu) lkhdnew <- lkhdnew-sum(wt*y[,i+1])/sum(wt)*log(exp(nu[i])-1)
-            if (nnu>1) {
-                for (i in 1:(nnu-1)) {
-                    tmp <- 0
-                    for (j in (i+1):nnu) tmp <- tmp+sum(wt*y[,j+1])/sum(wt)
-                    lkhdnew <- lkhdnew-tmp*nu[i]
+            lkhd.line <- function(x) {
+                ddnew <- dd-c(x)*dd.diff
+                nu <- exp(ddnew)
+                G <- c(0,cumsum(nu))
+                P <- exp(outer(eta,G,"+"))
+                lkhd <- 0
+                for (i in 1:(nnu+1))
+                    lkhd <- lkhd+sum(wt*(y[,i]+y[,i+1])*log(1+P[,i]))/sum(wt)
+                for (i in 1:nnu) lkhd <- lkhd-sum(wt*y[,i+1])/sum(wt)*log(exp(nu[i])-1)
+                if (nnu>1) {
+                    for (i in 1:(nnu-1)) {
+                        tmp <- 0
+                        for (j in (i+1):nnu) tmp <- tmp+sum(wt*y[,j+1])/sum(wt)
+                        lkhd <- lkhd-tmp*nu[i]
+                    }
                 }
+                lkhd
             }
-            if (!is.finite(lkhdnew)) {
-                grad <- grad/2
+            if (!is.finite(lkhdnew <- lkhd.line(1))) {
+                dd.diff <- dd.diff/2
                 next
             }
-            if (lkhdnew-lkhd<(1+abs(lkhd)*10*.Machine$double.eps)) break
-            grad <- grad/2
-            if (max(grad)/mumax<10*.Machine$double.eps) break
+            ddnew <- dd-dd.diff
+            if (lkhdnew-lkhd<(1+abs(lkhd))*10*.Machine$double.eps) break
+            z <- nlm0(lkhd.line,c(0,1))
+            ddnew <- dd-z$est*dd.diff
+            lkhdnew <- z$min
+            break
         }
         disc <- abs(lkhdnew-lkhd)/(1+abs(lkhd))
         disc <- max(disc,max(abs(dd-ddnew)/(1+abs(dd))))
diff -Nru r-cran-gss-2.2-0/R/project.sscden1.R r-cran-gss-2.2-2/R/project.sscden1.R
--- r-cran-gss-2.2-0/R/project.sscden1.R	2014-07-03 23:05:04.000000000 +0000
+++ r-cran-gss-2.2-2/R/project.sscden1.R	2020-05-16 13:21:37.000000000 +0000
@@ -162,7 +162,7 @@
             mu <- s.eta[id.s0]
         }
         nn <- length(mu)
-        z <- chol(v,pivot=TRUE)
+        suppressWarnings(z <- chol(v,pivot=TRUE))
         v <- z
         rkv <- attr(z,"rank")
         m.eps <- .Machine$double.eps
diff -Nru r-cran-gss-2.2-0/R/project.ssden1.R r-cran-gss-2.2-2/R/project.ssden1.R
--- r-cran-gss-2.2-0/R/project.ssden1.R	2014-07-03 23:05:04.000000000 +0000
+++ r-cran-gss-2.2-2/R/project.ssden1.R	2020-05-16 13:22:11.000000000 +0000
@@ -70,7 +70,7 @@
         v <- cbind(rbind(ss.wk,t(sr.wk)),rbind(sr.wk,rr.wk))
         mu <- c(s.eta[id.s],r.eta.wk)
         nn <- length(mu)
-        z <- chol(v,pivot=TRUE)
+        suppressWarnings(z <- chol(v,pivot=TRUE))
         v <- z
         rkv <- attr(z,"rank")
         m.eps <- .Machine$double.eps
diff -Nru r-cran-gss-2.2-0/R/project.sshzd1.R r-cran-gss-2.2-2/R/project.sshzd1.R
--- r-cran-gss-2.2-0/R/project.sshzd1.R	2014-07-03 23:05:04.000000000 +0000
+++ r-cran-gss-2.2-2/R/project.sshzd1.R	2020-05-16 13:22:51.000000000 +0000
@@ -137,7 +137,7 @@
         v <- cbind(rbind(ss.wk,t(sr.wk)),rbind(sr.wk,rr.wk))
         mu <- c(s.eta[id.s],r.eta.wk)
         nn <- length(mu)
-        z <- chol(v,pivot=TRUE)
+        suppressWarnings(z <- chol(v,pivot=TRUE))
         v <- z
         rkv <- attr(z,"rank")
         m.eps <- .Machine$double.eps
diff -Nru r-cran-gss-2.2-0/R/ssanova.R r-cran-gss-2.2-2/R/ssanova.R
--- r-cran-gss-2.2-0/R/ssanova.R	2014-07-03 23:05:04.000000000 +0000
+++ r-cran-gss-2.2-2/R/ssanova.R	2020-05-16 13:26:16.000000000 +0000
@@ -145,7 +145,7 @@
                 10^(2*ran.scal)*random$sigma$fun(lambda[-1],random$sigma$env)
         }
         if (qr.trace) {
-            qq.wk <- chol(q.wk,pivot=TRUE)
+            suppressWarnings(qq.wk <- chol(q.wk,pivot=TRUE))
             sr <- cbind(s,10^theta*r[,attr(qq.wk,"pivot")])
             sr <- rbind(sr,cbind(matrix(0,nxiz,nnull),qq.wk))
             sr <- qr(sr,tol=0)
@@ -324,7 +324,7 @@
             r.wk0 <- wt*r.wk0
         }
         if (qr.trace) {
-            qq.wk <- chol(q.wk,pivot=TRUE)
+            suppressWarnings(qq.wk <- chol(q.wk,pivot=TRUE))
             sr <- cbind(s.wk,r.wk0[,attr(qq.wk,"pivot")])
             sr <- rbind(sr,cbind(matrix(0,nxiz,nnull),qq.wk))
             sr <- qr(sr,tol=0)
@@ -505,7 +505,7 @@
     wk1[1:nnull,1:nnull] <- sms
     wk1[1:nnull,nnull+(1:rkq)] <- -t(t(dr)/val)
     wk1[nnull+(1:rkq),nnull+(1:rkq)] <- wk2
-    z <- chol(wk1,pivot=TRUE)
+    suppressWarnings(z <- chol(wk1,pivot=TRUE))
     wk1 <- z
     rkw <- attr(z,"rank")
     while (wk1[rkw,rkw]<wk1[1,1]*sqrt(.Machine$double.eps)) rkw <- rkw-1
diff -Nru r-cran-gss-2.2-0/R/sscomp.R r-cran-gss-2.2-2/R/sscomp.R
--- r-cran-gss-2.2-0/R/sscomp.R	1970-01-01 00:00:00.000000000 +0000
+++ r-cran-gss-2.2-2/R/sscomp.R	2020-05-26 03:28:36.000000000 +0000
@@ -0,0 +1,48 @@
+## Composition estimation with one sample
+sscomp <- function(x,wt=rep(1,length(x)),alpha=1.4)
+{
+    ## Check inputs
+    if ((nlvl <- length(x))<3)
+        stop("gss error in sscomp: length of x should be 3 or more")
+    if (length(x)!=length(wt))
+        stop("gss error in sscomp: x and wt mismatch in lengths")
+    ## Generate terms
+    cnt <- x
+    x <- as.factor(1:nlvl)
+    mf <- model.frame(~x)
+    term <- mkterm(mf,NULL)
+    rk <- term$x$rk
+    x.basis <- x[-nlvl]
+    r <- rk$fun(x.basis,x,nu=1,env=rk$env,out=TRUE)
+    q <- r[,-nlvl]
+    qd.wt <- as.vector(wt)
+    ## Fit the model
+    nt <- b.wt <- 1
+    t.wt <- matrix(1,nlvl,1)
+    bias0 <- list(nt=nt,wt=b.wt,qd.wt=t.wt)
+    z <- sspdsty(NULL,r,q,cnt,NULL,r,qd.wt,1e-7,30,alpha,bias0)
+    ## return fitted probabilities
+    fit <- exp(t(r)%*%z$c)*qd.wt
+    rownames(fit) <- rownames(x)
+    fit/sum(fit)
+}
+## Composition estimation with a matrix input
+sscomp2 <- function(x,alpha=1.4)
+{
+    if (!is.matrix(x))
+        stop("gss error in sscomp2: x should be a matrix")
+    if (min(x)<0)
+        stop("gss error in sscomp2: x should have non-negative entries")
+    if (any(apply(x,2,sum)==0))
+        stop("gss error in sscomp2: column totals of x must be positive")
+    nlvl <- dim(x)[1]
+    yy <- apply(x,1,sum)
+    p0 <- sscomp(yy)
+    fit <- NULL
+    for (i in 1:dim(x)[2]) {
+        fit <- cbind(fit,sscomp(x[,i],p0,alpha))
+    }
+    rownames(fit) <- rownames(x)
+    colnames(fit) <- colnames(x)
+    fit
+}
diff -Nru r-cran-gss-2.2-0/src/dchdc0.f r-cran-gss-2.2-2/src/dchdc0.f
--- r-cran-gss-2.2-0/src/dchdc0.f	2019-06-26 03:04:40.000000000 +0000
+++ r-cran-gss-2.2-2/src/dchdc0.f	1970-01-01 00:00:00.000000000 +0000
@@ -1,234 +0,0 @@
-      subroutine dchdc0(a,lda,p,work,jpvt,job,info)
-      integer lda,p,jpvt(1),job,info
-      double precision a(lda,1),work(1)
-c
-c     dchdc computes the cholesky decomposition of a positive definite
-c     matrix.  a pivoting option allows the user to estimate the
-c     condition of a positive definite matrix or determine the rank
-c     of a positive semidefinite matrix.
-c
-c     on entry
-c
-c         a      double precision(lda,p).
-c                a contains the matrix whose decomposition is to
-c                be computed.  onlt the upper half of a need be stored.
-c                the lower part of the array a is not referenced.
-c
-c         lda    integer.
-c                lda is the leading dimension of the array a.
-c
-c         p      integer.
-c                p is the order of the matrix.
-c
-c         work   double precision.
-c                work is a work array.
-c
-c         jpvt   integer(p).
-c                jpvt contains integers that control the selection
-c                of the pivot elements, if pivoting has been requested.
-c                each diagonal element a(k,k)
-c                is placed in one of three classes according to the
-c                value of jpvt(k).
-c
-c                   if jpvt(k) .gt. 0, then x(k) is an initial
-c                                      element.
-c
-c                   if jpvt(k) .eq. 0, then x(k) is a free element.
-c
-c                   if jpvt(k) .lt. 0, then x(k) is a final element.
-c
-c                before the decomposition is computed, initial elements
-c                are moved by symmetric row and column interchanges to
-c                the beginning of the array a and final
-c                elements to the end.  both initial and final elements
-c                are frozen in place during the computation and only
-c                free elements are moved.  at the k-th stage of the
-c                reduction, if a(k,k) is occupied by a free element
-c                it is interchanged with the largest free element
-c                a(l,l) with l .ge. k.  jpvt is not referenced if
-c                job .eq. 0.
-c
-c        job     integer.
-c                job is an integer that initiates column pivoting.
-c                if job .eq. 0, no pivoting is done.
-c                if job .ne. 0, pivoting is done.
-c
-c     on return
-c
-c         a      a contains in its upper half the cholesky factor
-c                of the matrix a as it has been permuted by pivoting.
-c
-c         jpvt   jpvt(j) contains the index of the diagonal element
-c                of a that was moved into the j-th position,
-c                provided pivoting was requested.
-c
-c         info   contains the index of the last positive diagonal
-c                element of the cholesky factor.
-c
-c     for positive definite matrices info = p is the normal return.
-c     for pivoting with positive semidefinite matrices info will
-c     in general be less than p.  however, info may be greater than
-c     the rank of a, since rounding error can cause an otherwise zero
-c     element to be positive. indefinite systems will always cause
-c     info to be less than p.
-c
-c     linpack. this version dated 08/14/78 .
-c     j.j. dongarra and g.w. stewart, argonne national laboratory and
-c     university of maryland.
-c
-c
-c     blas daxpy,dswap
-c     fortran dsqrt
-c
-c     internal variables
-c
-      integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl
-      double precision temp
-      double precision maxdia
-      logical swapk,negk
-c
-      pl = 1
-      pu = 0
-      info = p
-      if (job .eq. 0) go to 160
-c
-c        pivoting has been requested. rearrange the
-c        the elements according to jpvt.
-c
-         do 70 k = 1, p
-            swapk = jpvt(k) .gt. 0
-            negk = jpvt(k) .lt. 0
-            jpvt(k) = k
-            if (negk) jpvt(k) = -jpvt(k)
-            if (.not.swapk) go to 60
-               if (k .eq. pl) go to 50
-                  call dswap(pl-1,a(1,k),1,a(1,pl),1)
-                  temp = a(k,k)
-                  a(k,k) = a(pl,pl)
-                  a(pl,pl) = temp
-                  plp1 = pl + 1
-                  if (p .lt. plp1) go to 40
-                  do 30 j = plp1, p
-                     if (j .ge. k) go to 10
-                        temp = a(pl,j)
-                        a(pl,j) = a(j,k)
-                        a(j,k) = temp
-                     go to 20
-   10                continue
-                     if (j .eq. k) go to 20
-                        temp = a(k,j)
-                        a(k,j) = a(pl,j)
-                        a(pl,j) = temp
-   20                continue
-   30             continue
-   40             continue
-                  jpvt(k) = jpvt(pl)
-                  jpvt(pl) = k
-   50          continue
-               pl = pl + 1
-   60       continue
-   70    continue
-         pu = p
-         if (p .lt. pl) go to 150
-         do 140 kb = pl, p
-            k = p - kb + pl
-            if (jpvt(k) .ge. 0) go to 130
-               jpvt(k) = -jpvt(k)
-               if (pu .eq. k) go to 120
-                  call dswap(k-1,a(1,k),1,a(1,pu),1)
-                  temp = a(k,k)
-                  a(k,k) = a(pu,pu)
-                  a(pu,pu) = temp
-                  kp1 = k + 1
-                  if (p .lt. kp1) go to 110
-                  do 100 j = kp1, p
-                     if (j .ge. pu) go to 80
-                        temp = a(k,j)
-                        a(k,j) = a(j,pu)
-                        a(j,pu) = temp
-                     go to 90
-   80                continue
-                     if (j .eq. pu) go to 90
-                        temp = a(k,j)
-                        a(k,j) = a(pu,j)
-                        a(pu,j) = temp
-   90                continue
-  100             continue
-  110             continue
-                  jt = jpvt(k)
-                  jpvt(k) = jpvt(pu)
-                  jpvt(pu) = jt
-  120          continue
-               pu = pu - 1
-  130       continue
-  140    continue
-  150    continue
-  160 continue
-      do 270 k = 1, p
-c
-c        reduction loop.
-c
-         maxdia = a(k,k)
-         kp1 = k + 1
-         maxl = k
-c
-c        determine the pivot element.
-c
-         if (k .lt. pl .or. k .ge. pu) go to 190
-            do 180 l = kp1, pu
-               if (a(l,l) .le. maxdia) go to 170
-                  maxdia = a(l,l)
-                  maxl = l
-  170          continue
-  180       continue
-  190    continue
-c
-c        quit if the pivot element is not positive.
-c
-         if (maxdia .gt. 0.0d0) go to 200
-            info = k - 1
-c     ......exit
-            go to 280
-  200    continue
-         if (k .eq. maxl) go to 210
-c
-c           start the pivoting and update jpvt.
-c
-            km1 = k - 1
-            call dswap(km1,a(1,k),1,a(1,maxl),1)
-            a(maxl,maxl) = a(k,k)
-            a(k,k) = maxdia
-            jp = jpvt(maxl)
-            jpvt(maxl) = jpvt(k)
-            jpvt(k) = jp
-  210    continue
-c
-c        reduction step. pivoting is contained across the rows.
-c
-         work(k) = dsqrt(a(k,k))
-         a(k,k) = work(k)
-         if (p .lt. kp1) go to 260
-         do 250 j = kp1, p
-            if (k .eq. maxl) go to 240
-               if (j .ge. maxl) go to 220
-                  temp = a(k,j)
-                  a(k,j) = a(j,maxl)
-                  a(j,maxl) = temp
-               go to 230
-  220          continue
-               if (j .eq. maxl) go to 230
-                  temp = a(k,j)
-                  a(k,j) = a(maxl,j)
-                  a(maxl,j) = temp
-  230          continue
-  240       continue
-            a(k,j) = a(k,j)/work(k)
-            work(j) = a(k,j)
-            temp = -a(k,j)
-            call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1)
-  250    continue
-  260    continue
-  270 continue
-  280 continue
-      return
-      end
diff -Nru r-cran-gss-2.2-0/src/init.c r-cran-gss-2.2-2/src/init.c
--- r-cran-gss-2.2-0/src/init.c	2020-04-25 17:16:10.000000000 +0000
+++ r-cran-gss-2.2-2/src/init.c	2020-05-16 12:22:22.000000000 +0000
@@ -36,7 +36,6 @@
 extern void F77_NAME(reg)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
 extern void F77_NAME(regaux)(void *, void *, void *, void *, void *, void *, void *, void *, void *);
 extern void F77_NAME(dmcdc)(void *, void *, void *, void *, void *, void *);
-extern void F77_NAME(dchdc0)(void *, void *, void *, void *, void *, void *, void *);
 
 static const R_CMethodDef CEntries[] = {
     {"quad_smolyak", (DL_FUNC) &quad_smolyak, 4},
@@ -70,7 +69,6 @@
     {"reg",          (DL_FUNC) &F77_NAME(reg),          18},
     {"regaux",       (DL_FUNC) &F77_NAME(regaux),        9},
     {"dmcdc",        (DL_FUNC) &F77_NAME(dmcdc),         6},
-    {"dchdc0",        (DL_FUNC) &F77_NAME(dchdc0),       7},
     {NULL, NULL, 0}
 };