diff --git a/DESCRIPTION b/DESCRIPTION index 2bcddacd..ef2077d5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: Matrix -Version: 1.2-16 -Date: 2019-03-04 +Version: 1.2-17 +Date: 2019-03-11 Priority: recommended Title: Sparse and Dense Matrix Classes and Methods Contact: Doug and Martin @@ -34,6 +34,8 @@ BuildResaveData: no License: GPL (>= 2) | file LICENCE URL: http://Matrix.R-forge.R-project.org/ BugReports: https://r-forge.r-project.org/tracker/?group_id=61 +NeedsCompilation: yes +Packaged: 2019-03-20 21:37:04 UTC; maechler Author: Douglas Bates [aut], Martin Maechler [aut, cre] (), Timothy A. Davis [ctb] (SuiteSparse and 'cs' C libraries, notably @@ -45,9 +47,4 @@ Author: Douglas Bates [aut], Regents of the University of California), R Core Team [ctb] (base R matrix implementation) Repository: CRAN -Repository/R-Forge/Project: matrix -Repository/R-Forge/Revision: 3291 -Repository/R-Forge/DateTimeStamp: 2019-03-04 11:00:53 -Date/Publication: 2019-03-08 15:33:45 UTC -NeedsCompilation: yes -Packaged: 2019-03-04 11:32:46 UTC; rforge +Date/Publication: 2019-03-22 22:56:52 UTC diff --git a/MD5 b/MD5 index 1d833398..33e8d284 100644 --- a/MD5 +++ b/MD5 @@ -1,10 +1,10 @@ 675a3a48a357cb2ee05ed285387b3e85 *ChangeLog -1a6f6fa0b7de869ebb6b603fb2cc077d *DESCRIPTION +5e59faf86756d7dcf42826b8d318c0ec *DESCRIPTION 49c7113b5e0e4fbf269132f6a5ff5b06 *LICENCE c64f1d9e5087f242f8a60c4f135fbf10 *NAMESPACE 7e494fcaf80175fe73944c4d8732eb90 *R/AllClass.R b8ddeb8ff877f40a12519c3e965debb2 *R/AllGeneric.R -6b4d466dfdba769a0c6aeeb57b7c8e3f *R/Auxiliaries.R +f46239d6b3b7f1d7af26980777bacbe6 *R/Auxiliaries.R e184b3cf64c092b6ed0a0d22eb332a44 *R/CHMfactor.R 3d80552f2d6bbaef906bd95dc4b734f6 *R/Csparse.R d7497671f04e08d0d703b650fd3eacff *R/HBMM.R @@ -12,7 +12,7 @@ d7497671f04e08d0d703b650fd3eacff *R/HBMM.R 4d794b4715d6efa432fd7cf86f69da7a *R/KhatriRao.R 7cd5d5938844fa73c1a2b50fae7fdd19 *R/LU.R 639a65f83a4dccc663aad5299708071b *R/Math.R -9eb05df834157e3129ed758df50d879b *R/Matrix.R +8090100a5f06f969fc92a686e555800a *R/Matrix.R e19ed36934229e22f41552d4d40c5d5a *R/MatrixFactorization.R 6e6e090d031ce8f02c8b4a6cd5b6ab08 *R/Ops.R a738e970d3c972b9ce5bc9492024a724 *R/Rsparse.R @@ -76,9 +76,9 @@ b3d2c13789b1cc7ad4b38d1ee14a07bd *R/sparseMatrix.R 7fcad2d91e94cfac03a4276a6acb211c *R/sparseVector.R 97e508e9e124aac0fb3dc3d046de347f *R/symmetricMatrix.R a3265f0761ef916f610e9849c6677600 *R/triangularMatrix.R -97734dbeedf076622f5a20e8cd6c0937 *R/zzz.R -8531259eb72b7318c9028de5c7584cb5 *TODO -53cea45c7db3841b03c8aecb778ad256 *build/partial.rdb +7bb8f3b96f0cd3fc41bf5e4d75d085d9 *R/zzz.R +f3b1b137d928891510363a3a18f9cfa6 *TODO +77fc9f3d733ea9855b5b5026431a4eea *build/partial.rdb f6a48395281382c8210c31d00fb9d7a7 *build/vignette.rds f2ad5375e270deeb7b041272dd095032 *cleanup 53de12cc7d1449cb2fa136f09333286b *data/CAex.R @@ -86,20 +86,20 @@ f2ad5375e270deeb7b041272dd095032 *cleanup 3ed2e24d0f9b1af73f0b432f99644d64 *data/USCounties.R 8734f0b040c6292983d273d4251d250a *inst/Copyrights ac3697580885d17a062ce39df01338ad *inst/Doxyfile -e53099d1529073886025d47d1b81dd0e *inst/NEWS.Rd +76ebc2b80bf1698523d3c104a0f096ff *inst/NEWS.Rd ea747a35318c564002a7144ebd510e28 *inst/doc/Announce.txt 0eeea40a3fb6da7157d2862bba10713a *inst/doc/Comparisons.R 09e13507b9832a30aa92e8aff0dc96a0 *inst/doc/Comparisons.Rnw -7a81c839f873c82768bb567a9fdb71d2 *inst/doc/Comparisons.pdf +42d1e2ebceb4386a7e81b12d60ee5b6b *inst/doc/Comparisons.pdf f5541439d0f08315e57f86d026a1c40d *inst/doc/Design-issues.R 5a9e99ce004c325e01fb2773e58052f5 *inst/doc/Design-issues.Rnw -c209eeefbc76b7286a53031634725995 *inst/doc/Design-issues.pdf +e6e6f6671bf9b0075f6b688a6deec860 *inst/doc/Design-issues.pdf 6f66bd95aab28d7bf06a019e7fbc14b5 *inst/doc/Intro2Matrix.R a6d6903590a8f5ac036fdeb74b48bb89 *inst/doc/Intro2Matrix.Rnw -a4a3a928d25bee227577e54e8985c62c *inst/doc/Intro2Matrix.pdf +ea6a658bb1c0874563f623eabd3a0e38 *inst/doc/Intro2Matrix.pdf 188b839523ce6ec32e1f00ec4e2bc18e *inst/doc/Introduction.R c39a26dfe7ccaafd044e88468155b153 *inst/doc/Introduction.Rnw -bd180d3d5a75876b6201a8458d119316 *inst/doc/Introduction.pdf +7a0268bde26480ecdbbb0ffa47dbcbeb *inst/doc/Introduction.pdf 7d60c8defc517d361513bac03d04f13b *inst/doc/SuiteSparse/AMD.txt e4f8cd28fc8be8ab14106630610b9c5f *inst/doc/SuiteSparse/CHOLMOD.txt 0a80fb0db9cefb02e290ffd2f3431e6c *inst/doc/SuiteSparse/COLAMD.txt @@ -108,7 +108,7 @@ d75882d4bb768ba0ea352291652daaee *inst/doc/SuiteSparse/SPQR.txt 6d217288f5da4fe419afaa34988bf42d *inst/doc/SuiteSparse/UserGuides.txt 2ab394c4daa5c94c2498ba3fa0578462 *inst/doc/sparseModels.R 491d961ca64f8da569c817b82daa9c58 *inst/doc/sparseModels.Rnw -62882f881b274478bd63ea10085b763f *inst/doc/sparseModels.pdf +2e592e2bec481ff8fe86d04499bea367 *inst/doc/sparseModels.pdf dcd11f6947f910f743254824e930b2c7 *inst/external/CAex_slots.rda be886d6bb832210bb654b9ad064fe0ff *inst/external/KNex_slots.rda 90f019ec81e67d7f3542f7ca39bf3f2d *inst/external/USCounties_slots.rda @@ -220,7 +220,7 @@ f5127ffbf702e88219a61c8385c53de2 *man/lu.Rd 4a897749bfc1ee9afd46bb58f8fa8c05 *man/matrix-products.Rd 4a449d40a733322c55cf7146d8ddb43e *man/nMatrix-class.Rd 3e94bd26a125e628e83c7b6d4d8a045c *man/ndenseMatrix-class.Rd -46dc1c9089fefb669bf1b8137edd2d12 *man/nearPD.Rd +115b74fe54e0b5ada18182d9bde630b7 *man/nearPD.Rd 07c563eefd115c5f74e5ba28231ffb46 *man/ngeMatrix-class.Rd ffb57c3063559175113661c71b6f130d *man/nnzero.Rd 001c81a2b8215c310e817aa60ea2095b *man/norm.Rd @@ -383,31 +383,31 @@ e14e53d0b36e1d0647348a326e8b3b54 *src/COLAMD/Include/colamd.h b56b67fa4e3e8bc168f5ad98d27f207a *src/COLAMD/Source/Makefile 561ff6512783cc2af5b343ff3630267d *src/COLAMD/Source/colamd.c 8f9d70daa9adafe39e555573d974d17e *src/COLAMD/Source/colamd_global.c -e52babc8b8e37bd088058bc96f39d73e *src/Csparse.c +60251cd94325ee15b1502f9ea8bf1184 *src/Csparse.c 684b703246b23a66b77fce3dc0e5d73a *src/Csparse.h e879c68d033040e33a588710797323b5 *src/Makevars 739a9b145f3e90adf5b58f743e4d3e77 *src/Mutils.c -5ea701c24ba16afb99f5164e7f1f4907 *src/Mutils.h +6b01b5ee807043e6cf42191f68306af9 *src/Mutils.h 34a94c6db1f395da28ca4e94dd9cf9ca *src/SuiteSparse_config/Makefile e9b2be5ad0056f588187cf504cf60764 *src/SuiteSparse_config/SuiteSparse_config.c ff586155b86948cc5b4eec822b8670e2 *src/SuiteSparse_config/SuiteSparse_config.h d41d8cd98f00b204e9800998ecf8427e *src/SuiteSparse_config/SuiteSparse_config.mk -b08d523906ab962f8e0729cb4d3d2f15 *src/Syms.h +1cf7c17bb7f79d7d22577cf4770cfc37 *src/Syms.h 5dfee1d8667f9553469dea08cec731e6 *src/TMatrix_as.c 0ef00bab2b7ffb79f64b25311decd236 *src/TMatrix_as.h 6551aecedc2b54c2c1f4c8a0973861f2 *src/Tsparse.c 20c550b57138bf22290ea759dfcb1c3e *src/Tsparse.h 3714d37c755fbf0b5b2a88362b651390 *src/abIndex.c ea6bdd1b33f38c8c2913eeef6eb21be3 *src/abIndex.h -f6f9e84779e01f8dbc4272539b7b6e43 *src/chm_common.c +32f9854050f40537e61497a181be2f80 *src/chm_common.c b6ea652eed1988974dec16b9cb9a35c1 *src/chm_common.h 4893273047274b21bfa8542880de3fa6 *src/cs.c 2d48f8b8f4e1eb4598302396f70e8c89 *src/cs.h -0a56e3408129fe35387a17bfd35e29a6 *src/cs_utils.c +749274eb9cfbbb5cf0db83c0ae4300cd *src/cs_utils.c 0b6112e86dab4fcececbadd2e84582ea *src/cs_utils.h -0f9ca4d63c86b9ae9299a9bdf4fad9c2 *src/dense.c +e1b9bfdf21d644dffd2b3764c6514cd6 *src/dense.c f14d4d8cd3b94ab5e95c8c9bff9630f1 *src/dense.h -c51ca6433a1a0361c3d95064a6b1690d *src/dgCMatrix.c +ad83ec224514ddc5304264281ea43450 *src/dgCMatrix.c 2392fd374cf7d259fedbef27f0a13e7f *src/dgCMatrix.h 62db9316297b07b1a82897e558b8e81a *src/dgTMatrix.c 341e051982a88f9c4d9a02c1376978fc *src/dgTMatrix.h @@ -433,7 +433,7 @@ ca00a8f2c785b0b4e3ead5cb888dc330 *src/dtpMatrix.h feba4a84f226246a992cc103f3da8874 *src/dtrMatrix.h cc1dca4dce1c8f68cf41f991ad349236 *src/factorizations.c 37fd937251ae9f1c68d974bd6169ade3 *src/factorizations.h -1e240fe50368743af3baa91f8ec289bb *src/init.c +069b9c048677ed614b7bb99a97609740 *src/init.c 5258096040c5fc0b0374fe510dc46c7d *src/ldense.c d8522bb27a18eea16d17c411d8da9e38 *src/ldense.h 0b0a474b5b19cb09637063134a0ba203 *src/lgCMatrix.c @@ -443,7 +443,7 @@ f3b6b0ff9a8088a21a1c46c9d7797f4d *src/scripts/DEPS.mkf 0aba26556ef5ddd3adc0031d41ef0d4c *src/scripts/DEPS.mkf_make.sh 15eb76d731d373a415312098f6f707e5 *src/scripts/SOURCES_C.mkf f1d325ea608f75622777bef1fbe55611 *src/scripts/fixup-fn.R -c4f73d8e14e8268b7ec4e7ce548e109e *src/sparseQR.c +80235a6bc16762bcd571ff7a9f553d14 *src/sparseQR.c 68c429cbb92387ef9ce70dac75a6e934 *src/sparseQR.h 2e623510d32213208f40240fc50a0c57 *src/t_Csparse_subassign.c 2394f5def013c7e571231935320fa6c3 *src/t_Csparse_validate.c @@ -452,7 +452,7 @@ c55fe12890cabcf61d562a4969825c02 *src/t_Matrix_rle.c edb93b23eb7e6db1f454970b4c663f3a *src/t_matrix_to_Csp.c e6de21ff3bf0908eaca620d9b868b4f1 *src/t_sparseVector.c 3dcc9418cfcbd6ec989b16a91b2b945a *tests/Class+Meth.R -25d02c4cc55e2d072606224fe229af9b *tests/Simple.R +b7de3c2add7cf86ac1e651697cd5f25b *tests/Simple.R 4966b9ac5a23de46814439a84174973b *tests/abIndex-tsts.R ce85ee2aba7d78e88fef738a0d1c7118 *tests/base-matrix-fun.R 78d7038fa6dfd7401a6a20efae58ab1e *tests/bind.R diff --git a/R/Auxiliaries.R b/R/Auxiliaries.R index a04bddd9..64ffe7da 100644 --- a/R/Auxiliaries.R +++ b/R/Auxiliaries.R @@ -69,7 +69,7 @@ is.null.DN <- function(dn) { ##' return 'x' unless it is NULL where you'd use 'orElse' `%||%` <- function(x, orElse) if(!is.null(x)) x else orElse -## not %in% : +##' not %in% : `%nin%` <- function (x, table) is.na(match(x, table)) nonTRUEoption <- function(ch) is.null(v <- getOption(ch)) || !isTRUE(v) diff --git a/R/Matrix.R b/R/Matrix.R index 3fa0007b..f4c69b7f 100644 --- a/R/Matrix.R +++ b/R/Matrix.R @@ -218,8 +218,7 @@ Matrix <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, } if(is.null(sparse1 <- sparse) && (i.M || is(data, "matrix"))) sparse <- sparseDefault(data) - doDN <- TRUE - i.m <- is.matrix(data) + doDN <- TRUE # by default if (i.M) { if (!sV) { if(!missing(nrow) || !missing(ncol)|| !missing(byrow)) @@ -230,7 +229,7 @@ Matrix <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, ## else : convert dense <-> sparse -> at end } } - else if (!i.m) { ## cut & paste from "base::matrix" : + else if(!is.matrix(data)) { ## cut & paste from "base::matrix" : ## avoid copying to strip attributes in simple cases if (is.object(data) || !is.atomic(data)) data <- as.vector(data) if(length(data) == 1 && is0(data) && !identical(sparse, FALSE)) { diff --git a/R/zzz.R b/R/zzz.R index 9fba7f2f..b96ff60e 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -59,6 +59,11 @@ if(getRversion() >= "3.2.0") { lengths <- function (x, use.names = TRUE) vapply(x, length, 1L, USE.NAMES = use.names) } +if(getRversion() < "3.5.0") { + isFALSE <- function (x) is.logical(x) && length(x) == 1L && !is.na(x) && !x + isTRUE <- function (x) is.logical(x) && length(x) == 1L && !is.na(x) && x +} + .onUnload <- function(libpath) { library.dynam.unload("Matrix", libpath) diff --git a/TODO b/TODO index ffd664a8..e3b3beb6 100644 --- a/TODO +++ b/TODO @@ -9,6 +9,8 @@ * *Urgent* in some sense --------------------------------------------------- ** TODO (partly DONE via workaround "round up" to 100): print() / show() for small options(max.print=): --> tests/Simple.R {'max.print'} +** TODO API change: Should Matrix(diag(2), sparse=TRUE, doDiag=TRUE) not rather give "ddiMatrix" ?? + Why change? Originally "ddiMatrix" etc extended denseMatrix but now sparseMatrix ** DONE qr.coef() has *wrong* (column)names, even in full-rank case: see man/qr-methods.Rd ("FIXME"); maybe related to ** DONE qr.R(), qrR() etc have *wrong* currently *lose* column/row names {compared to base R's qr.R}, see, drop0(R. <- qr.R(qx), tol=1e-15) # columns are int b1 c1 c2 b2 c3 {in man/qr-methods.Rd} diff --git a/build/partial.rdb b/build/partial.rdb index 19bdd6bb..fb57c8b2 100644 Binary files a/build/partial.rdb and b/build/partial.rdb differ diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 92db93a3..d9cd9fb2 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -14,6 +14,22 @@ %% \item \code{as(, "sparseMatrix")} no longer fails %% when \code{prod(dim(.))} is larger than \eqn{2^{31} - 1}. % +\section{Changes in version 1.2-17 (2019-03-.., svn r....)}{ + \subsection{New Features}{ + \itemize{ + \item (none) + } + } + \subsection{Bug Fixes}{ + \itemize{ + \item Fix new \code{PROTECT()} warnings (bugs?) found by \command{rchk}. + + \item Provide \code{isFALSE()} for \R < 3.5.0 as now need it for + sparseMatrix printing. + } + } +} + \section{Changes in version 1.2-16 (2019-03-04, svn r3291)}{ \subsection{New Features}{ \itemize{ diff --git a/inst/doc/Comparisons.pdf b/inst/doc/Comparisons.pdf index cec28afa..e68e4fbf 100644 Binary files a/inst/doc/Comparisons.pdf and b/inst/doc/Comparisons.pdf differ diff --git a/inst/doc/Design-issues.pdf b/inst/doc/Design-issues.pdf index 78330d51..a1208e30 100644 Binary files a/inst/doc/Design-issues.pdf and b/inst/doc/Design-issues.pdf differ diff --git a/inst/doc/Intro2Matrix.pdf b/inst/doc/Intro2Matrix.pdf index a0267f22..99b8f439 100644 Binary files a/inst/doc/Intro2Matrix.pdf and b/inst/doc/Intro2Matrix.pdf differ diff --git a/inst/doc/Introduction.pdf b/inst/doc/Introduction.pdf index 51809d3b..4a601642 100644 Binary files a/inst/doc/Introduction.pdf and b/inst/doc/Introduction.pdf differ diff --git a/inst/doc/sparseModels.pdf b/inst/doc/sparseModels.pdf index e2c2cc96..199e017b 100644 Binary files a/inst/doc/sparseModels.pdf and b/inst/doc/sparseModels.pdf differ diff --git a/man/nearPD.Rd b/man/nearPD.Rd index 21302b8b..e24a45c6 100644 --- a/man/nearPD.Rd +++ b/man/nearPD.Rd @@ -32,14 +32,14 @@ nearPD(x, corr = FALSE, keepDiag = FALSE, do2eigen = TRUE, direct fixpoint iteration \eqn{Y_k = P_U(P_S(Y_{k-1}))}{Y(k) = P_U(P_S(Y(k-1)))}.} \item{only.values}{logical; if \code{TRUE}, the result is just the - vector of eigen values of the approximating matrix.} + vector of eigenvalues of the approximating matrix.} \item{ensureSymmetry}{logical; by default, \code{\link{symmpart}(x)} is used whenever \code{\link{isSymmetric}(x)} is not true. The user can explicitly set this to \code{TRUE} or \code{FALSE}, saving the symmetry test. \emph{Beware} however that setting it \code{FALSE} for an \bold{a}symmetric input \code{x}, is typically nonsense!} \item{eig.tol}{defines relative positiveness of eigenvalues compared - to largest one, \eqn{\lambda_1}. Eigen values \eqn{\lambda_k} are + to largest one, \eqn{\lambda_1}. Eigenvalues \eqn{\lambda_k} are treated as if zero when \eqn{\lambda_k / \lambda_1 \le eig.tol}.} \item{conv.tol}{convergence tolerance for Higham algorithm.} \item{posd.tol}{tolerance for enforcing positive definiteness (in the @@ -55,27 +55,28 @@ nearPD(x, corr = FALSE, keepDiag = FALSE, do2eigen = TRUE, \details{ This implements the algorithm of Higham (2002), and then (if \code{do2eigen} is true) forces positive definiteness using code from - \code{\link[sfsmisc]{posdefify}}. The algorithm of Knol DL and ten - Berge (1989) (not implemented here) is more general in (1) that it - allows constraints to fix some rows (and columns) of the matrix and (2) - to force the smallest eigenvalue to have a certain value. + \code{\link[sfsmisc]{posdefify}}. The algorithm of Knol and ten + Berge (1989) (not implemented here) is more general in that it + allows constraints to (1) fix some rows (and columns) of the matrix and + (2) force the smallest eigenvalue to have a certain value. Note that setting \code{corr = TRUE} just sets \code{diag(.) <- 1} within the algorithm. Higham (2002) uses Dykstra's correction, but the version by Jens - Oehlschlaegel did not use it (accidentally), and has still lead to - good results; this simplification, now only via \code{doDykstra = FALSE}, - was active in \code{nearPD()} upto Matrix version 0.999375-40. + Oehlschlaegel did not use it (accidentally), and still gave + reasonable results; this simplification, now only + used if \code{doDykstra = FALSE}, + was active in \code{nearPD()} up to Matrix version 0.999375-40. } \value{ - If \code{only.values = TRUE}, a numeric vector of eigen values of the + If \code{only.values = TRUE}, a numeric vector of eigenvalues of the approximating matrix; Otherwise, as by default, an S3 object of \code{\link{class}} \code{"nearPD"}, basically a list with components \item{mat}{a matrix of class \code{\linkS4class{dpoMatrix}}, the computed positive-definite matrix.} - \item{eigenvalues}{numeric vector of eigen values of \code{mat}.} + \item{eigenvalues}{numeric vector of eigenvalues of \code{mat}.} \item{corr}{logical, just the argument \code{corr}.} \item{normF}{the Frobenius norm (\code{\link{norm}(x-X, "F")}) of the difference between the original and the resulting matrix.} diff --git a/src/Csparse.c b/src/Csparse.c index 7162dd52..35999c6f 100644 --- a/src/Csparse.c +++ b/src/Csparse.c @@ -163,8 +163,8 @@ SEXP Csparse_to_dense(SEXP x, SEXP symm_or_tri) /* transp: */ FALSE); // -> a [dln]geMatrix if(is_sym) { // ==> want [dln]syMatrix - const char cl1 = class_P(ans)[0]; PROTECT(ans); + const char cl1 = class_P(ans)[0]; SEXP aa = PROTECT(NEW_OBJECT_OF_CLASS((cl1 == 'd') ? "dsyMatrix" : ((cl1 == 'l') ? "lsyMatrix" : "nsyMatrix"))); // No need to duplicate() as slots of ans are freshly allocated and ans will not be used @@ -176,8 +176,8 @@ SEXP Csparse_to_dense(SEXP x, SEXP symm_or_tri) return aa; } else if(is_tri) { // ==> want [dln]trMatrix - const char cl1 = class_P(ans)[0]; PROTECT(ans); + const char cl1 = class_P(ans)[0]; SEXP aa = PROTECT(NEW_OBJECT_OF_CLASS((cl1 == 'd') ? "dtrMatrix" : ((cl1 == 'l') ? "ltrMatrix" : "ntrMatrix"))); // No need to duplicate() as slots of ans are freshly allocated and ans will not be used diff --git a/src/Mutils.h b/src/Mutils.h index 87120077..4d0a5dee 100644 --- a/src/Mutils.h +++ b/src/Mutils.h @@ -300,10 +300,10 @@ Rboolean any_NA_in_x(SEXP obj) static R_INLINE SEXP inv_permutation(SEXP p_, SEXP zero_p, SEXP zero_res) { - int np = 0; + int np = 1; if(!isInteger(p_)) {p_ = PROTECT(coerceVector(p_, INTSXP)); np++; } int *p = INTEGER(p_), n = LENGTH(p_); - SEXP val = allocVector(INTSXP, n);// (not PROTECT()ing: no alloc from here on) + SEXP val = PROTECT(allocVector(INTSXP, n)); int *v = INTEGER(val), p_0 = asLogical(zero_p), r_0 = asLogical(zero_res); if(!p_0) v--; // ==> use 1-based indices // shorter (but not 100% sure if ok: is LHS always eval'ed *before* RHS ?) : diff --git a/src/Syms.h b/src/Syms.h index dbb17311..dfdcd447 100644 --- a/src/Syms.h +++ b/src/Syms.h @@ -7,6 +7,7 @@ SEXP Matrix_iSym, Matrix_jSym, Matrix_lengthSym, + Matrix_LSym, Matrix_RSym, Matrix_USym, Matrix_pSym, Matrix_permSym, Matrix_uploSym, diff --git a/src/chm_common.c b/src/chm_common.c index 1992b131..dd498ea1 100644 --- a/src/chm_common.c +++ b/src/chm_common.c @@ -48,27 +48,33 @@ void CHM_store_common() { } void CHM_restore_common() { - SEXP rho = chm_common_env; - c.dbound = asReal(findVarInFrame(rho, dboundSym)); - c.grow0 = asReal(findVarInFrame(rho, grow0Sym)); - c.grow1 = asReal(findVarInFrame(rho, grow1Sym)); - c.grow2 = asInteger(findVarInFrame(rho, grow2Sym)); - c.maxrank = asInteger(findVarInFrame(rho, maxrankSym)); - c.supernodal_switch = asReal(findVarInFrame(rho, supernodal_switchSym)); - c.supernodal = asLogical(findVarInFrame(rho, supernodalSym)); - c.final_asis = asLogical(findVarInFrame(rho, final_asisSym)); - c.final_super = asLogical(findVarInFrame(rho, final_superSym)); - c.final_ll = asLogical(findVarInFrame(rho, final_llSym)); - c.final_pack = asLogical(findVarInFrame(rho, final_packSym)); - c.final_monotonic = asLogical(findVarInFrame(rho, final_monotonicSym)); - c.final_resymbol = asLogical(findVarInFrame(rho, final_resymbolSym)); - c.prefer_zomplex = asLogical(findVarInFrame(rho, prefer_zomplexSym)); - c.prefer_upper = asLogical(findVarInFrame(rho, prefer_upperSym)); - c.quick_return_if_not_posdef = - asLogical(findVarInFrame(rho, quick_return_if_not_posdefSym)); - c.nmethods = asInteger(findVarInFrame(rho, nmethodsSym)); - c.method[0].ordering = asInteger(findVarInFrame(rho, m0_ordSym)); - c.postorder = asLogical(findVarInFrame(rho, postorderSym)); + SEXP rho = chm_common_env, var; + +#define SET_AS_FROM_FRAME(_V_, _KIND_, _SYM_) \ + var = PROTECT(findVarInFrame(rho, _SYM_)); \ + _V_ = _KIND_(var); \ + UNPROTECT(1) + + SET_AS_FROM_FRAME(c.dbound, asReal, dboundSym); + SET_AS_FROM_FRAME(c.grow0, asReal, grow0Sym); + SET_AS_FROM_FRAME(c.grow1, asReal, grow1Sym); + SET_AS_FROM_FRAME(c.grow2, asInteger, grow2Sym); + SET_AS_FROM_FRAME(c.maxrank,asInteger, maxrankSym); + SET_AS_FROM_FRAME(c.supernodal_switch, asReal, supernodal_switchSym); + SET_AS_FROM_FRAME(c.supernodal, asLogical, supernodalSym); + SET_AS_FROM_FRAME(c.final_asis, asLogical, final_asisSym); + SET_AS_FROM_FRAME(c.final_super, asLogical, final_superSym); + SET_AS_FROM_FRAME(c.final_ll, asLogical, final_llSym); + SET_AS_FROM_FRAME(c.final_pack, asLogical, final_packSym); + SET_AS_FROM_FRAME(c.final_monotonic,asLogical, final_monotonicSym); + SET_AS_FROM_FRAME(c.final_resymbol, asLogical, final_resymbolSym); + SET_AS_FROM_FRAME(c.prefer_zomplex, asLogical, prefer_zomplexSym); + SET_AS_FROM_FRAME(c.prefer_upper, asLogical, prefer_upperSym); + SET_AS_FROM_FRAME(c.quick_return_if_not_posdef, + asLogical, quick_return_if_not_posdefSym); + SET_AS_FROM_FRAME(c.nmethods, asInteger, nmethodsSym); + SET_AS_FROM_FRAME(c.method[0].ordering, asInteger, m0_ordSym); + SET_AS_FROM_FRAME(c.postorder, asLogical, postorderSym); } SEXP CHM_set_common_env(SEXP rho) { diff --git a/src/cs_utils.c b/src/cs_utils.c index e7aae357..8f41ab2b 100644 --- a/src/cs_utils.c +++ b/src/cs_utils.c @@ -207,8 +207,8 @@ csn *Matrix_as_csn(csn *ans, SEXP x) if (ctype < 0) error(_("invalid class of object to %s"), "Matrix_as_csn"); - ans->U = Matrix_as_cs(GET_SLOT(x, install("U"))); - ans->L = Matrix_as_cs(GET_SLOT(x, install("L"))); + ans->U = Matrix_as_cs(GET_SLOT(x, Matrix_USym)); + ans->L = Matrix_as_cs(GET_SLOT(x, Matrix_LSym)); switch(ctype) { case 0: ans->B = (double*) NULL; @@ -294,9 +294,9 @@ SEXP Matrix_csn_to_SEXP(csn *N, char *cl, int dofree, SEXP dn) ans = PROTECT(NEW_OBJECT_OF_CLASS(cl)); /* allocate and copy common slots */ /* FIXME: Use the triangular matrix classes for csn_LU */ - SET_SLOT(ans, install("L"), /* these are free'd later if requested */ + SET_SLOT(ans, Matrix_LSym, /* these are free'd later if requested */ Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // FIXME: dn - SET_SLOT(ans, install("U"), + SET_SLOT(ans, Matrix_USym, Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn)); // FIXME: dn switch(ctype) { case 0: diff --git a/src/dense.c b/src/dense.c index 2c8902c4..83ce0e8e 100644 --- a/src/dense.c +++ b/src/dense.c @@ -22,24 +22,23 @@ static int left_cyclic(double x[], int ldx, int j, int k, double cosines[], double sines[]) { - double *lastcol; - int i, jj; - if (j >= k) error(_("incorrect left cyclic shift, j (%d) >= k (%d)"), j, k); if (j < 0) error(_("incorrect left cyclic shift, j (%d) < 0"), j, k); if (ldx < k) error(_("incorrect left cyclic shift, k (%d) > ldx (%d)"), k, ldx); - lastcol = (double*) R_alloc(k+1, sizeof(double)); + + double *lastcol = (double*) R_alloc(k+1, sizeof(double)); + int i; /* keep a copy of column j */ for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx]; /* For safety, zero the rest */ for(i = j+1; i <= k; i++) lastcol[i] = 0.; - for(jj = j+1; jj <= k; jj++) { /* columns to be shifted */ - int diagind = jj*(ldx+1), ind = (jj-j) - 1; + for(int jj = j+1, ind = 0; jj <= k; jj++, ind++) { /* columns to be shifted */ + int diagind = jj*(ldx+1); // ind == (jj-j) - 1 double tmp = x[diagind], cc, ss; - /* Calculate the Givens rotation. */ + /* Calculate the Givens rotation. */ /* This modified the super-diagonal element */ F77_CALL(drotg)(x + diagind-1, &tmp, cosines + ind, sines + ind); cc = cosines[ind]; ss = sines[ind]; @@ -68,7 +67,7 @@ SEXP getGivens(double x[], int ldx, int jmin, int rank) SET_VECTOR_ELT(ans, 0, ScalarInteger(jmin)); SET_VECTOR_ELT(ans, 1, ScalarInteger(rank)); SET_VECTOR_ELT(ans, 2, cosines = allocVector(REALSXP, shiftlen)); - SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen)); + SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen)); setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 4)); SET_STRING_ELT(nms, 0, mkChar("jmin")); SET_STRING_ELT(nms, 1, mkChar("rank")); @@ -99,29 +98,29 @@ SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank) SEXP lsq_dense_Chol(SEXP X, SEXP y) { SEXP ans; - int info, n, p, k, *Xdims, *ydims; - double *xpx, d_one = 1., d_zero = 0.; + double d_one = 1., d_zero = 0.; if (!(isReal(X) & isMatrix(X))) error(_("X must be a numeric (double precision) matrix")); - Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP)); - n = Xdims[0]; - p = Xdims[1]; + int *Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP)), + n = Xdims[0], + p = Xdims[1]; if (!(isReal(y) & isMatrix(y))) error(_("y must be a numeric (double precision) matrix")); - ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP)); + int *ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP)); if (ydims[0] != n) error(_( "number of rows in y (%d) does not match number of rows in X (%d)"), ydims[0], n); - k = ydims[1]; + int k = ydims[1]; if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k); ans = PROTECT(allocMatrix(REALSXP, p, k)); F77_CALL(dgemm)("T", "N", &p, &k, &n, &d_one, REAL(X), &n, REAL(y), &n, &d_zero, REAL(ans), &p); - xpx = (double *) R_alloc(p * p, sizeof(double)); + double *xpx = (double *) R_alloc(p * p, sizeof(double)); F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero, xpx, &p); + int info; F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info); if (info) error(_("Lapack routine dposv returned error code %d"), info); UNPROTECT(1); @@ -131,34 +130,30 @@ SEXP lsq_dense_Chol(SEXP X, SEXP y) SEXP lsq_dense_QR(SEXP X, SEXP y) { - SEXP ans; - int info, n, p, k, *Xdims, *ydims, lwork; - double *work, tmp, *xvals; - if (!(isReal(X) & isMatrix(X))) error(_("X must be a numeric (double precision) matrix")); - Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP)); - n = Xdims[0]; - p = Xdims[1]; + int *Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP)), + n = Xdims[0], + p = Xdims[1]; if (!(isReal(y) & isMatrix(y))) error(_("y must be a numeric (double precision) matrix")); - ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP)); + int *ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP)); if (ydims[0] != n) error(_( "number of rows in y (%d) does not match number of rows in X (%d)"), ydims[0], n); - k = ydims[1]; + int k = ydims[1]; if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k); - xvals = (double *) Memcpy(R_alloc(n * p, sizeof(double)), REAL(X), n * p); - ans = PROTECT(duplicate(y)); - lwork = -1; + double tmp, *xvals = (double *) Memcpy(R_alloc(n * p, sizeof(double)), REAL(X), n * p); + SEXP ans = PROTECT(duplicate(y)); + int lwork = -1, info; F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n, &tmp, &lwork, &info); if (info) error(_("First call to Lapack routine dgels returned error code %d"), info); lwork = (int) tmp; - work = (double *) R_alloc(lwork, sizeof(double)); + double *work = (double *) R_alloc(lwork, sizeof(double)); F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n, work, &lwork, &info); if (info) @@ -448,14 +443,13 @@ SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test) SEXP ddense_symmpart(SEXP x) /* Class of the value will be dsyMatrix */ { - SEXP dx = dup_mMatrix_as_dgeMatrix(x); + SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x)); int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0]; if(n != adims[1]) { error(_("matrix is not square! (symmetric part)")); return R_NilValue; /* -Wall */ } else { - PROTECT(dx); SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("dsyMatrix")), dns, nms_dns; double *xx = REAL(GET_SLOT(dx, Matrix_xSym)); diff --git a/src/dgCMatrix.c b/src/dgCMatrix.c index 9571c347..727cb61a 100644 --- a/src/dgCMatrix.c +++ b/src/dgCMatrix.c @@ -284,7 +284,8 @@ SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames) } } else ALLOC_SLOT(ans, install("q"), INTSXP, 0); - SET_SLOT(ans, install("R"), Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn)); + SEXP R = PROTECT(Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn)); + SET_SLOT(ans, Matrix_RSym, R); UNPROTECT(1); // R if(do_dn) UNPROTECT(1); // dn cs_nfree(N); cs_sfree(S); @@ -338,7 +339,7 @@ SEXP dgCMatrix_SPQR(SEXP Ap, SEXP ordering, SEXP econ, SEXP tol) * may make sense if to be used in the "spqr_solve" routines .. ?? */ /* SET_VECTOR_ELT(ans, 1, */ /* chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); */ - SET_SLOT(ans, install("R"), + SET_SLOT(ans, Matrix_RSym, chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); cholmod_free_sparse(&Al, &cl); cholmod_free_sparse(&R, &cl); @@ -415,7 +416,7 @@ void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep SET_VECTOR_ELT(dn, 1, R_NilValue); // colnames(.) := NULL } } - SET_SLOT(ans, install("L"), + SET_SLOT(ans, Matrix_LSym, Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0, do_dn ? dn : R_NilValue)); if(keep_dimnms) { @@ -435,7 +436,7 @@ void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep SET_VECTOR_ELT(dn, 0, R_NilValue); // rownames(.) := NULL } } - SET_SLOT(ans, install("U"), + SET_SLOT(ans, Matrix_USym, Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0, do_dn ? dn : R_NilValue)); if(do_dn) UNPROTECT(1); // dn Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */ @@ -501,8 +502,8 @@ SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse) lu = get_factors(Ap, "LU"); } qslot = GET_SLOT(lu, install("q")); - L = AS_CSP__(GET_SLOT(lu, install("L"))); - U = AS_CSP__(GET_SLOT(lu, install("U"))); + L = AS_CSP__(GET_SLOT(lu, Matrix_LSym)); + U = AS_CSP__(GET_SLOT(lu, Matrix_USym)); R_CheckStack(); if (U->n != n) error(_("Dimensions of system to be solved are inconsistent")); diff --git a/src/init.c b/src/init.c index 1a87da3d..3d397e44 100644 --- a/src/init.c +++ b/src/init.c @@ -365,6 +365,9 @@ R_init_Matrix(DllInfo *dll) Matrix_permSym = install("perm"); Matrix_uploSym = install("uplo"); Matrix_xSym = install("x"); + Matrix_LSym = install("L"); + Matrix_RSym = install("R"); + Matrix_USym = install("U"); Matrix_VSym = install("V"); Matrix_NS = R_FindNamespace(mkString("Matrix")); diff --git a/src/sparseQR.c b/src/sparseQR.c index 8a685d1a..5ae73fab 100644 --- a/src/sparseQR.c +++ b/src/sparseQR.c @@ -3,7 +3,7 @@ SEXP sparseQR_validate(SEXP x) { CSP V = AS_CSP__(GET_SLOT(x, Matrix_VSym)), - R = AS_CSP__(GET_SLOT(x, install("R"))); + R = AS_CSP__(GET_SLOT(x, Matrix_RSym)); SEXP beta = GET_SLOT(x, Matrix_betaSym), p = GET_SLOT(x, Matrix_pSym), q = GET_SLOT(x, install("q")); @@ -145,7 +145,7 @@ SEXP sparseQR_qty(SEXP qr, SEXP y, SEXP trans, SEXP keep_dimnames) // Compute qr.coef(qr, y) := R^{-1} Q' y {modulo row and column permutations} SEXP sparseQR_coef(SEXP qr, SEXP y) { - SEXP qslot = GET_SLOT(qr, install("q")), R_ = GET_SLOT(qr, install("R")); + SEXP qslot = GET_SLOT(qr, install("q")), R_ = GET_SLOT(qr, Matrix_RSym); CSP R = AS_CSP__(R_); // FIXME: check n_R, M (= R->m) vs n, m int *q = INTEGER(qslot), lq = LENGTH(qslot), n_R = R->n; // = ncol(R) diff --git a/tests/Simple.R b/tests/Simple.R index dd5762ff..eb11d17f 100644 --- a/tests/Simple.R +++ b/tests/Simple.R @@ -31,18 +31,25 @@ if(interactive()) { # ^^^^^^ to show Matrix.msg()s ### Matrix() ''smartness'' -(d4 <- d40 <- Matrix(diag(4))) +(d40 <- Matrix( diag(4))) (z4 <- Matrix(0*diag(4))) (o4 <- Matrix(1+diag(4))) (tr <- Matrix(cbind(1,0:1))) (M4 <- Matrix(m4 <- cbind(0,rbind(6*diag(3),0)))) dM4 <- Matrix(M4, sparse = FALSE) d4. <- diag(4); dimnames(d4.) <- dns <- rep(list(LETTERS[1:4]), 2) +d4a <- diag(4); dimnames(d4a) <- dna <- list(LETTERS[1:4], letters[1:4])# "a"symmetric +m1a <- matrix(0, dimnames=list("A","b"))# "a"symmetric d4di<- as(d4., "diagonalMatrix") +d4da<- as(d4a, "diagonalMatrix") d4d <- as(d4., "denseMatrix") +d4aS <- Matrix(d4a, sparse=TRUE, doDiag=FALSE) +d1aS <- Matrix(m1a, sparse=TRUE, doDiag=FALSE) stopifnot(identical(d4di@x, numeric()), # was "named" unnecessarily - identical(dimnames(d4 <- Matrix(d4.)), dns), identical(unname(d4), d40), + identical(dimnames(d4 <- Matrix(d4.)), dns), + identical4(d40, Matrix(diag(4)), unname(d4), unname(d4da)), identical3(d4, as(d4., "Matrix"), as(d4., "diagonalMatrix")), + is(d4aS, "dtCMatrix"), # not "dsC*", as asymmetric dimnames is(d4d, "denseMatrix")) class(mN <- Matrix(NA, 3,4)) # NA *is* logical @@ -69,6 +76,7 @@ stopifnot(all(is.na(sL@x)), ## not yet: all(is.na(sL)), validObject(Matrix(c(NA,0), 4, 4))) stopifnotValid(Matrix(c(NA,0,0,0), 4, 4), "sparseMatrix") I <- i1 <- I1 <- Diagonal(1) +## TODO? stopifnot(identical(I, Matrix(1, sparse=TRUE))) # doDiag=TRUE default I1[1,1] <- i1[1, ] <- I [ ,1] <- NA stopifnot(identical3(I,i1,I1)) image(d4) # gave infinite recursion