diff -Nru rmatrix-1.1-3/debian/changelog rmatrix-1.1-4/debian/changelog --- rmatrix-1.1-3/debian/changelog 2014-06-21 12:31:13.000000000 +0000 +++ rmatrix-1.1-4/debian/changelog 2014-06-21 12:31:13.000000000 +0000 @@ -1,8 +1,16 @@ -rmatrix (1.1-3-1precise0) precise; urgency=low +rmatrix (1.1-4-1precise0) precise; urgency=low * Compilation for Ubuntu 12.04.4 LTS - -- Michael Rutter Tue, 01 Apr 2014 02:55:36 +0000 + -- Michael Rutter Sat, 21 Jun 2014 12:29:10 +0000 + +rmatrix (1.1-4-1) unstable; urgency=low + + * New upstream release + + * debian/control: Set Build-Depends: to current R version + + -- Dirk Eddelbuettel Sun, 15 Jun 2014 15:42:27 -0500 rmatrix (1.1-3-1) unstable; urgency=low diff -Nru rmatrix-1.1-3/debian/control rmatrix-1.1-4/debian/control --- rmatrix-1.1-3/debian/control 2014-06-21 12:31:13.000000000 +0000 +++ rmatrix-1.1-4/debian/control 2014-06-21 12:31:13.000000000 +0000 @@ -2,7 +2,7 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper (>= 7.0.0), cdbs, r-base-dev (>= 3.0.3), r-cran-lattice (>= 0.12-11.1) +Build-Depends: debhelper (>= 7.0.0), cdbs, r-base-dev (>= 3.1.0), r-cran-lattice (>= 0.12-11.1) Standards-Version: 3.9.5 Package: r-cran-matrix diff -Nru rmatrix-1.1-3/DESCRIPTION rmatrix-1.1-4/DESCRIPTION --- rmatrix-1.1-3/DESCRIPTION 2014-03-31 10:51:10.000000000 +0000 +++ rmatrix-1.1-4/DESCRIPTION 2014-06-15 18:26:41.000000000 +0000 @@ -1,6 +1,6 @@ Package: Matrix -Version: 1.1-3 -Date: 2014-03-26 +Version: 1.1-4 +Date: 2014-06-14 Priority: recommended Title: Sparse and Dense Matrix Classes and Methods Author: Douglas Bates and Martin Maechler @@ -23,7 +23,7 @@ Davis. All sections of that code are covered by the GPL or LGPL licenses. See the directory doc/UFsparse for details. URL: http://Matrix.R-forge.R-project.org/ -Packaged: 2014-03-30 16:34:00 UTC; maechler +Packaged: 2014-06-14 16:22:31 UTC; maechler NeedsCompilation: yes Repository: CRAN -Date/Publication: 2014-03-31 12:51:10 +Date/Publication: 2014-06-15 20:26:41 Binary files /tmp/dq1Sb5R1g0/rmatrix-1.1-3/inst/doc/Comparisons.pdf and /tmp/1Zkg5R_7HI/rmatrix-1.1-4/inst/doc/Comparisons.pdf differ Binary files /tmp/dq1Sb5R1g0/rmatrix-1.1-3/inst/doc/Design-issues.pdf and /tmp/1Zkg5R_7HI/rmatrix-1.1-4/inst/doc/Design-issues.pdf differ Binary files /tmp/dq1Sb5R1g0/rmatrix-1.1-3/inst/doc/Intro2Matrix.pdf and /tmp/1Zkg5R_7HI/rmatrix-1.1-4/inst/doc/Intro2Matrix.pdf differ Binary files /tmp/dq1Sb5R1g0/rmatrix-1.1-3/inst/doc/Introduction.pdf and /tmp/1Zkg5R_7HI/rmatrix-1.1-4/inst/doc/Introduction.pdf differ Binary files /tmp/dq1Sb5R1g0/rmatrix-1.1-3/inst/doc/sparseModels.pdf and /tmp/1Zkg5R_7HI/rmatrix-1.1-4/inst/doc/sparseModels.pdf differ diff -Nru rmatrix-1.1-3/inst/test-tools-Matrix.R rmatrix-1.1-4/inst/test-tools-Matrix.R --- rmatrix-1.1-3/inst/test-tools-Matrix.R 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/inst/test-tools-Matrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -206,7 +206,9 @@ r } -mkLDL <- function(n, density = 1/3) { +mkLDL <- function(n, density = 1/3, + d0 = 10, ## use '10' to be "different" from L entries + rcond = (n < 99), condest = (n >= 100)) { ## Purpose: make nice artificial A = L D L' (with exact numbers) decomp ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 15 Mar 2008 @@ -218,11 +220,13 @@ L[sample(n*n, nnz)] <- seq_len(nnz) L <- tril(L,-1) diag(L) <- 1 - d.half <- sample(10*(n:1))# random permutation ; use '10*' to be "different" from L entries + d.half <- d0 * sample.int(n)# random permutation D <- Diagonal(x = d.half * d.half) A <- tcrossprod(L * rep(d.half, each=n)) ## = as(L %*% D %*% t(L), "symmetricMatrix") - list(A = A, L = L, d.half = d.half, D = D) + list(A = A, L = L, d.half = d.half, D = D, + rcond.A = if(rcond) rcond(A, useInv=TRUE), + cond.A = if(condest) condest(A)$est) } eqDeterminant <- function(m1, m2, NA.Inf.ok=FALSE, tol=.Machine$double.eps^0.5, ...) diff -Nru rmatrix-1.1-3/man/CHMfactor-class.Rd rmatrix-1.1-4/man/CHMfactor-class.Rd --- rmatrix-1.1-3/man/CHMfactor-class.Rd 2014-03-11 18:33:16.000000000 +0000 +++ rmatrix-1.1-4/man/CHMfactor-class.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -12,6 +12,7 @@ % \alias{coerce,CHMfactor,Matrix-method} \alias{coerce,CHMfactor,sparseMatrix-method} +\alias{coerce,CHMfactor,triangularMatrix-method} \alias{coerce,CHMfactor,pMatrix-method} \alias{expand,CHMfactor-method} %\alias{solve,CHMfactor,...}%--> solve-methods.Rd @@ -97,20 +98,22 @@ \item{isLDL}{\code{(x)} returns a \code{\link{logical}} indicating if \code{x} is an \eqn{LDL'} decomposition or (when \code{FALSE}) an \eqn{LL'} one.} - \item{coerce}{\code{signature(from = "CHMfactor", to = "sparseMatrix")} - Returns the lower triangular factor \eqn{L} from the \eqn{LL'} - form of the Cholesky factorization. Note that (currently) the - factor from the \eqn{LL'} form is always returned, even if the - \code{"CHMfactor"} object represents an \eqn{LDL'} decomposition. + \item{coerce}{\code{signature(from = "CHMfactor", to = "sparseMatrix")} + (or equivalently, \code{to = "Matrix"} or \code{to = "triangularMatrix"}) + + \code{as(*, "sparseMatrix")} returns the lower triangular factor + \eqn{L} from the \eqn{LL'} form of the Cholesky factorization. + Note that (currently) the factor from the \eqn{LL'} form is always + returned, even if the \code{"CHMfactor"} object represents an + \eqn{LDL'} decomposition. Furthermore, this is the factor after any fill-reducing permutation has been applied. See the \code{expand} method for obtaining both the permutation matrix, \eqn{P}, and the lower Cholesky factor, \eqn{L}.} - \item{coerce}{\code{signature(from = "CHMfactor", to = - "pMatrix")} Returns the permutation matrix, \eqn{P}, - representing the fill-reducing permutation used in the - decomposition.} - \item{expand}{\code{signature(x = "CHMfactor")} Returns a list with + \item{coerce}{\code{signature(from = "CHMfactor", to = "pMatrix")} + returns the permutation matrix \eqn{P}, representing the + fill-reducing permutation used in the decomposition.} + \item{expand}{\code{signature(x = "CHMfactor")} returns a list with components \code{P}, the matrix representing the fill-reducing permutation, and \code{L}, the lower triangular Cholesky factor. The original positive-definite matrix \eqn{A} corresponds to the product @@ -118,7 +121,7 @@ product may apparently have more non-zeros than the original matrix, even after applying \code{\link{drop0}} to it. However, the extra "non-zeros" should be very small in magnitude.} - \item{image}{\code{signature(x = "CHMfactor")} Plot the image of the + \item{image}{\code{signature(x = "CHMfactor"):} Plot the image of the lower triangular factor, \eqn{L}, from the decomposition. This method is equivalent to \code{image(as(x, "sparseMatrix"))} so the comments in the above description of the \code{coerce} method @@ -135,7 +138,7 @@ to produce \code{a}. Analogously, \code{system = "L"} returns the solution \eqn{x}, of \eqn{Lx = b}. Similarly, for all system codes \bold{but} \code{"P"} and \code{"Pt"} - where, e.g., \code{x <- solve(a, b,system="P")} is equivalent to + where, e.g., \code{x <- solve(a, b, system="P")} is equivalent to \code{x <- P \%*\% b}. See also \code{\link{solve-methods}}. diff -Nru rmatrix-1.1-3/man/Cholesky.Rd rmatrix-1.1-4/man/Cholesky.Rd --- rmatrix-1.1-3/man/Cholesky.Rd 2014-02-28 22:50:04.000000000 +0000 +++ rmatrix-1.1-4/man/Cholesky.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -7,12 +7,13 @@ \alias{Cholesky,Matrix-method}% <- "good" error message \alias{.SuiteSparse_version} \title{Cholesky Decomposition of a Sparse Matrix} +\concept{Choleski}% alternative English form \usage{ Cholesky(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, \dots) } \description{ - Computes the Cholesky decomposition of a sparse, symmetric, - positive-definite matrix. However, typically \code{\link{chol}()} + Computes the Cholesky (aka \dQuote{Choleski}) decomposition of a + sparse, symmetric, positive-definite matrix. However, typically \code{\link{chol}()} should rather be used unless you are interested in the different kinds of sparse Cholesky decompositions. } @@ -45,8 +46,9 @@ In other words, the result of \code{Cholesky()} is \emph{not} a matrix, and if you want one, you should probably rather use - \code{\link{chol}()}. + \code{\link{chol}()}, see Details. } +%FAILS WHY??\newcommand{\tR}{\tilde{R}}% both \tR{} and simple '\tR' fails \details{ This is a generic function with special methods for different types of matrices. Use \code{\link{showMethods}("Cholesky")} to list all @@ -56,8 +58,20 @@ --- the only one available currently --- is based on functions from the CHOLMOD library. - Again: If you just want the Cholesky decomposition of a matrix, you - should probably rather use \code{\link{chol}(.)}. + Again: If you just want the Cholesky decomposition of a matrix in a + straightforward way, you should probably rather use \code{\link{chol}(.)}. + + Note that if \code{perm=TRUE} (default), the decomposition is + \deqn{A = P' \tilde{L} D \tilde{L}' P = P' L L' P,}{A = P' L~ D L~' P = P' L L' P,} + where \eqn{L} can be extracted by \code{as(*, "Matrix")}, \eqn{P} by + \code{as(*, "pMatrix")} and both by \code{\link{expand}(*)}, see the + class \code{\linkS4class{CHMfactor}} documentation. + + Note that consequently, you cannot easily get the \dQuote{traditional} + cholesky factor \eqn{R}, from this decomposition, as + \deqn{R'R = A = P'LL'P = P'\tilde{R}'\tilde{R} P = (\tilde{R}P)' (\tilde{R}P),}{ + R'R = A = P'LL'P = P' R~' R~ P = (R~ P)' (R~ P),} + but \eqn{\tilde{R}P}{R~ P} is \emph{not} triangular even though \eqn{\tilde{R}}{R~} is. } \references{ Tim Davis (2005) diff -Nru rmatrix-1.1-3/man/chol.Rd rmatrix-1.1-4/man/chol.Rd --- rmatrix-1.1-3/man/chol.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/chol.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -27,16 +27,33 @@ \usage{ chol(x, \dots) \S4method{chol}{dsCMatrix}(x, pivot = FALSE, \dots) +\S4method{chol}{dsparseMatrix}(x, pivot = FALSE, cache = TRUE, \dots) } \arguments{ \item{x}{a (sparse or dense) square matrix, here inheriting from class \code{\linkS4class{Matrix}}; if \code{x} is not positive definite, an error is signalled.} - \item{pivot}{logical indicating if pivoting is to be used.} + \item{pivot}{logical indicating if pivoting is to be used. Currently, + this is \emph{not} made use of for dense matrices.} + \item{cache}{logical indicating if the result should be cashed in + \code{x@factors}; note that this argument is experimental and only + available for some sparse matrices.} \item{\dots}{potentially further arguments passed to methods.} } -% \details{ -% } +\details{ + Note that these Cholesky factorizations are typically \emph{cached} with + \code{x} currently, and these caches are available in + \code{x@factors}, which may be useful for the sparse case when + \code{pivot = TRUE}, where the permutation can be retrieved; see also + the examples. + + However, this should not be considered part of the API and made + use of. Rather consider \code{\link{Cholesky}()} in such situations, + since \code{chol(x, pivot=TRUE)} uses the same algorithm (but not the + same return value!) as \code{\link{Cholesky}(x, LDL=FALSE)} and + \code{chol(x)} corresponds to + \code{\link{Cholesky}(x, perm=FALSE, LDL=FALSE)}. +} \section{Methods}{ Use \code{\link{showMethods}(chol)} to see all; some are worth mentioning here: @@ -54,14 +71,16 @@ \item{chol}{\code{signature(x = "dsCMatrix", pivot = "logical")}: Returns (and stores) the Cholesky decomposition of \code{x}. If - \code{pivot} is \code{TRUE} (the default) the Approximate Minimal - Degree (AMD) algorithm is used to create a reordering of the rows - and columns of \code{x} so as to reduce fill-in.} + \code{pivot} is true, the Approximate Minimal Degree (AMD) + algorithm is used to create a reordering of the rows and columns + of \code{x} so as to reduce fill-in.} } } \value{ a matrix of class \code{\linkS4class{Cholesky}}, - i.e., upper triangular: \eqn{R} such that \eqn{R'R = x}. + i.e., upper triangular: \eqn{R} such that \eqn{R'R = x} (if + \code{pivot=FALSE}) \emph{or} \eqn{P' R'R P = x} (if + \code{pivot=TRUE} and \eqn{P} is the corresponding permutation matrix). } \references{ Timothy A. Davis (2006) @@ -73,7 +92,8 @@ \emph{SIAM J. Matrix Analysis and Applications}, \bold{17}, 4, 886--905. } -\seealso{The default from \pkg{base}, \code{\link[base]{chol}}. +\seealso{The default from \pkg{base}, \code{\link[base]{chol}}; for more + flexibility (but not returning a matrix!) \code{\link{Cholesky}}. } \examples{ showMethods(chol, inherited = FALSE) # show different methods @@ -103,11 +123,16 @@ (M2 <- toeplitz(as(c(1,.5, rep(0,12), -.1), "sparseVector"))) c2 <- chol(M2) C2 <- chol(M2, pivot=TRUE) -## check the caching of the factorizations: -spd <- as(M2@factors[["spdCholesky"]], "Matrix") -sPd <- as(M2@factors[["sPdCholesky"]], "Matrix") -stopifnot(identical(t(spd), c2), - all.equal(t(sPd), C2, tolerance=0))#%% FIXME -- why not identical()? +## For the experts, check the caching of the factorizations: +ff <- M2@factors[["spdCholesky"]] +FF <- M2@factors[["sPdCholesky"]] +L1 <- as(ff, "Matrix")# pivot=FALSE: no perm. +L2 <- as(FF, "Matrix"); P2 <- as(FF, "pMatrix") +stopifnot(identical(t(L1), c2), + all.equal(t(L2), C2, tolerance=0),#-- why not identical()? + all.equal(M2, tcrossprod(L1)), # M = LL' + all.equal(M2, crossprod(crossprod(L2, P2)))# M = P'L L'P + ) } \keyword{algebra} \keyword{array} diff -Nru rmatrix-1.1-3/man/condest.Rd rmatrix-1.1-4/man/condest.Rd --- rmatrix-1.1-3/man/condest.Rd 2011-09-17 09:34:27.000000000 +0000 +++ rmatrix-1.1-4/man/condest.Rd 2014-04-12 21:35:47.000000000 +0000 @@ -5,7 +5,7 @@ \description{ \dQuote{Estimate}, i.e. compute approximately the CONDition number of a (potentially large, often sparse) matrix \code{A}. - It works by apply a fast approximation of the 1-norm, + It works by apply a fast \emph{randomized} approximation of the 1-norm, \code{norm(A,"1")}, through \code{onenormest(.)}. } \usage{ @@ -32,21 +32,30 @@ respectively.} \item{n}{\code{ == nrow(A)}, only needed when \code{A} is not specified.} \item{iter.max}{maximal number of iterations for the 1-norm estimator.} - \item{eps}{the relaive change that is deemed irrelevant.} + \item{eps}{the relative change that is deemed irrelevant.} } -% \details{ -% %% ~~ If necessary, more details than the description above ~~ -% } -\value{Both functions return a \code{\link{list}}; - \code{onenormest()} with components, - \item{est}{a number \eqn{> 0}, the estimated \code{norm(A, "1")}.} - \item{v}{the maximal \eqn{A X} column.} +\details{ + \code{\link{condest}()} calls \code{\link{lu}(A)}, and subsequently + \code{onenormest(A.x = , At.x = )} to compute an approximate norm of + the \emph{inverse} of \code{A}, \eqn{A^{-1}}, in a way which + keeps using sparse matrices efficiently when \code{A} is sparse. - The function \code{condest()} returns a list with components, - \item{est}{a number \eqn{> 0}, the estimated condition number + Note that \code{onenormest()} uses random vectors and hence + \emph{both} functions' results are random, i.e., depend on the random + seed, see, e.g., \code{\link{set.seed}()}. +} +\value{Both functions return a \code{\link{list}}; + \code{condest()} with components, + \item{est}{a number \eqn{> 0}, the estimated (1-norm) condition number \eqn{\hat\kappa}{k.}; when \eqn{r :=}\code{rcond(A)}, \eqn{1/\hat\kappa \approx r}{1/k. ~= r}.} - \item{v}{integer vector length \code{n}, with an \code{1} at the index + \item{v}{the maximal \eqn{A x} column, scaled to norm(v) = 1. + Consequently, \eqn{norm(A v) = norm(A) / est}; + when \code{est} is large, \code{v} is an approximate null vector.} + + The function \code{onenormest()} returns a list with components, + \item{est}{a number \eqn{> 0}, the estimated \code{norm(A, "1")}.} + \item{v}{0-1 integer vector length \code{n}, with an \code{1} at the index \code{j} with maximal column \code{A[,j]} in \eqn{A}.} \item{w}{numeric vector, the largest \eqn{A x} found.} \item{iter}{the number of iterations used.} @@ -76,9 +85,19 @@ data(KNex) mtm <- with(KNex, crossprod(mm)) system.time(ce <- condest(mtm)) +sum(abs(ce$v)) ## || v ||_1 == 1 +## Prove that || A v || = || A || / est (as ||v|| = 1): +stopifnot(all.equal(norm(mtm \%*\% ce$v), + norm(mtm) / ce$est)) + ## reciprocal 1 / ce$est system.time(rc <- rcond(mtm)) # takes ca 3 x longer rc all.equal(rc, 1/ce$est) # TRUE -- the approxmation was good + +one <- onenormest(mtm) +str(one) ## est = 12.3 +## the maximal column: +which(one$v == 1) # mostly 4, rarely 1, depending on random seed } diff -Nru rmatrix-1.1-3/man/dgeMatrix-class.Rd rmatrix-1.1-4/man/dgeMatrix-class.Rd --- rmatrix-1.1-3/man/dgeMatrix-class.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/dgeMatrix-class.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -92,16 +92,14 @@ \item{eigen}{\code{signature(x = "dgeMatrix", only.values= "missing")}: ...} \item{norm}{\code{signature(x = "dgeMatrix", type = "character")}: ... } \item{norm}{\code{signature(x = "dgeMatrix", type = "missing")}: ... } - \item{rcond}{\code{signature(x = "dgeMatrix", norm = "character")}: ... } - \item{rcond}{\code{signature(x = "dgeMatrix", norm = "missing")}: ... } + \item{rcond}{\code{signature(x = "dgeMatrix", norm = "character")} + or \code{norm = "missing"}: + the reciprocal condition number, \code{\link{rcond}()}.} \item{rowMeans}{\code{signature(x = "dgeMatrix")}: rowwise means (averages)} \item{rowSums}{\code{signature(x = "dgeMatrix")}: rowwise sums} - \item{t}{\code{signature(x = "dgeMatrix")}: ... } + \item{t}{\code{signature(x = "dgeMatrix")}: matrix transpose} } } -%\references{} -%\author{} -%\note{} \seealso{ Classes \code{\linkS4class{Matrix}}, \code{\linkS4class{dtrMatrix}}, and \code{\linkS4class{dsyMatrix}}. diff -Nru rmatrix-1.1-3/man/Diagonal.Rd rmatrix-1.1-4/man/Diagonal.Rd --- rmatrix-1.1-3/man/Diagonal.Rd 2012-10-05 22:18:37.000000000 +0000 +++ rmatrix-1.1-4/man/Diagonal.Rd 2014-04-12 22:20:38.000000000 +0000 @@ -25,13 +25,13 @@ from this argument, either \code{"U"} or \code{"L"}. Only rarely will it make sense to change this from the default.} \item{shape}{string of 1 character, one of \code{c("t","s","g")}, to - chose a triangular, symmetric or general result matrix.} + choose a triangular, symmetric or general result matrix.} \item{unitri}{optional logical indicating if a triangular result should be \dQuote{unit-triangular}, i.e., with \code{diag = "U"} slot, if possible. The default, \code{\link{missing}}, is the same as \code{\link{TRUE}}.} \item{kind}{string of 1 character, one of \code{c("d","l","n")}, to - chose the storage mode of the result, from classes + choose the storage mode of the result, from classes \code{\linkS4class{dsparseMatrix}}, \code{\linkS4class{lsparseMatrix}}, or \code{\linkS4class{nsparseMatrix}}, respectively.} diff -Nru rmatrix-1.1-3/man/dsyMatrix-class.Rd rmatrix-1.1-4/man/dsyMatrix-class.Rd --- rmatrix-1.1-3/man/dsyMatrix-class.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/dsyMatrix-class.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -75,10 +75,12 @@ \item{coerce}{\code{signature(from = "dsyMatrix", to = "matrix")}} \item{coerce}{\code{signature(from = "dsyMatrix", to = "dspMatrix")}} \item{coerce}{\code{signature(from = "dspMatrix", to = "dsyMatrix")}} - \item{norm}{\code{signature(x = "dspMatrix", type = "character")}} - \item{norm}{\code{signature(x = "dsyMatrix", type = "character")}} - \item{norm}{\code{signature(x = "dspMatrix", type = "missing")}} - \item{norm}{\code{signature(x = "dsyMatrix", type = "missing")}} + \item{norm}{\code{signature(x = "dspMatrix", type = "character")}, or + \code{x = "dsyMatrix"} or \code{type = "missing"}: Computes the + matrix norm of the desired type, see, \code{\link{norm}}.} + \item{rcond}{\code{signature(x = "dspMatrix", type = "character")}, or + \code{x = "dsyMatrix"} or \code{type = "missing"}: Computes the + reciprocal condition number, \code{\link{rcond}()}.} \item{solve}{\code{signature(a = "dspMatrix", b = "....")}, and} \item{solve}{\code{signature(a = "dsyMatrix", b = "....")}: \code{x <- solve(a,b)} solves \eqn{A x = b} for \eqn{x}; see @@ -90,7 +92,6 @@ } } %\references{} -%\author{} \seealso{ Classes \code{\linkS4class{dgeMatrix}} and \code{\linkS4class{Matrix}}; \code{\link[base]{solve}}, \code{\link{norm}}, \code{\link{rcond}}, diff -Nru rmatrix-1.1-3/man/image-methods.Rd rmatrix-1.1-4/man/image-methods.Rd --- rmatrix-1.1-3/man/image-methods.Rd 2013-12-23 14:09:05.000000000 +0000 +++ rmatrix-1.1-4/man/image-methods.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -64,9 +64,10 @@ the case of rastering.} \item{col.regions}{vector of gradually varying colors; see \code{\link[lattice]{levelplot}}.} - \item{lwd}{non-negative number or \code{NULL} (default), specifying the - line-width of the rectangles of each non-zero matrix entry (drawn by - \code{\link[grid]{grid.rect}}). The default depends on the matrix + \item{lwd}{(only used when \code{useRaster} is false:) non-negative + number or \code{NULL} (default), specifying the line-width of the + rectangles of each non-zero matrix entry (drawn by + \code{\link[grid]{grid.rect}}). The default depends on the matrix dimension and the device size.} \item{\dots}{further arguments passed to methods and \code{\link[lattice]{levelplot}}, notably \code{at} for specifying diff -Nru rmatrix-1.1-3/man/KhatriRao.Rd rmatrix-1.1-4/man/KhatriRao.Rd --- rmatrix-1.1-3/man/KhatriRao.Rd 2013-08-31 21:51:18.000000000 +0000 +++ rmatrix-1.1-4/man/KhatriRao.Rd 2014-04-26 19:14:30.000000000 +0000 @@ -78,4 +78,3 @@ } \keyword{methods} \keyword{array} - diff -Nru rmatrix-1.1-3/man/lu.Rd rmatrix-1.1-4/man/lu.Rd --- rmatrix-1.1-3/man/lu.Rd 2013-08-08 21:06:22.000000000 +0000 +++ rmatrix-1.1-4/man/lu.Rd 2014-04-26 19:14:30.000000000 +0000 @@ -4,6 +4,7 @@ \alias{lu,matrix-method} \alias{lu,dgeMatrix-method} \alias{lu,dgCMatrix-method} +\alias{lu,dtCMatrix-method} \usage{ lu(x, \dots) \S4method{lu}{matrix}(x, warnSing = TRUE, \dots) @@ -27,7 +28,7 @@ returns \code{\link{NA}} instead of an LU decomposition. No warning is signalled and the useR should be careful in that case. } - \item{order}{logical or integer, used to chose which will-reducing + \item{order}{logical or integer, used to choose which fill-reducing permutation technique will be used internally. Do not change unless you know what you are doing.} \item{tol}{positive number indicating the pivoting tolerance used in diff -Nru rmatrix-1.1-3/man/matrix-products.Rd rmatrix-1.1-4/man/matrix-products.Rd --- rmatrix-1.1-3/man/matrix-products.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/matrix-products.Rd 2014-05-08 10:31:52.000000000 +0000 @@ -7,24 +7,17 @@ %%-- %*% ------------------------------------- \alias{\%*\%,dgeMatrix,dgeMatrix-method} \alias{\%*\%,dgeMatrix,matrix-method} -\alias{\%*\%,dgeMatrix,numeric-method} +\alias{\%*\%,dgeMatrix,numLike-method} \alias{\%*\%,matrix,dgeMatrix-method} -\alias{\%*\%,numeric,dgeMatrix-method} +\alias{\%*\%,numLike,dgeMatrix-method} \alias{\%*\%,numLike,CsparseMatrix-method} \alias{\%*\%,CsparseMatrix,CsparseMatrix-method} \alias{\%*\%,CsparseMatrix,ddenseMatrix-method} \alias{\%*\%,CsparseMatrix,matrix-method} -\alias{\%*\%,CsparseMatrix,numeric-method} +\alias{\%*\%,CsparseMatrix,numLike-method} \alias{\%*\%,ddenseMatrix,CsparseMatrix-method} \alias{\%*\%,matrix,CsparseMatrix-method} -\alias{\%*\%,numeric,CsparseMatrix-method} \alias{\%*\%,ddenseMatrix,ddenseMatrix-method} -\alias{\%*\%,dgCMatrix,dgeMatrix-method} -\alias{\%*\%,dgCMatrix,matrix-method} -\alias{\%*\%,numeric,dgCMatrix-method} -\alias{\%*\%,ddenseMatrix,dgCMatrix-method} -\alias{\%*\%,dgeMatrix,dgCMatrix-method} -\alias{\%*\%,matrix,dgCMatrix-method} \alias{\%*\%,dgeMatrix,diagonalMatrix-method} \alias{\%*\%,matrix,diagonalMatrix-method} \alias{\%*\%,diagonalMatrix,dgeMatrix-method} @@ -104,23 +97,20 @@ \alias{crossprod-methods} \alias{crossprod,dgeMatrix,dgeMatrix-method} \alias{crossprod,dgeMatrix,matrix-method} -\alias{crossprod,dgeMatrix,numeric-method} +\alias{crossprod,dgeMatrix,numLike-method} \alias{crossprod,dgeMatrix,missing-method} \alias{crossprod,matrix,dgeMatrix-method} \alias{crossprod,numLike,dgeMatrix-method} -\alias{crossprod,numLike,dgeMatrix-method} -\alias{crossprod,numeric,CsparseMatrix-method} +\alias{crossprod,numLike,CsparseMatrix-method} \alias{crossprod,CsparseMatrix,missing-method} \alias{crossprod,CsparseMatrix,CsparseMatrix-method} \alias{crossprod,CsparseMatrix,ddenseMatrix-method} \alias{crossprod,CsparseMatrix,matrix-method} -\alias{crossprod,CsparseMatrix,numeric-method} +\alias{crossprod,CsparseMatrix,numLike-method} \alias{crossprod,ddenseMatrix,CsparseMatrix-method} \alias{crossprod,matrix,CsparseMatrix-method} \alias{crossprod,ddenseMatrix,missing-method} \alias{crossprod,ddenseMatrix,dgCMatrix-method} -\alias{crossprod,dgCMatrix,missing-method} -\alias{crossprod,dgCMatrix,matrix-method} \alias{crossprod,dgCMatrix,dgeMatrix-method} \alias{crossprod,CsparseMatrix,diagonalMatrix-method} \alias{crossprod,diagonalMatrix,CsparseMatrix-method} @@ -181,9 +171,9 @@ \alias{crossprod,Matrix,TsparseMatrix-method} \alias{crossprod,TsparseMatrix,TsparseMatrix-method} \alias{crossprod,TsparseMatrix,missing-method} -\alias{crossprod,dgTMatrix,missing-method} -\alias{crossprod,dgTMatrix,matrix-method} -\alias{crossprod,dgTMatrix,numeric-method} +\alias{crossprod,symmetricMatrix,Matrix-method} +\alias{crossprod,symmetricMatrix,missing-method} +\alias{crossprod,symmetricMatrix,ANY-method} %%-- tcrossprod ------------------------------ \alias{tcrossprod-methods} \alias{tcrossprod,dgeMatrix,missing-method} @@ -192,18 +182,12 @@ \alias{tcrossprod,dgeMatrix,numLike-method} \alias{tcrossprod,matrix,dgeMatrix-method} \alias{tcrossprod,numLike,dgeMatrix-method} -\alias{tcrossprod,dgCMatrix,missing-method} -\alias{tcrossprod,dgTMatrix,missing-method} \alias{tcrossprod,CsparseMatrix,ddenseMatrix-method} \alias{tcrossprod,CsparseMatrix,matrix-method} \alias{tcrossprod,CsparseMatrix,numLike-method} -\alias{tcrossprod,Matrix,numLike-method} \alias{tcrossprod,ddenseMatrix,CsparseMatrix-method} -\alias{tcrossprod,dgeMatrix,numLike-method} \alias{tcrossprod,matrix,CsparseMatrix-method} \alias{tcrossprod,numLike,CsparseMatrix-method} -\alias{tcrossprod,numLike,Matrix-method} -\alias{tcrossprod,numLike,dgeMatrix-method} \alias{tcrossprod,CsparseMatrix,CsparseMatrix-method} \alias{tcrossprod,CsparseMatrix,missing-method} \alias{tcrossprod,ddenseMatrix,missing-method} @@ -218,13 +202,13 @@ \alias{tcrossprod,matrix,diagonalMatrix-method} \alias{tcrossprod,sparseMatrix,diagonalMatrix-method} \alias{tcrossprod,ANY,Matrix-method} -\alias{tcrossprod,Matrix,numeric-method} +\alias{tcrossprod,numLike,Matrix-method} +\alias{tcrossprod,Matrix,numLike-method} \alias{tcrossprod,Matrix,ANY-method} \alias{tcrossprod,Matrix,missing-method} \alias{tcrossprod,Matrix,Matrix-method} \alias{tcrossprod,Matrix,matrix-method} \alias{tcrossprod,matrix,Matrix-method} -\alias{tcrossprod,numeric,Matrix-method} \alias{tcrossprod,ddenseMatrix,dtrMatrix-method} \alias{tcrossprod,dtrMatrix,dtrMatrix-method} \alias{tcrossprod,matrix,dtrMatrix-method} @@ -254,6 +238,15 @@ \alias{tcrossprod,Matrix,TsparseMatrix-method} \alias{tcrossprod,TsparseMatrix,TsparseMatrix-method} \alias{tcrossprod,TsparseMatrix,missing-method} +\alias{tcrossprod,Matrix,symmetricMatrix-method} +\alias{tcrossprod,ANY,symmetricMatrix-method} +% (special class methods, added only for speedup:) +\alias{tcrossprod,ddenseMatrix,dsCMatrix-method} +\alias{tcrossprod,ddenseMatrix,lsCMatrix-method} +\alias{tcrossprod,ddenseMatrix,nsCMatrix-method} +\alias{tcrossprod,matrix,dsCMatrix-method} +\alias{tcrossprod,matrix,lsCMatrix-method} +\alias{tcrossprod,matrix,nsCMatrix-method} %%--------------- %- these three are needed for R CMD check: \alias{\%*\%} @@ -320,9 +313,9 @@ ditto for several other signatures, use \code{showMethods("crossprod", class = "dgeMatrix")}, matrix crossproduct, an efficient version of \code{t(x) \%*\% y}.} - \item{crossprod}{\code{signature(x = "dgCMatrix", y = "missing")} + \item{crossprod}{\code{signature(x = "CsparseMatrix", y = "missing")} returns \code{t(x) \%*\% x} as an \code{dsCMatrix} object.} - \item{crossprod}{\code{signature(x = "dgTMatrix", y = "missing")} + \item{crossprod}{\code{signature(x = "TsparseMatrix", y = "missing")} returns \code{t(x) \%*\% x} as an \code{dsCMatrix} object.} \item{crossprod,tcrossprod}{\code{signature(x = "dtrMatrix", y = diff -Nru rmatrix-1.1-3/man/norm.Rd rmatrix-1.1-4/man/norm.Rd --- rmatrix-1.1-3/man/norm.Rd 2010-05-01 20:49:55.000000000 +0000 +++ rmatrix-1.1-4/man/norm.Rd 2014-04-12 21:35:47.000000000 +0000 @@ -43,6 +43,10 @@ \code{dlange}, \code{dlansy}, \code{dlantr}, \code{zlange}, \code{zlansy}, and \code{zlantr}. } +\seealso{ + \code{\link{onenormest}()}, an \emph{approximate} randomized estimate + of the 1-norm condition number, efficient for large sparse matrices. +} \references{ Anderson, E., et al. (1994). \emph{LAPACK User's Guide,} diff -Nru rmatrix-1.1-3/man/pMatrix-class.Rd rmatrix-1.1-4/man/pMatrix-class.Rd --- rmatrix-1.1-3/man/pMatrix-class.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/pMatrix-class.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -89,12 +89,6 @@ or column permutation. Here are the four different cases for an arbitrary matrix \eqn{M} and a permutation matrix \eqn{P} (where we assume matching dimensions): - %% \describe{ - %% \item{\eqn{MP = }}{\code{M \%*\% P} ..} - %% \item{\eqn{PM = }}{\code{P \%*\% M} ..} - %% \item{\eqn{P'M = }}{\code{crossprod(P,M)} (less efficiently, \code{t(P) \%*\% M}) ..} - %% \item{\eqn{MP' = }}{\code{trossprod(M,P)} (less efficiently, \code{M \%*\% t(P)}) ..} - %% } \tabular{lclcl}{ \eqn{MP }\tab= \tab\code{M \%*\% P} \tab= \tab\code{M[, i(p)]}\cr \eqn{PM }\tab= \tab\code{P \%*\% M} \tab= \tab\code{M[ p , ]} \cr diff -Nru rmatrix-1.1-3/man/rcond.Rd rmatrix-1.1-4/man/rcond.Rd --- rmatrix-1.1-3/man/rcond.Rd 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/rcond.Rd 2014-04-12 21:35:47.000000000 +0000 @@ -1,9 +1,8 @@ \name{rcond} \title{Estimate the Reciprocal Condition Number} \alias{rcond} -% most methods are documented in Matrix-class.Rd +% most methods are "documented" in Matrix-class.Rd \alias{rcond,ANY,missing-method} -\alias{rcond,matrix,character-method}%<- no longer for R >= 2.7.0 \alias{rcond,Matrix,character-method} \alias{rcond,ldenseMatrix,character-method} \alias{rcond,ndenseMatrix,character-method} @@ -22,11 +21,12 @@ } \arguments{ \item{x}{an \R object that inherits from the \code{Matrix} class.} - \item{norm}{ - Character indicating the type of norm to be used in the estimate. - The default is \code{"O"} for the 1-norm (\code{"O"} is equivalent - to \code{"1"}). The other possible value is \code{"I"} for the - infinity norm, see also \code{\link{norm}}. + \item{norm}{character string indicating the type of norm to be used in + the estimate. The default is \code{"O"} for the 1-norm (\code{"O"} is + equivalent to \code{"1"}). For sparse matrices, when \code{useInv=TRUE}, + \code{norm} can be any of the \code{kind}s allowed for \code{\link{norm}}; + otherwise, the other possible value is \code{"I"} for the infinity + norm, see also \code{\link{norm}}. } \item{useInv}{logical (or \code{"Matrix"} containing \code{\link{solve}(x)}). If not false, compute the reciprocal @@ -82,7 +82,7 @@ \eqn{p=2} (Euclidean) \code{\link{norm}}. \code{\link[base]{solve}}. - \code{\link{condest}}, a newer approximate estimate of + \code{\link{condest}}, a newer \emph{approximate} estimate of the (1-norm) condition number, particularly efficient for large sparse matrices. } @@ -117,12 +117,15 @@ 36*(iM <- solve(M)) # still sparse MM <- kronecker(Diagonal(10), kronecker(Diagonal(5),kronecker(m,M))) dim(M3 <- kronecker(bdiag(M,M),MM)) # 12'800 ^ 2 -if(interactive()) ## takes about 2 seconds +if(interactive()) ## takes about 2 seconds if you have >= 8 GB RAM system.time(r <- rcond(M3)) ## whereas this is *fast* even though it computes solve(M3) system.time(r. <- rcond(M3, useInv=TRUE)) if(interactive()) ## the values are not the same c(r, r.) # 0.05555 0.013888 +## for all 4 norms available for sparseMatrix : +cbind(rr <- sapply(c("1","I","F","M"), + function(N) rcond(M3, norm=N, useInv=TRUE))) \dontshow{stopifnot(all.equal(r., 1/72, tolerance=1e-12))} } \keyword{array} diff -Nru rmatrix-1.1-3/man/rsparsematrix.Rd rmatrix-1.1-4/man/rsparsematrix.Rd --- rmatrix-1.1-3/man/rsparsematrix.Rd 1970-01-01 00:00:00.000000000 +0000 +++ rmatrix-1.1-4/man/rsparsematrix.Rd 2014-04-26 19:45:28.000000000 +0000 @@ -0,0 +1,65 @@ +\name{rsparsematrix} +\alias{rsparsematrix} +\title{Random Sparse Matrix} +\description{ + Generate a Random Sparse Matrix Efficiently. +} +\usage{ +rsparsematrix(nrow, ncol, density, nnz = round(density * maxE), + symmetric = FALSE, + rand.x = function(n) signif(rnorm(nnz), 2), \dots) +} +\arguments{ + \item{nrow, ncol}{number of rows and columns, i.e., the matrix + dimension (\code{\link{dim}}).} + \item{density}{optional number in \eqn{[0,1]}, the density is the + proportion of non-zero entries among all matrix entries. If + specified it determines the default for \code{nnz}, otherwise + \code{nnz} needs to be specified.} + \item{nnz}{number of non-zero entries, for a sparse matrix typically + considerably smaller than \code{nrow*ncol}. Must be specified if + \code{density} is not.} + \item{symmetric}{logical indicating if result should be a matrix of + class \code{\linkS4class{symmetricMatrix}}. Note that in the symmetric + case, \code{nnz} denotes the number of non zero entries of the upper + (or lower) part of the matrix, including the diagonal.} + \item{rand.x}{the random number generator for the \code{x} slot, a + \code{\link{function}} such that \code{rand.x(n)} generates a + numeric vector of length \code{n}. Typical examples are + \code{rand.x = rnorm}, or \code{rand.x = runif}; the default is nice + for didactical purposes.} + \item{\dots}{optionally further arguments passed to + \code{\link{sparseMatrix}()}, notably \code{giveCsparse}.} +} +\details{ + The algorithm first samples \dQuote{encoded} \eqn{(i,j)}s without + replacement, via one dimensional indices, if not \code{symmetric} + \code{\link{sample.int}(nrow*ncol, nnz)}, then gets + \code{x <- rand.x(nnz)} and calls + \code{\link{sparseMatrix}(i=i, j=j, x=x, ..)}. +} +\value{ + a \code{\linkS4class{sparseMatrix}}, say \code{M} of dimension (nrow, + ncol), i.e., with \code{dim(M) == c(nrow, ncol)}, if \code{symmetric} + is not true, with \code{nzM <- \link{nnzero}(M)} fulfilling + \code{nzM <= nnz} and typically, \code{nzM == nnz}. +} +\author{Martin Maechler} +\examples{ +set.seed(17)# to be reproducible +M <- rsparsematrix(8, 12, nnz = 30) # small example, not very sparse +M +M1 <- rsparsematrix(1000, 20, nnz = 123, rand.x = runif) +summary(M1) + +## a random *symmetric* Matrix +(S9 <- rsparsematrix(9, 9, nnz = 10, symmetric=TRUE)) # dsCMatrix +nnzero(S9)# ~ 20: as 'nnz' only counts one "triangle" + +## a [T]riplet representation sparseMatrix: +T2 <- rsparsematrix(40, 12, nnz = 99, giveCsparse=FALSE) +head(T2) +} +\keyword{array} +\keyword{distribution} + diff -Nru rmatrix-1.1-3/man/sparseMatrix.Rd rmatrix-1.1-4/man/sparseMatrix.Rd --- rmatrix-1.1-3/man/sparseMatrix.Rd 2014-03-25 12:41:40.000000000 +0000 +++ rmatrix-1.1-4/man/sparseMatrix.Rd 2014-04-26 19:14:30.000000000 +0000 @@ -6,7 +6,7 @@ matrix, inheriting from \code{\link{class}} \code{\linkS4class{CsparseMatrix}} (or \code{\linkS4class{TsparseMatrix}} if \code{giveCsparse} is false), - from locations (and values) of its nonzero entries. + from locations (and values) of its non-zero entries. This is the recommended user interface rather than direct \code{\link{new}("***Matrix", ....)} calls. @@ -90,7 +90,7 @@ \code{"\linkS4class{TsparseMatrix}"} class, unless \code{use.last.ij} is set to true. %% - By default, when \code{giveCsparse} is true, the \code{\linkS4class{CsparseMatrix}} + By default, when \code{giveCsparse} is true, the \code{\linkS4class{CsparseMatrix}} derived from this triplet form is returned. The reason for returning a \code{\linkS4class{CsparseMatrix}} object @@ -135,9 +135,9 @@ (A1 <- sparseMatrix(i[perm], j[perm], x = x[perm])) stopifnot(identical(A, A1)) -## The slots are 0-index based, so +## The slots are 0-index based, so try( sparseMatrix(i=A@i, p=A@p, x= seq_along(A@x)) ) -## fails and you should say so: 1-indexing is FALSE: +## fails and you should say so: 1-indexing is FALSE: sparseMatrix(i=A@i, p=A@p, x= seq_along(A@x), index1 = FALSE) ## the (i,j) pairs can be repeated, in which case the x's are summed @@ -214,4 +214,3 @@ }% donttest }% example \keyword{array} - diff -Nru rmatrix-1.1-3/man/sparseVector.Rd rmatrix-1.1-4/man/sparseVector.Rd --- rmatrix-1.1-3/man/sparseVector.Rd 2012-09-03 07:34:21.000000000 +0000 +++ rmatrix-1.1-4/man/sparseVector.Rd 2014-04-26 19:14:30.000000000 +0000 @@ -5,7 +5,7 @@ User friendly construction sparse vectors, i.e., objects inheriting from \code{\link{class}} \code{\linkS4class{sparseVector}}, from indices and values of its - nonzero entries. + non-zero entries. } \details{ zero entries in \code{x} are dropped automatically, analogously as diff -Nru rmatrix-1.1-3/man/USCounties.Rd rmatrix-1.1-4/man/USCounties.Rd --- rmatrix-1.1-3/man/USCounties.Rd 2013-07-03 10:13:08.000000000 +0000 +++ rmatrix-1.1-4/man/USCounties.Rd 2014-04-12 21:37:37.000000000 +0000 @@ -41,21 +41,23 @@ system.time(MJ <- sapply(rho, function(x) determinant(IM - x * USCounties, logarithm = TRUE)$modulus)) +## can be done faster, by update()ing the Cholesky factor: nWC <- -USCounties C1 <- Cholesky(nWC, Imult = 2) system.time(MJ1 <- n * log(rho) + sapply(rho, function(x) - c(determinant(update(C1, nWC, 1/x))$modulus))) + 2 * c(determinant(update(C1, nWC, 1/x))$modulus))) all.equal(MJ, MJ1) +\dontshow{stopifnot( all.equal(MJ, MJ1) )} C2 <- Cholesky(nWC, super = TRUE, Imult = 2) system.time(MJ2 <- n * log(rho) + sapply(rho, function(x) - c(determinant(update(C2, nWC, 1/x))$modulus))) -all.equal(MJ, MJ2) + 2 * c(determinant(update(C2, nWC, 1/x))$modulus))) +all.equal(MJ, MJ2) \dontshow{stopifnot(all.equal(MJ, MJ2))} system.time(MJ3 <- n * log(rho) + Matrix:::ldetL2up(C1, nWC, 1/rho)) -all.equal(MJ, MJ3) +stopifnot(all.equal(MJ, MJ3)) system.time(MJ4 <- n * log(rho) + Matrix:::ldetL2up(C2, nWC, 1/rho)) -all.equal(MJ, MJ4) +stopifnot(all.equal(MJ, MJ4)) } \keyword{datasets} diff -Nru rmatrix-1.1-3/MD5 rmatrix-1.1-4/MD5 --- rmatrix-1.1-3/MD5 2014-03-31 10:51:11.000000000 +0000 +++ rmatrix-1.1-4/MD5 2014-06-15 18:26:41.000000000 +0000 @@ -1,21 +1,21 @@ cbb51e0c5dddc4dfa7f38664828ec67f *ChangeLog -5c949a4bb7238e81556209482d0dc148 *DESCRIPTION -cffe04628af18a46857699979fd93afc *NAMESPACE -76221218654d33eb611ba3d717a2eeaf *R/AllClass.R +cf2358020be83f6b0b8bfed5fb30b82b *DESCRIPTION +354815bf59af7396a92af141be3081fe *NAMESPACE +fb15af8ffff4f48935acd88cacc0f2bb *R/AllClass.R 17f3ee5af7800f04cbdd44ae706720db *R/AllGeneric.R -1e84512b1b5c1a943399734c69342516 *R/Auxiliaries.R -83a9a959b3b605b931840aeed827ec0b *R/CHMfactor.R -bd1eaca6c99d42e640e8ad3b9dcb2700 *R/Csparse.R +41f415257f9e7f45d36ca02289b24772 *R/Auxiliaries.R +7025e2165b918d54fc398d70f41035c2 *R/CHMfactor.R +3758c175b4a44cc89c98dc8dc034ff8f *R/Csparse.R c900a4080317d3a6b766324be4e64189 *R/HBMM.R 3f37d43942f54523314df6f5406cf488 *R/Hilbert.R 3eacab76e65bc0493686a2dab2646e78 *R/KhatriRao.R 7cd5d5938844fa73c1a2b50fae7fdd19 *R/LU.R -a89714a5731f72637323113c6d304cca *R/Matrix.R +510c7aceb9e4f7be0d8a2ed1718a8139 *R/Matrix.R e19ed36934229e22f41552d4d40c5d5a *R/MatrixFactorization.R -cc57f1182b998bf25a87411b631ae8ce *R/Ops.R +4276c67690f81171fa178ee214bc5c1f *R/Ops.R 5f24daa36457fcc57724cd57adebfca6 *R/Rsparse.R 6433ca318f74c40760e68b1d56466d09 *R/SparseM-conv.R -8180622d77731ef6419b057ff09ad5f5 *R/Tsparse.R +68b6185aede9b53929b190695040cc75 *R/Tsparse.R 9cf8de63fa0b6f64df31cd32d89eb12a *R/abIndex.R 4b467ca6296fe89da23404d61644cd2d *R/bandSparse.R 9df38d1b4f2dff615379b62e72220564 *R/bind2.R @@ -25,24 +25,24 @@ 2ccd198e6a1bd31debbcc7d07ae320f4 *R/dMatrix.R 1f35e9187df902b674a21d02c1bb7228 *R/ddenseMatrix.R 8b7796b05cfd393b2fb069cf9e4852fe *R/denseMatrix.R -7b809fd532b95dee1640290602357b9f *R/dgCMatrix.R -7049a824e11e25381aabf2803213ad3f *R/dgTMatrix.R +0cec552e2e2715fc1c718caa205193a1 *R/dgCMatrix.R +1db86f9b20fbd9c15386b17134eb30ef *R/dgTMatrix.R 6c60e41bcea1653b0310a4c70947c1a8 *R/dgeMatrix.R -0d4d855785a372e711af1735e60ff007 *R/diagMatrix.R -cc13b7e9cc4a37c38253ed868f02d862 *R/dpoMatrix.R -b9cedc3fdadf705cdee22603d3b4e89f *R/dppMatrix.R -f0bd64abec232f7b0e38fadea05aa1c2 *R/dsCMatrix.R +b15d610f1132bad18181a833fca8a947 *R/diagMatrix.R +be76b734bbed865c5a76f90ecd874448 *R/dpoMatrix.R +239dd04632141e119924f4b3a15e2188 *R/dppMatrix.R +90547d24e15ddb234093b1c4254f6e05 *R/dsCMatrix.R 3ee77774e2a81752173065ef0f4100db *R/dsTMatrix.R 6914bb8994c5854ba34deb23bc1648ea *R/dspMatrix.R 4049e470e739d985bfb0b4422f474049 *R/dsparseMatrix.R 65ad0d532603c978905cf71ab9738ebe *R/dsyMatrix.R -304128d12d4644b5238efc73d917085e *R/dtCMatrix.R +98a8cbf502792ff5b19aefce4de9bfaa *R/dtCMatrix.R e826da941555bb348d828af69c77b7dc *R/dtTMatrix.R 784978f7c2816709d791e5b03f770b27 *R/dtpMatrix.R -65434e63c4bd207eeae394a284687466 *R/dtrMatrix.R +fbc6a432ecd9457f2994eefd1701838a *R/dtrMatrix.R 952b2122422438fe96f3e8f08056b99a *R/eigen.R 333accc94fd2c5ba991ce2473e4cd510 *R/expm.R -c1307d2637d83c8ed1d8703601be64ad *R/indMatrix.R +b803b46feaca3381d6dd7f72df561ba0 *R/indMatrix.R 233e9458cb509610a5d8efde8cc0e4ce *R/kronecker.R 1435034568896038f41377240b744929 *R/lMatrix.R d7af31c7887001750716a0c09d3d7469 *R/ldenseMatrix.R @@ -65,16 +65,16 @@ c90a8c1c8f6328c3d257ce73ee8622dd *R/ntCMatrix.R f12c8b3ac756a91e8c0dfdf270a38b98 *R/ntTMatrix.R b8166a9d06aab03a53432358a59ad9e5 *R/pMatrix.R -a1d1122776dca542c9df6d3848695073 *R/products.R +7a7c2aa08e21b247fe60d09e2aafdfb7 *R/products.R 011ff226817062c7d05be602b2043b0c *R/rankMatrix.R 8dc90941d51a42e4270cd2ae717f3cb5 *R/spModels.R -d717d01c1e689833ec1a033e8a341681 *R/sparseMatrix.R +e22a414110bd441be51a8d3450ddd8da *R/sparseMatrix.R 9091e186403012bf94f7db9c01157ad6 *R/sparseQR.R 2cf18ea0dcddf37822b97909345ba14f *R/sparseVector.R 382304430f9b8df603f5b43e072c57c9 *R/symmetricMatrix.R de53e77d9edbc9130fc4b5c5d51996af *R/triangularMatrix.R 9b6dbb2197d9c745528e5964991c266c *R/zzz.R -7166f37c5cea86be8cb0b84ea7ed608d *TODO +b168a14dbc322902d9c72168868950b6 *TODO f2ad5375e270deeb7b041272dd095032 *cleanup d3528bc0cd82f248821809f253bffe35 *data/CAex.R 7586556183d8d0a3c7979702254f93c1 *data/KNex.R @@ -82,17 +82,17 @@ 8734f0b040c6292983d273d4251d250a *inst/Copyrights cd003f462d78d4471ff1e295f79f7b9a *inst/Doxyfile ea747a35318c564002a7144ebd510e28 *inst/doc/Announce.txt -005948925e8ce13edb44d9a5aabcbd47 *inst/doc/Comparisons.pdf -03e04274c4a70ea3478fc211fbab1bd5 *inst/doc/Design-issues.pdf -8d3be28a72c9bff8c7dfddab52fdf3aa *inst/doc/Intro2Matrix.pdf -f2e1a391a5b2ce4ba1b7fc153cca6d83 *inst/doc/Introduction.pdf +8db52169cd6161b437c7ecd604731747 *inst/doc/Comparisons.pdf +0b2a0d01a167ee2c985e1756c02ec15f *inst/doc/Design-issues.pdf +e6b15a992a718a0296f9e099e8974bfe *inst/doc/Intro2Matrix.pdf +6e0383e34836b3a2f3b538640f471a39 *inst/doc/Introduction.pdf 7d60c8defc517d361513bac03d04f13b *inst/doc/SuiteSparse/AMD.txt e4f8cd28fc8be8ab14106630610b9c5f *inst/doc/SuiteSparse/CHOLMOD.txt 0a80fb0db9cefb02e290ffd2f3431e6c *inst/doc/SuiteSparse/COLAMD.txt d75882d4bb768ba0ea352291652daaee *inst/doc/SuiteSparse/SPQR.txt 29429afaf74eaac98cb1957ddd4e8a38 *inst/doc/SuiteSparse/SuiteSparse_config.txt 6d217288f5da4fe419afaa34988bf42d *inst/doc/SuiteSparse/UserGuides.txt -9c520ea62a2485ffe96821e8a634fe96 *inst/doc/sparseModels.pdf +eb37f8e496657e01eda1a58b3c2a1269 *inst/doc/sparseModels.pdf dcd11f6947f910f743254824e930b2c7 *inst/external/CAex_slots.rda be886d6bb832210bb654b9ad064fe0ff *inst/external/KNex_slots.rda 90f019ec81e67d7f3542f7ca39bf3f2d *inst/external/USCounties_slots.rda @@ -115,18 +115,18 @@ 027e86d29cc531a651ad0ad45bf2b92d *inst/po/pl/LC_MESSAGES/Matrix.mo af8b2f8436f9bc41ecca82c1b8612b6c *inst/po/pl/LC_MESSAGES/R-Matrix.mo 977cf42988f0a24c9db3e4c902f2eced *inst/test-tools-1.R -668e5c5f3682ace3dc4690965f830508 *inst/test-tools-Matrix.R +2296a0ebafd3f94edfd436d0d3b67104 *inst/test-tools-Matrix.R 384ad19fd9e35dfa867ff36dc6d502f4 *inst/test-tools.R 008b2c14ce99ff19c7c5fa9e8c8736b1 *man/BunchKaufman-methods.Rd 3f382d978f8947e4e17047f95d20ed84 *man/CAex.Rd -80e9d96f95b44539ad8d9f90c4a88aad *man/CHMfactor-class.Rd +5b10b3332f494d1fc81c004541234d8b *man/CHMfactor-class.Rd f39cd7f9daaa240217e1d6ce0f810b47 *man/Cholesky-class.Rd -b02f4b300be2cb5c13a66f73ac221551 *man/Cholesky.Rd +5d5811453fa47a25fef4ac4d973385e3 *man/Cholesky.Rd 4ef1b053e942af5a8a2289f4e252d323 *man/CsparseMatrix-class.Rd -4e990802e5691f02b99bd439682626cc *man/Diagonal.Rd +808e4ca0630256acf52598f4a34d77d7 *man/Diagonal.Rd 74df4ccf77c8797a57520d36afe98a0f *man/Hilbert.Rd 2ad7cc5b8c267f23bf7012c5248f87a5 *man/KNex.Rd -8bff8b5e9ce36e848a959c76133ae4a1 *man/KhatriRao.Rd +ff1568bd735efff9a9c2c463c35c7a4a *man/KhatriRao.Rd b79c16e807bb0dc031ddb41644d605f6 *man/LU-class.Rd 464bc4f6abcbbb2d0caa36d151edac10 *man/Matrix-class.Rd 0241d3d91abca6086cb03e1292cb4432 *man/Matrix.Rd @@ -137,7 +137,7 @@ 601a206196fb38ff2c3eb6a696aebeef *man/SparseM-conv.Rd b946bf426edf45858278e6eaacb14f4a *man/Subassign-methods.Rd 3a229d88e2e56d2d8061fe509536e12a *man/TsparseMatrix-class.Rd -a4cba35d3e8296c63c2cd85fedef9c3a *man/USCounties.Rd +9e48553229ec7a335544a30511e0344c *man/USCounties.Rd 0502c8aba3321a819cd00dd1e553ab5d *man/Xtrct-methods.Rd 7abbeae6546c97bd4a59cabf15abe143 *man/abIndex-class.Rd 0611d202a92b93597a926a452ac5f525 *man/abIseq.Rd @@ -148,11 +148,11 @@ 8c5d951ae5615c8f93e5894f7e6b0581 *man/bandSparse.Rd e01f3ae7320eea9a3d51fa6776b8daa5 *man/bdiag.Rd a6cad9af0b78e935e26195af7a314190 *man/cBind.Rd -34de6b2dc0fcef0f5e1e71f371e19a1b *man/chol.Rd +8c7312211f0224f12b21c3342de442c7 *man/chol.Rd 56c41276d74b7373b67fe98947b62834 *man/chol2inv-methods.Rd c8f1619647e27a9f350230cf22ba4832 *man/colSums.Rd 30e5e49902a94ad96a9b343f1805a14e *man/compMatrix-class.Rd -2a7107c0cd0a6933ba15384f3557f018 *man/condest.Rd +bc004ec546340928a245be773286408d *man/condest.Rd 39cc7e7c6b303f13a79f9ea8ab50ba52 *man/dMatrix-class.Rd 67fcf5d02576b9c5a2d0b54237c0ec26 *man/ddenseMatrix-class.Rd fea7d6aa90395d29ab34ebe90603aec4 *man/ddiMatrix-class.Rd @@ -160,7 +160,7 @@ feec03c7616aabc27d3f7b5f0958e393 *man/dgCMatrix-class.Rd 33311e5379364ad7d6fa816ce187506d *man/dgRMatrix-class.Rd bb72c1b1e9f6837dc354b3c6aea3e14d *man/dgTMatrix-class.Rd -e139cd3b9a34a66157852af42e325c97 *man/dgeMatrix-class.Rd +2bbda14fe19c769958294b8006c83448 *man/dgeMatrix-class.Rd 7b339b007d683c02526d2ead6043bf31 *man/diagU2N.Rd 795d0536fd8de31cd75c30babad85c22 *man/diagonalMatrix-class.Rd cab5f0fa9bd7b9c584d9025ace66e4e6 *man/dpoMatrix-class.Rd @@ -168,7 +168,7 @@ 0c14415ce7caa0d8b5b8c2c03d004fbd *man/dsCMatrix-class.Rd e6232ded0c7b441674776381b8e25efa *man/dsRMatrix-class.Rd 5a84601e5f8ed336d3cc947e6703d271 *man/dsparseMatrix-class.Rd -df9c058ff1c06ba43dc972c7bdc30033 *man/dsyMatrix-class.Rd +1f52e18e5bc70601f1c49abb4ca5658e *man/dsyMatrix-class.Rd e264e1fe36ffd4b34d2290d9af6f37ae *man/dtCMatrix-class.Rd 9e88fe051b42ffdfc4fdc08361372a21 *man/dtRMatrix-class-def.Rd 3f894c8473e90abe89ca43c9d8cbe465 *man/dtpMatrix-class.Rd @@ -181,7 +181,7 @@ 829bc7814aa9d26fd79c0f6c7a04aaa0 *man/formatSparseM.Rd 10534177267141a74ee6f1403862bf3c *man/generalMatrix-class.Rd 0bbdeba42ac962a61cdb9610897c3e6f *man/graph2T.Rd -9c9398fda6fc2709cf3ae78ec101809d *man/image-methods.Rd +a766a9c8090bd54484c4a0a7f5dacb1c *man/image-methods.Rd 49d38ed3606eb5247ffaf67a28d66051 *man/indMatrix-class.Rd 828999cf8aaac364bf35e98b84993b62 *man/index-class.Rd 469733dd46acbb038a6fee7165b96d40 *man/invPerm.Rd @@ -196,35 +196,36 @@ 7a21658e580247b8c3f5f8b4ca1eaad8 *man/lsparseMatrix-classes.Rd f1e1cfefe3650a6c19b1c63b7d6052dc *man/lsyMatrix-class.Rd fc78ae1f310100dc8dbe8cff430a2093 *man/ltrMatrix-class.Rd -026c79df06995b80995e8d81828a33e0 *man/lu.Rd -6c41a352760ce9142fe8dd184584ecde *man/matrix-products.Rd +9a889abad6ddae7c62ca17e920c1391c *man/lu.Rd +6f2cc7c6fd26515085dc52d3a2f8a64f *man/matrix-products.Rd 0b78b90513390c5816114d3058f0777c *man/nMatrix-class.Rd 2d86741969a8e728299f64d762dc398a *man/ndenseMatrix-class.Rd 46dc1c9089fefb669bf1b8137edd2d12 *man/nearPD.Rd 84f3a5d496768cfe6f139e3fab1ef042 *man/ngeMatrix-class.Rd b0eb88e96c037d58bdee94dbc2a92d6b *man/nnzero.Rd -603f32761ee0519acaa2ab808a2b0f6b *man/norm.Rd +baf1f784c556584df45b0b3ff7152b18 *man/norm.Rd 3a14b23ea37d2d9e9ad6e8b62ead68fb *man/nsparseMatrix-classes.Rd b299d889bba52b4d4e558b7d121eb97a *man/nsyMatrix-class.Rd 11890bb8c3bdd03d0d0ca39ab95a6c08 *man/ntrMatrix-class.Rd b44976fc24c22cc0d7d0bf722a951a3b *man/number-class.Rd -8135acf1f7666ec879d5cb547e49df1a *man/pMatrix-class.Rd +358232f60f227a3908652716415f6b5c *man/pMatrix-class.Rd 6e9825c54dc87baf423237b3ff934f5b *man/printSpMatrix.Rd f61b7ca27c49d694c6b4b7b1e71d753b *man/qr-methods.Rd 89d4ef7a94e5eec8c8ae7a3e23218e99 *man/rankMatrix.Rd -4bd5ce50393913c61b06b55833d33dd2 *man/rcond.Rd +e6a10cf788588d6a28eb8b8ed328a790 *man/rcond.Rd eaa76f30797997d7d49d45c2f9c55c4f *man/rep2abI.Rd b6fedde6a4c71d315554de6eb2eba367 *man/replValue-class.Rd 12d2c9caa3c91beecd1f6c60903478f6 *man/rleDiff-class.Rd +110d2e592b5c2cea7aeed0ef69f02665 *man/rsparsematrix.Rd 0cc917f4625559f00b18b296608abe7a *man/solve-methods.Rd a1bf541febf7b913a88fea7dd9c597bf *man/spMatrix.Rd 47b173d7c7515d8616daab061f58df12 *man/sparse.model.matrix.Rd f5ebe1e0fff8a192a4130b16dfbdd6ee *man/sparseLU-class.Rd 9686d870370daac8f885bcda3ac027f7 *man/sparseMatrix-class.Rd -370b2589b760e512d48767e5594a15ad *man/sparseMatrix.Rd +716e1e656de5517cc960c190ccc588e6 *man/sparseMatrix.Rd df61fe96a74790c4fdde5e7997db2dbf *man/sparseQR-class.Rd e64c2259f1b8a9d65bc9ca986417a4bd *man/sparseVector-class.Rd -f3553d085b3b3d43c601c7b6904e1038 *man/sparseVector.Rd +e7bf1f347f18c79f54acb50f237c9ee0 *man/sparseVector.Rd 401173e5302fdc743b1a4478e80a1121 *man/symmetricMatrix-class.Rd 516edf3833b0d2f0d483dbb06c16800d *man/symmpart.Rd 94ceb2403c2ddf02133bc4aa9de9b3ca *man/triangularMatrix-class.Rd @@ -359,11 +360,11 @@ a40eda898c98a8b193aae541a40e16b8 *src/COLAMD/Source/Makefile ed584482056a57c97c84f0f2701f76a6 *src/COLAMD/Source/colamd.c 8f9d70daa9adafe39e555573d974d17e *src/COLAMD/Source/colamd_global.c -39aebdd9d2fd790e98d257c0864c5ddd *src/Csparse.c -2cf6a37fa7139ed3dde40b022f037d3b *src/Csparse.h +34babed4fed53dd234450c4ee9ceedad *src/Csparse.c +11927e2daadce9f293e11269c6f8f4e0 *src/Csparse.h 8f61a05750779ce84f2ed402f5596488 *src/Makevars -a30a54f8115095d49d435854ca6b2fef *src/Mutils.c -44b953d1c59f2762fb9710d2e795ad52 *src/Mutils.h +7d5fcc8e430cbc70bc439f34fbb1ef74 *src/Mutils.c +616bc4038562ee33f14f9f4b2e031e98 *src/Mutils.h 1579a754fcbf8ebabf9b30c31289fdfe *src/SuiteSparse_config/Makefile e9b2be5ad0056f588187cf504cf60764 *src/SuiteSparse_config/SuiteSparse_config.c 5be5db8ccb99f88f234639b97d8227c5 *src/SuiteSparse_config/SuiteSparse_config.h @@ -387,17 +388,17 @@ 65b0e724e0bb67bcd82b5f3f57ea2b58 *src/dgCMatrix.h 130569749628885030c543d404f7995d *src/dgTMatrix.c b79dce2247ffabcb51997b8fe591b41b *src/dgTMatrix.h -dd482ee965db53ef71376da44f42bc23 *src/dgeMatrix.c +ccd85dbef0fd48e8b7fc259143e4cd40 *src/dgeMatrix.c abb3c73f4a912de19114010daa9608ab *src/dgeMatrix.h 41d0e269a44f8ea536a97c5fea0aa8c1 *src/dpoMatrix.c b58f1f64d6141e98de9eaf6d953e7f6d *src/dpoMatrix.h cf1e54516cd78c8d4093b5ffd500ba27 *src/dppMatrix.c 4cb570fdeeab6cd03ba14db58efeddce *src/dppMatrix.h -745e074d389834cdcceb968b6ef9190e *src/dsCMatrix.c -55fb2c2cb09bf7adff038e1e6c48bb20 *src/dsCMatrix.h +17b6231982b96f025929e095261f4f63 *src/dsCMatrix.c +dadce6c58bc388ff5f6501f4c3cafab1 *src/dsCMatrix.h 417baab279fa6ac1bbc3679201cebb56 *src/dspMatrix.c 461e5e7cafb6a6cc4e466f124c0ef54e *src/dspMatrix.h -7bcf5e2a660ebf3d37ffcc97ee2cad6d *src/dsyMatrix.c +8cdfdb9cff2e8f123a2f2a92fc6147cb *src/dsyMatrix.c 6bdd2595738cb9fc23d88a54e3a9cc2d *src/dsyMatrix.h 92c8eec17007585bce9e8762f51ff346 *src/dtCMatrix.c 0eab9c4475b8630bb8f2a3ebb782530e *src/dtCMatrix.h @@ -409,7 +410,7 @@ 27b670c0d4e3495642f5d55601e256aa *src/dtrMatrix.h f3c82d4a7faeaf3d820ea313093a563c *src/factorizations.c 23c4dbba1f4d18ef2ff794ad129e8174 *src/factorizations.h -b9a5e2718827c707ce96caac4f507c1e *src/init.c +a84fa6cb280a67accb2fdcaa35882bf3 *src/init.c 2582450ecf52d4e9d9fece3695fa7af2 *src/ldense.c d8522bb27a18eea16d17c411d8da9e38 *src/ldense.h 0b0a474b5b19cb09637063134a0ba203 *src/lgCMatrix.c @@ -427,19 +428,19 @@ 4399d656e1aec423a0b43d2b5ce9e1a7 *src/t_gCMatrix_colSums.c e6de21ff3bf0908eaca620d9b868b4f1 *src/t_sparseVector.c e0f570b1cb86d81cb278f2cc719c2046 *tests/Class+Meth.R -e6d1dcf18839f6614e22f85fb3fc82d5 *tests/Simple.R +29c9f24b0b47f403c17f37c488ecc252 *tests/Simple.R 2af48af7c30c89de75a11153cd6aa367 *tests/abIndex-tsts.R ce85ee2aba7d78e88fef738a0d1c7118 *tests/base-matrix-fun.R 1e92c4f780ea8228568b3134e3e968e9 *tests/bind.R db1c3877f76089d81185235015a8c6fd *tests/bind.Rout.save 21024c2ae5eb4f5b1e073a7568171d89 *tests/dg_Matrix.R -cec0c2addf89c41fa06f8db635b1dd29 *tests/dpo-test.R +1f5b88045e501047f8adf08c4cb542a8 *tests/dpo-test.R 25542600eadca4e8af5261e98e3b3c5d *tests/dtpMatrix.R -479f216457efd9c471cc19ed1d1ef1dd *tests/factorizing.R +5f523cc4ed7d6ad711c36261b9227c2b *tests/factorizing.R 37e0c707594f93706ae106721f5d8bbd *tests/group-methods.R 353eff02f0ba5bfd4788845935f73ad0 *tests/indexing.R 8ad1a5ba71c94350d489718164e974f5 *tests/indexing.Rout.save -680500b65db259673c4aef94731e87b8 *tests/matprod.R +85916bde090155b7ff4f7c09ebaff471 *tests/matprod.R 06d9de87597574aff0aa03776011bfa8 *tests/matr-exp.R add6449b56011b1b8f8db4f29cdb7306 *tests/other-pkgs.R 2787ae0732f209eee2a0965c196c3363 *tests/spModel.matrix.R diff -Nru rmatrix-1.1-3/NAMESPACE rmatrix-1.1-4/NAMESPACE --- rmatrix-1.1-3/NAMESPACE 2014-03-26 18:17:09.000000000 +0000 +++ rmatrix-1.1-4/NAMESPACE 2014-04-26 19:14:30.000000000 +0000 @@ -32,7 +32,7 @@ "Hilbert", "KhatriRao", "Matrix", "spMatrix", - "sparseMatrix", + "sparseMatrix", "rsparsematrix", "Schur", "abIseq", "abIseq1", "rep2abI", "band", diff -Nru rmatrix-1.1-3/R/AllClass.R rmatrix-1.1-4/R/AllClass.R --- rmatrix-1.1-3/R/AllClass.R 2013-10-21 10:02:21.000000000 +0000 +++ rmatrix-1.1-4/R/AllClass.R 2014-04-26 19:14:30.000000000 +0000 @@ -635,7 +635,7 @@ validity = function(object) .Call(LU_validate, object)) setClass("sparseLU", contains = "LU", - representation(L = "dgCMatrix", U = "dgCMatrix", + representation(L = "dtCMatrix", U = "dtCMatrix", p = "integer", q = "integer")) ##--- QR --- @@ -645,7 +645,7 @@ p = "integer", R = "dgCMatrix", q = "integer"), validity = function(object) .Call(sparseQR_validate, object)) -##-- "SPQR" ---> ./spqr.R for noew +##-- "SPQR" ---> ./spqr.R for now ## "denseQR" -- ? (``a version of'' S3 class "qr") diff -Nru rmatrix-1.1-3/R/Auxiliaries.R rmatrix-1.1-4/R/Auxiliaries.R --- rmatrix-1.1-3/R/Auxiliaries.R 2013-12-19 13:52:00.000000000 +0000 +++ rmatrix-1.1-4/R/Auxiliaries.R 2014-06-14 15:50:25.000000000 +0000 @@ -589,48 +589,73 @@ ## nr= nrow: since i in {0,1,.., nrow-1} these are 1L "decimal" encodings: ## Further, these map to and from the usual "Fortran-indexing" (but 0-based) -## __ care against integer overflow __ encodeInd <- function(ij, di) { stopifnot(length(di) == 2) nr <- di[1L] + ## __ care against integer overflow __ if(prod(di) >= .Machine$integer.max) nr <- as.double(nr) ij[,1] + ij[,2] * nr } encodeInd2 <- function(i,j, di) { stopifnot(length(di) == 2) nr <- di[1L] + ## __ care against integer overflow __ if(prod(di) >= .Machine$integer.max) nr <- as.double(nr) i + j * nr } -}## no more needed +} else { +##' Encode Matrix index (i,j) |--> i + j * nrow {i,j : 0-origin} +##' +##' @param ij 2-column integer matrix +##' @param dim dim(.), i.e. length 2 integer vector +##' @param checkBnds logical indicating 0 <= ij[,k] < di[k] need to be checked. +##' +##' @return encoded index; integer if prod(dim) is small; double otherwise +encodeInd <- function(ij, dim, orig1=FALSE, checkBnds=TRUE) + .Call(m_encodeInd, ij, dim, orig1, checkBnds) +## --> in ../src/Mutils.c : m_encodeInd(ij, di, orig_1, chk_bnds) +## ~~~~~~~~~~~ + +##' Here, 1-based indices (i,j) are default: +encodeInd2 <- function(i, j, dim, orig1=TRUE, checkBnds=TRUE) + .Call(m_encodeInd2, i,j, dim, orig1, checkBnds) + +} -## 'code' is 0-based; as.integer(.) because encode*() may produce double +##' Decode "encoded" (i,j) indices back to cbind(i,j) +##' This is the inverse of encodeInd(.) +##' +##' @title Decode "Encoded" (i,j) Indices +##' @param code integer in 0:((n x m - 1) <==> encodeInd(.) result +##' @param nr the number of rows +##' @return +##' @author Martin Maechler decodeInd <- function(code, nr) cbind(as.integer(code %% nr), - as.integer(code %/% nr)) + as.integer(code %/% nr), deparse.level=0L) -complementInd <- function(ij, dim, checkBnds=FALSE) +complementInd <- function(ij, dim, orig1=FALSE, checkBnds=FALSE) { ## Purpose: Compute the complement of the 2-column 0-based ij-matrix ## but as 1-based indices n <- prod(dim) if(n == 0) return(integer(0)) - seq_len(n)[-(1L + .Call(m_encodeInd, ij, dim, checkBnds))] + seq_len(n)[-(1L + .Call(m_encodeInd, ij, dim, orig1, checkBnds))] } unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2)) -intersectInd <- function(ij1, ij2, di, checkBnds=FALSE) { +intersectInd <- function(ij1, ij2, di, orig1=FALSE, checkBnds=FALSE) { ## from 2-column (i,j) matrices where i in {0,.., nrow-1}, ## return only the *common* entries - decodeInd(intersect(.Call(m_encodeInd, ij1, di, checkBnds), - .Call(m_encodeInd, ij2, di, checkBnds)), nr=di[1]) + decodeInd(intersect(.Call(m_encodeInd, ij1, di, orig1, checkBnds), + .Call(m_encodeInd, ij2, di, orig1, checkBnds)), nr=di[1]) } -WhichintersectInd <- function(ij1, ij2, di, checkBnds=FALSE) { +WhichintersectInd <- function(ij1, ij2, di, orig1=FALSE, checkBnds=FALSE) { ## from 2-column (i,j) matrices where i \in {0,.., nrow-1}, ## find *where* common entries are in ij1 & ij2 - m1 <- match(.Call(m_encodeInd, ij1, di, checkBnds), - .Call(m_encodeInd, ij2, di, checkBnds)) + m1 <- match(.Call(m_encodeInd, ij1, di, orig1, checkBnds), + .Call(m_encodeInd, ij2, di, orig1, checkBnds)) ni <- !is.na(m1) list(which(ni), m1[ni]) } @@ -689,11 +714,11 @@ ## is 'x' a uniq Tsparse Matrix ? is_not_uniqT <- function(x, di = dim(x)) - is.unsorted(x@j) || anyDuplicated(.Call(m_encodeInd2, x@i, x@j, di, FALSE)) + is.unsorted(x@j) || anyDuplicated(.Call(m_encodeInd2, x@i, x@j, di, FALSE, FALSE)) ## is 'x' a TsparseMatrix with duplicated entries (to be *added* for uniq): is_duplicatedT <- function(x, di = dim(x)) - anyDuplicated(.Call(m_encodeInd2, x@i, x@j, di, FALSE)) + anyDuplicated(.Call(m_encodeInd2, x@i, x@j, di, FALSE, FALSE)) t_geMatrix <- function(x) { @@ -1312,7 +1337,8 @@ .dgC.0.factors(x) } -.set.factors <- function(x, name, value) .Call(R_set_factors, x, value, name) +.set.factors <- function(x, name, value, warn.no.slot=FALSE) + .Call(R_set_factors, x, value, name, warn.no.slot) ### Fast, much simplified version of tapply() tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) { @@ -1416,16 +1442,20 @@ int = if(n1 < n2) y[m1] else x[m2]) } -##' @title Warn about extraneous arguments in the "..." (of its caller) -## TODO NOTE: R/src/library/base/R/seq.R uses a simpler approach, -## *and* ngettext(.) {-> singular/plural *and* translatable} -##' @return -##' @author Martin Maechler, June 2012 -chk.s <- function(...) { - if(length(list(...))) - warning(gettextf("arguments %s are disregarded in\n %s", - sub(")$", '', sub("^list\\(", '', deparse(list(...), control=c()))), - deparse(sys.call(-1), control=c())), call. = FALSE, domain=NA) +##' @title Warn about extraneous arguments in the "..." (of its caller). +##' A merger of my approach and the one in seq.default() +##' @author Martin Maechler, June 2012, May 2014 +##' @param ... +##' @param which.call passed to sys.call(). A caller may use -2 if the message should +##' mention *its* caller +chk.s <- function(..., which.call = -1) { + if(nx <- length(list(...))) + warning(sprintf(ngettext(nx, + "extra argument %s will be disregarded in\n %s", + "extra arguments %s will be disregarded in\n %s"), + sub(")$", '', sub("^list\\(", '', deparse(list(...), control=c()))), + deparse(sys.call(which.call), control=c())), + call. = FALSE, domain=NA) } ##' *Only* to be used as function in diff -Nru rmatrix-1.1-3/R/CHMfactor.R rmatrix-1.1-4/R/CHMfactor.R --- rmatrix-1.1-3/R/CHMfactor.R 2013-09-02 15:48:30.000000000 +0000 +++ rmatrix-1.1-4/R/CHMfactor.R 2014-06-14 15:50:25.000000000 +0000 @@ -1,10 +1,10 @@ -### FIXME: We really want the separate parts (P,L,D) of A = P' L D L' P' -### ----- --> ~/R/MM/Pkg-ex/Matrix/chol-ex.R ---------------- -setAs("CHMfactor", "sparseMatrix", - function(from) .Call(CHMfactor_to_sparse, from)) - -setAs("CHMfactor", "Matrix", function(from) as(from, "sparseMatrix")) +### TODO: We really want the separate parts (P,L,D) of A = P' L D L' P +### --- --> ~/R/MM/Pkg-ex/Matrix/chol-ex.R --------------- +## but we currently only get A = P' L L' P --- now documented in ../man/Cholesky.Rd +setAs("CHMfactor", "sparseMatrix", function(from) .Call(CHMfactor_to_sparse, from)) +setAs("CHMfactor", "triangularMatrix", function(from) .Call(CHMfactor_to_sparse, from)) +setAs("CHMfactor", "Matrix", function(from) .Call(CHMfactor_to_sparse, from)) setAs("CHMfactor", "pMatrix", function(from) as(from@perm + 1L, "pMatrix")) @@ -33,7 +33,7 @@ system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), ...) { - chk.s(...) + chk.s(..., which.call=-2) sysDef <- eval(formals()$system) .Call(CHMfactor_solve, ##-> cholmod_solve() in ../src/CHOLMOD/Cholesky/cholmod_solve.c @@ -58,7 +58,7 @@ function(a, b, system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), ...) { - chk.s(...) + chk.s(..., which.call=-2) sysDef <- eval(formals()$system) .Call(CHMfactor_spsolve, #--> cholmod_spsolve() in ../src/CHOLMOD/Cholesky/cholmod_spsolve.c a, as(as(b, "CsparseMatrix"), "dgCMatrix"), @@ -73,7 +73,7 @@ function(a, b, system = c("A", "LDLt", "LD","DLt", "L","Lt", "D", "P","Pt"), ...) { - chk.s(...) + chk.s(..., which.call=-2) sysDef <- eval(formals()$system) system <- match.arg(system, sysDef) i.sys <- match(system, sysDef, nomatch = 0L) @@ -94,7 +94,7 @@ setMethod("chol2inv", signature(x = "CHMfactor"), function (x, ...) { - chk.s(...) + chk.s(..., which.call=-2) solve(x, system = "A") }) @@ -121,7 +121,7 @@ !is.null(v <- getOption("Matrix.verbose")) && v >= 1) message(gettextf("Quadratic matrix '%s' (=: A) is not formally\n symmetric. Will be treated as A A' ", "parent"), domain=NA) - chk.s(...) + chk.s(..., which.call=-2) .Call(CHMfactor_update, object, parent, mult) }) ##' fast version, somewhat hidden; here parent *must* be 'd[sg]CMatrix' diff -Nru rmatrix-1.1-3/R/Csparse.R rmatrix-1.1-4/R/Csparse.R --- rmatrix-1.1-3/R/Csparse.R 2014-03-26 18:16:16.000000000 +0000 +++ rmatrix-1.1-4/R/Csparse.R 2014-04-26 19:14:30.000000000 +0000 @@ -348,8 +348,8 @@ if(has.x && sum(sel) == lenRepl) { ## all entries to be replaced are non-zero: ## need indices instead of just 'sel', for, e.g., A[2:1, 2:1] <- v non0 <- cbind(match(x@i[sel], i1), - match(xj [sel], i2)) - 1L - iN0 <- 1L + .Call(m_encodeInd, non0, di = dind, checkBounds = FALSE) + match(xj [sel], i2), deparse.level=0L) + iN0 <- 1L + .Call(m_encodeInd, non0, di = dind, orig1=TRUE, checkBounds=FALSE) has0 <- if(spV) length(value@i) < lenV else any(value[!is.na(value)] == 0) diff -Nru rmatrix-1.1-3/R/dgCMatrix.R rmatrix-1.1-4/R/dgCMatrix.R --- rmatrix-1.1-3/R/dgCMatrix.R 2013-09-14 17:09:49.000000000 +0000 +++ rmatrix-1.1-4/R/dgCMatrix.R 2014-06-14 15:50:25.000000000 +0000 @@ -114,7 +114,7 @@ if(sparse) .solve.sparse.dgC(a, b, tol=tol) else .Call(dgCMatrix_matrix_solve, a, b, FALSE) .solve.dgC.mat <- function(a, b, sparse=FALSE, tol = .Machine$double.eps, ...) { - chk.s(...) + chk.s(..., which.call=-2) if(sparse) .solve.sparse.dgC(a, b, tol=tol) else .Call(dgCMatrix_matrix_solve, a, b, FALSE) } @@ -123,11 +123,11 @@ setMethod("solve", signature(a = "dgCMatrix", b = "dsparseMatrix"), function(a, b, sparse=NA, tol = .Machine$double.eps, ...) { - chk.s(...) + chk.s(..., which.call=-2) if(is.na(sparse)) { if(isSymmetric(a)) ## TODO: fast cholmod_symmetric() for Cholesky - return(solve(forceCspSymmetric(a, isTri=FALSE), b)) + return(solve(forceCspSymmetric(a, isTri=FALSE), b, tol=tol)) #-> sparse result ## else sparse <- FALSE # (old default) @@ -140,7 +140,7 @@ ## (MM: a bit less dumb now with possibility of staying sparse) setMethod("solve", signature(a = "dgCMatrix", b = "missing"), function(a, b, sparse=NA, tol = .Machine$double.eps, ...) { - chk.s(...) + chk.s(..., which.call=-2) if(is.na(sparse)) { if(isSymmetric(a)) ## TODO: fast cholmod_symmetric() for Cholesky diff -Nru rmatrix-1.1-3/R/dgTMatrix.R rmatrix-1.1-4/R/dgTMatrix.R --- rmatrix-1.1-3/R/dgTMatrix.R 2014-03-12 18:44:33.000000000 +0000 +++ rmatrix-1.1-4/R/dgTMatrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -62,7 +62,7 @@ ## "[<-" methods { setReplaceMethod()s } too ... -setMethod("image", "dgTMatrix", +setMethod("image", "dgTMatrix", ## *The* real one function(x, xlim = c(1, di[2]), ylim = c(di[1], 1), aspect = "iso", diff -Nru rmatrix-1.1-3/R/diagMatrix.R rmatrix-1.1-4/R/diagMatrix.R --- rmatrix-1.1-3/R/diagMatrix.R 2014-03-26 18:17:09.000000000 +0000 +++ rmatrix-1.1-4/R/diagMatrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -682,7 +682,7 @@ dx <- dim(x) dy <- dim(y <- .Call(Csparse_diagU2N, y)) if(dx[2] != dy[1]) stop("non-matching dimensions") - if(x@diag == "N") { + if(x@diag == "N") { if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix")) y <- as(y, "generalMatrix") y@x <- y@x * x@x[y@i + 1L] @@ -965,8 +965,7 @@ f0 <- callGeneric(0, e2) if(all0(f0)) { # remain diagonal L1 <- (le <- length(e2)) == 1L - E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE - ## E <- copyClass(e1, "ddiMatrix", check=FALSE) + E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE) ## storage.mode(E@x) <- "double" if(e1@diag == "U") { if(any((r <- callGeneric(1, e2)) != 1)) { @@ -990,8 +989,7 @@ f0 <- callGeneric(e1, 0) if(all0(f0)) { # remain diagonal L1 <- (le <- length(e1)) == 1L - E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE - ## E <- copyClass(e2, "ddiMatrix", check=FALSE) + E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE) ## storage.mode(E@x) <- "double" if(e2@diag == "U") { if(any((r <- callGeneric(e1, 1)) != 1)) { @@ -1023,8 +1021,7 @@ f0 <- callGeneric(0, e2) if(all0(f0)) { # remain diagonal L1 <- (le <- length(e2)) == 1L - E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE - ## E <- copyClass(e1, "ldiMatrix", check=FALSE) + E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE) ## storage.mode(E@x) <- "logical" if(e1@diag == "U") { if(any((r <- callGeneric(1, e2)) != 1)) { diff -Nru rmatrix-1.1-3/R/dpoMatrix.R rmatrix-1.1-4/R/dpoMatrix.R --- rmatrix-1.1-3/R/dpoMatrix.R 2013-09-23 20:36:59.000000000 +0000 +++ rmatrix-1.1-4/R/dpoMatrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -4,12 +4,12 @@ function(from) copyClass(.Call(dsyMatrix_as_dspMatrix, from), "dppMatrix", - sNames = c("x", "Dim", "Dimnames", "uplo", "factors"))) + sNames = c("x", "Dim", "Dimnames", "uplo", "factors")))#FIXME , check=FALSE setAs("dpoMatrix", "corMatrix", function(from) { if(!is.null(cm <- from@factors$correlation)) return(cm) - ## else + ## else sd <- sqrt(diag(from)) if(is.null(names(sd)) && !is.null(nms <- from@Dimnames[[1]])) names(sd) <- nms diff -Nru rmatrix-1.1-3/R/dppMatrix.R rmatrix-1.1-4/R/dppMatrix.R --- rmatrix-1.1-3/R/dppMatrix.R 2013-08-16 14:09:01.000000000 +0000 +++ rmatrix-1.1-4/R/dppMatrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -4,7 +4,7 @@ function(from) copyClass(.Call(dspMatrix_as_dsyMatrix, from), "dpoMatrix", - sNames = c("x", "Dim", "Dimnames", "uplo", "factors"))) + sNames = c("x", "Dim", "Dimnames", "uplo", "factors")))#FIXME , check=FALSE dpp2sC <- function(from) as(.Call(dspMatrix_as_dsyMatrix, from), "dsCMatrix") ## setAs("dppMatrix", "dsCMatrix", dpp2sC) setAs("dppMatrix", "CsparseMatrix", dpp2sC) @@ -27,7 +27,7 @@ stop("not a positive definite matrix") ## else copyClass(from, "dppMatrix", - sNames = c("x", "Dim", "Dimnames", "uplo", "factors")) + sNames = c("x", "Dim", "Dimnames", "uplo", "factors"))#FIXME , check=FALSE }) diff -Nru rmatrix-1.1-3/R/dsCMatrix.R rmatrix-1.1-4/R/dsCMatrix.R --- rmatrix-1.1-3/R/dsCMatrix.R 2013-09-14 17:09:49.000000000 +0000 +++ rmatrix-1.1-4/R/dsCMatrix.R 2014-04-12 21:37:37.000000000 +0000 @@ -75,8 +75,8 @@ else triu(as(x, "dgCMatrix"), k = k, ...) }) -solve.dsC.mat <- function(a,b, tol = .Machine$double.eps) { - r <- tryCatch(.Call(dsCMatrix_matrix_solve, a, b), +solve.dsC.mat <- function(a,b, LDL = NA, tol = .Machine$double.eps) { + r <- tryCatch(.Call(dsCMatrix_matrix_solve, a, b, LDL), error=function(e)NULL, warning=function(w)NULL) if(is.null(r)) { ## cholmod factorization was not ok Matrix.msg("solve.dsC.mat(): Cholmod factorization unsuccessful --> using LU()") @@ -86,8 +86,8 @@ } ## ``Fully-sparse'' solve() {different Cholmod routine, otherwise "the same"}: -solve.dsC.dC <- function(a,b, tol = .Machine$double.eps) { - r <- tryCatch(.Call(dsCMatrix_Csparse_solve, a, b), +solve.dsC.dC <- function(a,b, LDL = NA, tol = .Machine$double.eps) { + r <- tryCatch(.Call(dsCMatrix_Csparse_solve, a, b, LDL), error=function(e)NULL, warning=function(w)NULL) if(is.null(r)) { ## cholmod factorization was not ok Matrix.msg("solve.dsC.dC(): Cholmod factorization unsuccessful --> using LU()") @@ -99,37 +99,38 @@ ## . ------------------------ setMethod("solve", signature(a = "dsCMatrix", b = "ddenseMatrix"), - function(a, b, ...) { + function(a, b, LDL = NA, tol = .Machine$double.eps, ...) { if (class(b) != "dgeMatrix") b <- .Call(dup_mMatrix_as_dgeMatrix, b) - solve.dsC.mat(a,b) + solve.dsC.mat(a,b, LDL=LDL, tol=tol) }, valueClass = "dgeMatrix") setMethod("solve", signature(a = "dsCMatrix", b = "denseMatrix"), ## only triggers for diagonal*, ldense*.. (but *not* ddense: above) - function(a, b, ...) - solve.dsC.mat(a, as(.Call(dup_mMatrix_as_geMatrix, b), "dgeMatrix"))) + function(a, b, LDL = NA, tol = .Machine$double.eps, ...) + solve.dsC.mat(a, as(.Call(dup_mMatrix_as_geMatrix, b), "dgeMatrix"), + LDL=LDL, tol=tol)) setMethod("solve", signature(a = "dsCMatrix", b = "matrix"), - function(a, b, ...) - solve.dsC.mat(a, .Call(dup_mMatrix_as_dgeMatrix, b)), + function(a, b, LDL = NA, tol = .Machine$double.eps, ...) + solve.dsC.mat(a, .Call(dup_mMatrix_as_dgeMatrix, b), LDL=LDL, tol=tol), valueClass = "dgeMatrix") setMethod("solve", signature(a = "dsCMatrix", b = "numeric"), - function(a, b, ...) - solve.dsC.mat(a, .Call(dup_mMatrix_as_dgeMatrix, b)), + function(a, b, LDL = NA, tol = .Machine$double.eps, ...) + solve.dsC.mat(a, .Call(dup_mMatrix_as_dgeMatrix, b), LDL=LDL, tol=tol), valueClass = "dgeMatrix") ## . ------------------------ setMethod("solve", signature(a = "dsCMatrix", b = "dsparseMatrix"), - function(a, b, ...) { + function(a, b, LDL = NA, tol = .Machine$double.eps, ...) { cb <- getClassDef(class(b)) if (!extends(cb, "CsparseMatrix")) cb <- getClassDef(class(b <- as(b, "CsparseMatrix"))) if (extends(cb, "symmetricMatrix")) ## not supported (yet) by cholmod_spsolve b <- as(b, "dgCMatrix") - solve.dsC.dC(a,b) + solve.dsC.dC(a,b, LDL=LDL, tol=tol) }) setMethod("solve", signature(a = "dsCMatrix", b = "missing"), @@ -154,13 +155,16 @@ ### These two are very similar, the first one has the advantage to be applicable to 'Chx' directly: +## "used" currently only in ../tests/factorizing.R .diag.dsC <- function(x, Chx = Cholesky(x, LDL=TRUE), res.kind = "diag") { force(Chx) if(!missing(Chx)) stopifnot(.isLDL(Chx), is.integer(Chx@p), is.double(Chx@x)) - .Call(diag_tC, Chx@p, Chx@x, Chx@perm, res.kind) + .Call(diag_tC, Chx, res.kind) ## ^^^^^^^ from ../src/Csparse.c + ## => res.kind in ("trace", "sumLog", "prod", "min", "max", "range", "diag", "diagBack") } +## nowhere used/tested (FIXME?) ## here, we *could* allow a 'mult = 0' factor : .CHM.LDL.D <- function(x, perm = TRUE, res.kind = "diag") { .Call(dsCMatrix_LDL_D, x, perm, res.kind) @@ -181,7 +185,7 @@ ## these are slightly faster (ca. 3 to 4 %): ldet2.dsC <- function(x, ...) { Ch <- Cholesky(x, super = FALSE, ...) - .Call(diag_tC, Ch@p, Ch@x, Ch@perm, "sumLog") + .Call(diag_tC, Ch, "sumLog") } ## only very slightly ( ~ < 1% ) faster (than "ldet2"): ldet3.dsC <- function(x, perm = TRUE) @@ -206,7 +210,7 @@ if(is.null(Chx)) ## we do *not* have a positive definite matrix detSparseLU(x, logarithm) else { - d <- .Call(diag_tC, Chx@p, Chx@x, Chx@perm, res.kind = "diag") + d <- .Call(diag_tC, Chx, res.kind = "diag") mkDet(d, logarithm=logarithm) } }) diff -Nru rmatrix-1.1-3/R/dtCMatrix.R rmatrix-1.1-4/R/dtCMatrix.R --- rmatrix-1.1-3/R/dtCMatrix.R 2013-08-08 21:06:22.000000000 +0000 +++ rmatrix-1.1-4/R/dtCMatrix.R 2014-04-26 19:14:30.000000000 +0000 @@ -106,3 +106,29 @@ function(a, b, ...) .Call(dtCMatrix_matrix_solve, a, as.matrix(as.double(b)), FALSE), valueClass = "dgeMatrix") + +if(FALSE)## still not working +setMethod("diag", "dtCMatrix", + function(x, nrow, ncol) .Call(diag_tC, x, "diag")) + + +## no pivoting here, use L or U +setMethod("lu", "dtCMatrix", + function(x, ...) { + n <- (d <- x@Dim)[1L] + p <- 0:(n-1L) + if(x@uplo == "U") + new("sparseLU", + L = .trDiagonal(n, uplo="L"), + U = x, + p = p, q = p, Dim = d) + else { ## "L" : x = L = L I + d <- diag(x) + new("sparseLU", + L = x %*% Diagonal(n, 1/d), + U = .trDiagonal(n, x = d), + p = p, q = p, Dim = d) + } + }) + + diff -Nru rmatrix-1.1-3/R/dtrMatrix.R rmatrix-1.1-4/R/dtrMatrix.R --- rmatrix-1.1-3/R/dtrMatrix.R 2012-07-24 14:02:28.000000000 +0000 +++ rmatrix-1.1-4/R/dtrMatrix.R 2014-06-14 15:50:25.000000000 +0000 @@ -71,7 +71,7 @@ setMethod("chol2inv", signature(x = "dtrMatrix"), function (x, ...) { - chk.s(...) + chk.s(..., which.call=-2) if (x@diag != "N") x <- diagU2N(x) .Call(dtrMatrix_chol2inv, x) }) @@ -79,31 +79,31 @@ setMethod("solve", signature(a = "dtrMatrix", b="missing"), function(a, b, ...) { ## warn, as e.g. CHMfactor have 'system' as third argument - chk.s(...) + chk.s(..., which.call=-2) .Call(dtrMatrix_solve, a) }, valueClass = "dtrMatrix") setMethod("solve", signature(a = "dtrMatrix", b="ddenseMatrix"), function(a, b, ...) { - chk.s(...) + chk.s(..., which.call=-2) .Call(dtrMatrix_matrix_solve, a, b) }, valueClass = "dgeMatrix") setMethod("solve", signature(a = "dtrMatrix", b="dMatrix"), function(a, b, ...) { - chk.s(...) + chk.s(..., which.call=-2) .Call(dtrMatrix_matrix_solve, a, as(b,"denseMatrix")) }, valueClass = "dgeMatrix") setMethod("solve", signature(a = "dtrMatrix", b="Matrix"), function(a, b, ...) { - chk.s(...) + chk.s(..., which.call=-2) .Call(dtrMatrix_matrix_solve, a, as(as(b, "dMatrix"), "denseMatrix")) }, valueClass = "dgeMatrix") setMethod("solve", signature(a = "dtrMatrix", b="matrix"), function(a, b, ...) { - chk.s(...) + chk.s(..., which.call=-2) .Call(dtrMatrix_matrix_solve, a, b) }, valueClass = "dgeMatrix") diff -Nru rmatrix-1.1-3/R/indMatrix.R rmatrix-1.1-4/R/indMatrix.R --- rmatrix-1.1-3/R/indMatrix.R 2013-09-04 08:30:53.000000000 +0000 +++ rmatrix-1.1-4/R/indMatrix.R 2014-06-14 16:21:38.000000000 +0000 @@ -108,10 +108,16 @@ setMethod("crossprod", signature(x = "indMatrix", y = "indMatrix"), function(x, y) { mmultCheck(x,y, 2L) - xy <- interaction(x@perm, y@perm) + ## xy <- interaction(x@perm, y@perm) + ## this is wrong if any of the columns in X or Y are empty because interaction() + ## drops non-occuring levels from a non-factor. Explicitly defining a factor with + ## levels 1:ncol() avoids that. + nx <- x@Dim[2L] + ny <- y@Dim[2L] + xy <- interaction(factor(x@perm, levels=seq_len(nx)), + factor(y@perm, levels=seq_len(ny))) Matrix(data= tabulate(xy, nbins=nlevels(xy)), - nrow= x@Dim[2L], - ncol= y@Dim[2L]) + nrow = nx, ncol = ny) }) setMethod("tcrossprod", signature(x = "matrix", y = "indMatrix"), @@ -131,9 +137,18 @@ setMethod("kronecker", signature(X = "indMatrix", Y = "indMatrix"), function (X, Y, FUN = "*", make.dimnames = FALSE, ...) { if (FUN != "*") stop("kronecker method must use default 'FUN'") - perm <- as.integer(interaction(rep(X@perm, each =Y@Dim[1]), - rep(Y@perm, times=X@Dim[1]), - lex.order=TRUE)) +## perm <- as.integer(interaction(rep(X@perm, each =Y@Dim[1]), +## rep(Y@perm, times=X@Dim[1]), +## lex.order=TRUE)) + ## this is wrong if any of the columns in X or Y are empty because interaction() + ## drops non-occuring levels from a non-factor. Explicitly defining a factor + ## with levels 1:ncol(.) avoids that. + perm <- as.integer(interaction(factor(rep(X@perm, each =Y@Dim[1]), + levels=seq_len(X@Dim[2])), + factor(rep.int(Y@perm, times=X@Dim[1]), + levels=seq_len(Y@Dim[2])), + lex.order=TRUE)) + new("indMatrix", perm=perm, Dim=X@Dim*Y@Dim) }) diff -Nru rmatrix-1.1-3/R/Matrix.R rmatrix-1.1-4/R/Matrix.R --- rmatrix-1.1-3/R/Matrix.R 2014-03-27 09:30:57.000000000 +0000 +++ rmatrix-1.1-4/R/Matrix.R 2014-06-14 15:50:25.000000000 +0000 @@ -294,12 +294,12 @@ function (x, ...) chol2inv(as(as(x, "dMatrix"), "dtrMatrix"), ...)) setMethod("chol2inv", signature(x = "diagonalMatrix"), function (x, ...) { - chk.s(...) + chk.s(..., which.call=-2) tcrossprod(solve(x)) }) setMethod("chol2inv", signature(x = "sparseMatrix"), function (x, ...) { - chk.s(...) + chk.s(..., which.call=-2) ## for now: tcrossprod(solve(as(x,"triangularMatrix"))) }) @@ -483,8 +483,8 @@ message(gettextf("in Summary(, .): %s(<%s>, <%s>,...)\n", .Generic, class(x), class(a[[1]])), domain = NA) else - message(sprintf("in Summary(, .): %s(<%s>, <%s>)\n", - .Generic, class(x), class(a[[1]])), domain = NA) + message(gettextf("in Summary(, .): %s(<%s>, <%s>)\n", + .Generic, class(x), class(a[[1]])), domain = NA) do.call(.Generic, c(x, a, list(na.rm=na.rm))) }}) @@ -633,8 +633,8 @@ ij.x <- non0.i(x, cld) } - m1 <- .Call(m_encodeInd, ij.x, di, checkBounds = FALSE) - m2 <- .Call(m_encodeInd, ij -1L, di, checkBounds = TRUE) + m1 <- .Call(m_encodeInd, ij.x, di, orig1=FALSE, checkBounds=FALSE) + m2 <- .Call(m_encodeInd, ij, di, orig1= TRUE, checkBounds= TRUE) mi <- match(m2, m1, nomatch=0) mmi <- mi != 0L ## == (m2 %in% m1) ## Result: all FALSE or 0 apart from where we match non-zero entries diff -Nru rmatrix-1.1-3/R/Ops.R rmatrix-1.1-4/R/Ops.R --- rmatrix-1.1-3/R/Ops.R 2014-03-26 18:16:16.000000000 +0000 +++ rmatrix-1.1-4/R/Ops.R 2014-04-26 19:14:30.000000000 +0000 @@ -240,7 +240,7 @@ ## Here, we assume that 'r' and the indices align (!) encI <- .Call(m_encodeInd, non0ind(e1, cl1, uniqT=FALSE, xtendSymm=FALSE), - di = d, checkBounds = FALSE) + di = d, orig1=FALSE, checkBounds=FALSE) rx[1L + encI] <- r r <- new(lClass, x = rx, Dim = d, Dimnames = dimnames(e1)) } @@ -744,7 +744,7 @@ ## Here, we assume that 'r' and the indices align (!) encI <- .Call(m_encodeInd, non0ind(e1, cl1, uniqT=FALSE, xtendSymm=FALSE), - di = d, checkBounds = FALSE) + di = d, orig1=FALSE, checkBounds=FALSE) rx[1L + encI] <- r r <- new(lClass, x = rx, Dim = d, Dimnames = dimnames(e1)) } @@ -1115,6 +1115,7 @@ setMethod("Arith", signature(e1 = "CsparseMatrix", e2 = "numeric"), .Arith.CM.atom) setMethod("Arith", signature(e1 = "numeric", e2 = "CsparseMatrix"), .Arith.atom.CM) +##' compute indices for recycling of length 'len' to match sparseMatrix 'spM' .Ops.recycle.ind <- function(spM, len) { n <- prod(d <- dim(spM)) if(n < len) stop("vector too long in Matrix - vector operation") @@ -1122,7 +1123,7 @@ warning("longer object length\n\tis not a multiple of shorter object length") ## TODO(speedup!): construction of [1L + in0 %%len] via one .Call() in0 <- .Call(m_encodeInd, .Call(compressed_non_0_ij, spM, TRUE), - d, FALSE) + d, FALSE, FALSE) 1L + in0 %% len } diff -Nru rmatrix-1.1-3/R/products.R rmatrix-1.1-4/R/products.R --- rmatrix-1.1-3/R/products.R 2014-03-04 07:35:12.000000000 +0000 +++ rmatrix-1.1-4/R/products.R 2014-06-14 15:49:51.000000000 +0000 @@ -148,7 +148,7 @@ function(x, y) .Call(Csparse_dense_prod, x, y)) setMethod("%*%", signature(x = "CsparseMatrix", y = "matrix"), function(x, y) .Call(Csparse_dense_prod, x, y)) # was x %*% Matrix(y) -setMethod("%*%", signature(x = "CsparseMatrix", y = "numeric"), +setMethod("%*%", signature(x = "CsparseMatrix", y = "numLike"), function(x, y) .Call(Csparse_dense_prod, x, y)) @@ -237,7 +237,7 @@ setMethod("crossprod", signature(x = "dgeMatrix", y = "matrix"), function(x, y = NULL) .Call(dgeMatrix_matrix_crossprod, x, y, FALSE), valueClass = "dgeMatrix") -setMethod("crossprod", signature(x = "dgeMatrix", y = "numeric"), +setMethod("crossprod", signature(x = "dgeMatrix", y = "numLike"), function(x, y = NULL) .Call(dgeMatrix_matrix_crossprod, x, y, FALSE), valueClass = "dgeMatrix") setMethod("crossprod", signature(x = "matrix", y = "dgeMatrix"), @@ -310,13 +310,12 @@ function(x, y = NULL) .Call(Csparse_Csparse_crossprod, x, y, trans = FALSE)) -## FIXME: Generalize the class of y. This specific method is to replace one -## in dgCMatrix.R +## FIXME: Generalize the class of y. (?? still ??) setMethod("crossprod", signature(x = "CsparseMatrix", y = "ddenseMatrix"), function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y)) setMethod("crossprod", signature(x = "CsparseMatrix", y = "matrix"), function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y)) -setMethod("crossprod", signature(x = "CsparseMatrix", y = "numeric"), +setMethod("crossprod", signature(x = "CsparseMatrix", y = "numLike"), function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y)) @@ -391,15 +390,15 @@ function(x, y) t(.Call(Csparse_dense_crossprod, y, x))) setMethod("crossprod", signature(x = "matrix", y = "CsparseMatrix"), function(x, y) t(.Call(Csparse_dense_crossprod, y, x))) -setMethod("crossprod", signature(x = "numeric", y = "CsparseMatrix"), +setMethod("crossprod", signature(x = "numLike", y = "CsparseMatrix"), function(x, y) t(.Call(Csparse_dense_crossprod, y, x))) ## "Matrix" : cbind(), rbind() do names -> dimnames setMethod("crossprod", signature(x = "Matrix", y = "numLike"), - function(x, y) crossprod(x, cbind(y,deparse.level=0))) -setMethod("crossprod", signature(x = "numLike", y = "Matrix"), - function(x, y) crossprod(rbind(x,deparse.level=0), y)) + function(x, y) crossprod(x, cbind(y, deparse.level=0))) +setMethod("crossprod", signature(x = "numLike", y = "Matrix"), + function(x, y) crossprod(cbind(x, deparse.level=0), y)) setMethod("crossprod", signature(x = "Matrix", y = "matrix"), function(x, y) crossprod(x, Matrix(y))) @@ -409,20 +408,37 @@ ## sparseVector setMethod("crossprod", signature(x = "Matrix", y = "sparseVector"), function(x, y) crossprod(x, .sparseV2Mat(y))) -setMethod("crossprod", signature(x = "sparseVector", y = "Matrix"), +setMethod("crossprod", signature(x = "sparseVector", y = "Matrix"), function(x, y) crossprod(spV2M(x, nrow = length(x), ncol = 1L, check = FALSE), y)) setMethod("crossprod", signature(x = "sparseVector", y = "sparseVector"), sp.x.sp) setMethod("crossprod", signature(x = "sparseVector", y = "missing"), function(x, y=NULL) sp.x.sp(x,x)) +## Fallbacks -- symmetric LHS --> saving a t(.): +## {FIXME: want the method to be `%*%` -- but primitives are not allowed as methods} +setMethod("crossprod", signature(x = "symmetricMatrix", y = "Matrix"), function(x,y) x %*% y) +setMethod("crossprod", signature(x = "symmetricMatrix", y = "missing"), function(x,y) x %*% x) +setMethod("crossprod", signature(x = "symmetricMatrix", y = "ANY"), function(x,y) x %*% y) +## ## cheap fallbacks setMethod("crossprod", signature(x = "Matrix", y = "Matrix"), - function(x, y) t(x) %*% y) + function(x, y) { + Matrix.msg(sprintf( + "potentially suboptimal crossprod(\"%s\",\"%s\") as t(.) %s y", + class(x), class(y), "%*%")) + t(x) %*% y }) setMethod("crossprod", signature(x = "Matrix", y = "missing"), - function(x, y) t(x) %*% x) + function(x, y) { + Matrix.msg(paste0( + "potentially suboptimal crossprod(<",class(x),">) as t(.) %*% . ")) + t(x) %*% x }) setMethod("crossprod", signature(x = "Matrix", y = "ANY"), - function(x, y) t(x) %*% y) + function(x, y) { + Matrix.msg(sprintf( + "potentially suboptimal crossprod(\"%s\", <%s>[=]) as t(.) %s y", + class(x), class(y), "%*%")) + t(x) %*% y }) setMethod("crossprod", signature(x = "ANY", y = "Matrix"), function(x, y) t(x) %*% y) @@ -515,9 +531,17 @@ function(x, y) .Call(Csparse_dense_prod, x, rbind(y, deparse.level=0))) +### -- xy' = (yx')' -------------------- ### FIXME (speed): Csparse_dense_crossprod should also get a 'trans = TRUE' ## so we could have one less t() in this: -## -- xy' = (yx')' +for(T in c("d", "l", "n")) { ## speedup for *symmetric* RHS + .sCMatrix <- paste0(T, "sCMatrix") + setMethod("tcrossprod", signature(x = "ddenseMatrix", y = .sCMatrix), + function(x, y) t(.Call(Csparse_dense_prod, y, x))) + setMethod("tcrossprod", signature(x = "matrix", y = .sCMatrix), + function(x, y) t(.Call(Csparse_dense_prod, y, x))) +} +rm(T, .sCMatrix) setMethod("tcrossprod", signature(x = "ddenseMatrix", y = "CsparseMatrix"), function(x, y) t(.Call(Csparse_dense_prod, y, t(x)))) setMethod("tcrossprod", signature(x = "matrix", y = "CsparseMatrix"), @@ -550,8 +574,8 @@ ## "Matrix" -setMethod("tcrossprod", signature(x = "Matrix", y = "numLike"), - function(x, y) x %*% rbind(y,deparse.level=0)) +setMethod("tcrossprod", signature(x = "Matrix", y = "numLike"), + function(x, y) x %*% rbind(y, deparse.level=0)) setMethod("tcrossprod", signature(x = "numLike", y = "Matrix"), .v.Mt) setMethod("tcrossprod", signature(x = "Matrix", y = "matrix"), function(x, y = NULL) tcrossprod(x, Matrix(y))) @@ -559,7 +583,7 @@ function(x, y = NULL) tcrossprod(Matrix(x), y)) ## sparseVector -setMethod("tcrossprod", signature(x = "Matrix", y = "sparseVector"), +setMethod("tcrossprod", signature(x = "Matrix", y = "sparseVector"), function(x, y) tcrossprod(x, .sparseV2Mat(y))) setMethod("tcrossprod", signature(x = "sparseVector", y = "Matrix"), .v.Mt) setMethod("tcrossprod", signature(x = "sparseMatrix", y = "sparseVector"), @@ -574,15 +598,31 @@ spV2M(x, nrow=1L, ncol=length(x), check=FALSE)) +## Fallbacks -- symmetric RHS --> saving a t(.): +## {FIXME: want the method to be `%*%` -- but primitives are not allowed as methods} +setMethod("tcrossprod", signature(x = "Matrix", y = "symmetricMatrix"), function(x,y) x %*% y) +setMethod("tcrossprod", signature(x = "ANY", y = "symmetricMatrix"), function(x,y) x %*% y) +## ## cheap fallbacks setMethod("tcrossprod", signature(x = "Matrix", y = "Matrix"), - function(x, y = NULL) x %*% t(y)) + function(x, y = NULL) { + Matrix.msg(sprintf( + "potentially suboptimal tcrossprod(\"%s\",\"%s\") as x %s t(y)", + class(x), class(y), "%*%")) + x %*% t(y) }) setMethod("tcrossprod", signature(x = "Matrix", y = "missing"), - function(x, y = NULL) x %*% t(x)) + function(x, y = NULL) { + Matrix.msg(paste0( + "potentially suboptimal crossprod(<",class(x), ">) as . %*% t(.)")) + x %*% t(x) }) setMethod("tcrossprod", signature(x = "Matrix", y = "ANY"), function(x, y = NULL) x %*% t(y)) setMethod("tcrossprod", signature(x = "ANY", y = "Matrix"), - function(x, y = NULL) x %*% t(y)) + function(x, y = NULL) { + Matrix.msg(sprintf( + "potentially suboptimal tcrossprod(<%s>[=], \"%s\") as x %s t(y)", + class(x), class(y), "%*%")) + x %*% t(y) }) ## Local variables: ## mode: R diff -Nru rmatrix-1.1-3/R/sparseMatrix.R rmatrix-1.1-4/R/sparseMatrix.R --- rmatrix-1.1-3/R/sparseMatrix.R 2014-01-29 09:52:48.000000000 +0000 +++ rmatrix-1.1-4/R/sparseMatrix.R 2014-04-26 19:14:30.000000000 +0000 @@ -174,7 +174,7 @@ ## lens <- vapply(e, length, 0)) } - + if(use.weights) { eWts <- graph::edgeWeights(from); names(eWts) <- NULL i <- edge2i(eWts) @@ -225,7 +225,7 @@ from <- uniqTsparse(from) if(is.null(edgemode)) - edgemode <- + edgemode <- if(isSymmetric(from)) { # either "symmetricMatrix" or otherwise ##-> undirected graph: every edge only once! if(!is(from, "symmetricMatrix")) { @@ -418,7 +418,7 @@ dimnames(cx) <- dn } if(missing(iN0)) - iN0 <- 1L + .Call(m_encodeInd, non0ind(x, cld), di = d, FALSE) + iN0 <- 1L + .Call(m_encodeInd, non0ind(x, cld), di = d, FALSE, FALSE) ## ne <- length(iN0) if(asLogical) { @@ -441,7 +441,7 @@ ij <- rbind(ij, ij[notdiag, 2:1], deparse.level=0) F. <- c(F., F.[notdiag]) } - iN0 <- 1L + .Call(m_encodeInd, ij, di = d, FALSE) + iN0 <- 1L + .Call(m_encodeInd, ij, di = d, FALSE, FALSE) cx[iN0[F.]] <- ":" # non-structural FALSE (or "o", "," , "-" or "f")? } } @@ -521,7 +521,7 @@ ## -> cannot use 'm' alone d <- dim(cx) ne <- length(iN0 <- 1L + .Call(m_encodeInd, non0ind(x, cld), - di = d, FALSE)) + di = d, FALSE, FALSE)) if(0 < ne && (logi || ne < prod(d))) { cx <- formatSparseM(x, zero.print, align, m=m, asLogical=logi, digits=digits, cx=cx, iN0=iN0, dn=dn) @@ -844,3 +844,61 @@ ### --- sparse model matrix, fac2sparse, etc ----> ./spModels.R ### xtabs(*, sparse = TRUE) ---> part of standard package 'stats' since R 2.10.0 + +##' @title Random Sparse Matrix +##' @param nrow, +##' @param ncol number of rows and columns, i.e., the matrix dimension +##' @param nnz number of non-zero entries +##' @param rand.x random number generator for 'x' slot +##' @param ... optionally further arguments passed to sparseMatrix() +##' @return a sparseMatrix of dimension (nrow, ncol) +##' @author Martin Maechler +##' @examples M1 <- rSparseMatrix(1000, 20, nnz = 200) +##' summary(M1) +if(FALSE) ## better version below +rsparsematrix <- function(nrow, ncol, nnz, + rand.x = function(n) signif(rnorm(nnz), 2), + warn.nnz = TRUE, ...) +{ + maxi.sample <- 2^31 # maximum n+1 for which sample(n) returns integer + stopifnot((nnz <- as.integer(nnz)) >= 0, + nrow >= 0, ncol >= 0, nnz <= nrow * ncol, + nrow < maxi.sample, ncol < maxi.sample) + ## to ensure that nnz is strictly followed, must act on duplicated (i,j): + i <- sample.int(nrow, nnz, replace = TRUE) + j <- sample.int(ncol, nnz, replace = TRUE) + dim <- c(nrow, ncol) + it <- 0 + while((it <- it+1) < 100 && + anyDuplicated(n.ij <- encodeInd2(i, j, dim, checkBnds=FALSE))) { + m <- length(k.dup <- which(duplicated(n.ij))) + Matrix.msg(sprintf("%3g duplicated (i,j) pairs", m), .M.level = 2) + if(runif(1) <= 1/2) + i[k.dup] <- sample.int(nrow, m, replace = TRUE) + else + j[k.dup] <- sample.int(ncol, m, replace = TRUE) + } + if(warn.nnz && it == 100 && anyDuplicated(encodeInd2(i, j, dim, checkBnds=FALSE))) + warning("number of non zeros is smaller than 'nnz' because of duplicated (i,j)s") + sparseMatrix(i = i, j = j, x = rand.x(nnz), dims = dim, ...) +} + +## No warn.nnz needed, as we sample the encoded (i,j) with*out* replacement: +rsparsematrix <- function(nrow, ncol, density, + nnz = round(density * maxE), symmetric = FALSE, + rand.x = function(n) signif(rnorm(nnz), 2), ...) +{ + maxE <- if(symmetric) nrow*(nrow+1)/2 else nrow*ncol + stopifnot((nnz <- as.integer(nnz)) >= 0, + nrow >= 0, ncol >= 0, nnz <= maxE) + ## sampling with*out* replacement (replace=FALSE !): + ijI <- -1L + + if(symmetric) sample(indTri(nrow, diag=TRUE), nnz) + else sample.int(maxE, nnz) + ## i,j below correspond to ij <- decodeInd(code, nr) : + sparseMatrix(i = ijI %% nrow, + j = ijI %/% nrow, + index1 = FALSE, symmetric = symmetric, + x = rand.x(nnz), dims = c(nrow, ncol), ...) +} + diff -Nru rmatrix-1.1-3/R/Tsparse.R rmatrix-1.1-4/R/Tsparse.R --- rmatrix-1.1-3/R/Tsparse.R 2014-03-26 18:17:09.000000000 +0000 +++ rmatrix-1.1-4/R/Tsparse.R 2014-04-26 19:14:30.000000000 +0000 @@ -339,10 +339,10 @@ cl," to ",class(x)) } nr <- di[1] - x.i <- .Call(m_encodeInd2, x@i, x@j, di=di, FALSE) + x.i <- .Call(m_encodeInd2, x@i, x@j, di=di, FALSE, FALSE) if(anyDuplicated(x.i)) { ## == if(is_duplicatedT(x, di = di)) x <- uniqTsparse(x) - x.i <- .Call(m_encodeInd2, x@i, x@j, di=di, FALSE) + x.i <- .Call(m_encodeInd2, x@i, x@j, di=di, FALSE, FALSE) } n <- prod(di) @@ -616,9 +616,10 @@ if(any(sel)) { ## the 0-based indices of non-zero entries -- WRT to submatrix - non0 <- cbind(match(x@i[sel], i1), - match(x@j[sel], i2)) - 1L - iN0 <- 1L + .Call(m_encodeInd, non0, di = dind, FALSE) + iN0 <- 1L + .Call(m_encodeInd2, + match(x@i[sel], i1), + match(x@j[sel], i2), + di = dind, orig1=TRUE, FALSE) ## 1a) replace those that are already non-zero with non-0 values vN0 <- isN0(value[iN0]) @@ -782,7 +783,7 @@ clDx <- getClassDef(clx <- class(x)) } - ii.v <- .Call(m_encodeInd, i - 1L, di, checkBounds = TRUE)# 0-indexing + ii.v <- .Call(m_encodeInd, i, di, orig1=TRUE, checkBounds = TRUE) if(id <- anyDuplicated(ii.v, fromLast=TRUE)) { Matrix.msg("duplicate ij-entries in 'Matrix[ ij ] <- value'; using last", .M.level = 1) @@ -794,7 +795,7 @@ value <- value[nd] } } - ii.x <- .Call(m_encodeInd2, x@i, x@j, di, FALSE) + ii.x <- .Call(m_encodeInd2, x@i, x@j, di, FALSE, FALSE) m1 <- match(ii.v, ii.x) i.repl <- !is.na(m1) # those that need to be *replaced* diff -Nru rmatrix-1.1-3/src/Csparse.c rmatrix-1.1-4/src/Csparse.c --- rmatrix-1.1-3/src/Csparse.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/Csparse.c 2014-06-14 16:22:31.000000000 +0000 @@ -631,40 +631,49 @@ * * @return a SEXP, either a (double) number or a length n-vector of diagonal entries */ -SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind) +SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm, /* ^^^^^^ FIXME[Generalize] to int / ... */ + SEXP resultKind) { const char* res_ch = CHAR(STRING_ELT(resultKind,0)); - enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log + enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range } res_kind = ((!strcmp(res_ch, "trace")) ? trace : ((!strcmp(res_ch, "sumLog")) ? sum_log : ((!strcmp(res_ch, "prod")) ? prod : - ((!strcmp(res_ch, "diag")) ? diag : - ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : - -1))))); - int i, n_x, i_from = 0; + ((!strcmp(res_ch, "min")) ? min : + ((!strcmp(res_ch, "max")) ? max : + ((!strcmp(res_ch, "range")) ? range : + ((!strcmp(res_ch, "diag")) ? diag : + ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : + -1)))))))); + int i, n_x, i_from; SEXP ans = PROTECT(allocVector(REALSXP, /* ^^^^ FIXME[Generalize] */ (res_kind == diag || - res_kind == diag_backpermuted) ? n : 1)); + res_kind == diag_backpermuted) ? n : + (res_kind == range ? 2 : 1))); double *v = REAL(ans); /* ^^^^^^ ^^^^ FIXME[Generalize] */ + i_from = (is_U ? -1 : 0); + #define for_DIAG(v_ASSIGN) \ - for(i = 0; i < n; i++, i_from += n_x) { \ + for(i = 0; i < n; i++) { \ /* looking at i-th column */ \ n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \ + if( is_U) i_from += n_x; \ v_ASSIGN; \ + if(!is_U) i_from += n_x; \ } /* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix * for uplo = "U" (makes sense with a "dtCMatrix" !), - * should use x_x[i_from + (nx - 1)] instead of x_x[i_from], - * where nx = (x_p[i+1] - x_p[i]) + * should use x_x[i_from + (n_x - 1)] instead of x_x[i_from], + * where n_x = (x_p[i+1] - x_p[i]) */ switch(res_kind) { - case trace: + case trace: // = sum v[0] = 0.; for_DIAG(v[0] += x_x[i_from]); break; @@ -679,6 +688,23 @@ for_DIAG(v[0] *= x_x[i_from]); break; + case min: + v[0] = R_PosInf; + for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]); + break; + + case max: + v[0] = R_NegInf; + for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]); + break; + + case range: + v[0] = R_PosInf; + v[1] = R_NegInf; + for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]; + if(v[1] < x_x[i_from]) v[1] = x_x[i_from]); + break; + case diag: for_DIAG(v[i] = x_x[i_from]); break; @@ -708,6 +734,7 @@ * Extract the diagonal entries from *triangular* Csparse matrix __or__ a * cholmod_sparse factor (LDL = TRUE). * + * @param obj -- now a cholmod_sparse factor or a dtCMatrix * @param pslot 'p' (column pointer) slot of Csparse matrix/factor * @param xslot 'x' (non-zero entries) slot of Csparse matrix/factor * @param perm_slot 'perm' (= permutation vector) slot of corresponding CHMfactor; @@ -716,17 +743,27 @@ * * @return a SEXP, either a (double) number or a length n-vector of diagonal entries */ -SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind) +SEXP diag_tC(SEXP obj, SEXP resultKind) { + + SEXP + pslot = GET_SLOT(obj, Matrix_pSym), + xslot = GET_SLOT(obj, Matrix_xSym); + Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) && + *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U'); int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */ - *x_p = INTEGER(pslot), - *perm = INTEGER(perm_slot); + *x_p = INTEGER(pslot), pp = -1, *perm; double *x_x = REAL(xslot); /* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/ - return diag_tC_ptr(n, x_p, x_x, perm, resultKind); + if(R_has_slot(obj, Matrix_permSym)) + perm = INTEGER(GET_SLOT(obj, Matrix_permSym)); + else perm = &pp; + + return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind); } + /** * Create a Csparse matrix object from indices and/or pointers. * diff -Nru rmatrix-1.1-3/src/Csparse.h rmatrix-1.1-4/src/Csparse.h --- rmatrix-1.1-3/src/Csparse.h 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/Csparse.h 2014-06-14 16:22:31.000000000 +0000 @@ -41,8 +41,9 @@ SEXP Rsparse_validate(SEXP x); -SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind); -SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind); +SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm, + SEXP resultKind); +SEXP diag_tC(SEXP obj, SEXP resultKind); // FIXME: these are nowhere used (are they?) SEXP create_Csparse(char* cls, int* i, int* j, int* p, int np, diff -Nru rmatrix-1.1-3/src/dgeMatrix.c rmatrix-1.1-4/src/dgeMatrix.c --- rmatrix-1.1-3/src/dgeMatrix.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/dgeMatrix.c 2014-06-14 16:22:31.000000000 +0000 @@ -157,12 +157,14 @@ double one = 1.0, zero = 0.0; Rboolean y_has_dimNames; - if (isInteger(y)) { - y = PROTECT(coerceVector(y, REALSXP)); - nprot++; + if (!isReal(y)) { + if(isInteger(y) || isLogical(y)) { + y = PROTECT(coerceVector(y, REALSXP)); + nprot++; + } + else + error(_("Argument y must be numeric, integer or logical")); } - if (!isReal(y)) - error(_("Argument y must be numeric or integer")); if(isMatrix(y)) { yDims = INTEGER(getAttrib(y, R_DimSymbol)); yDnms = getAttrib(y, R_DimNamesSymbol); diff -Nru rmatrix-1.1-3/src/dsCMatrix.c rmatrix-1.1-4/src/dsCMatrix.c --- rmatrix-1.1-3/src/dsCMatrix.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/dsCMatrix.c 2014-06-14 16:22:31.000000000 +0000 @@ -185,6 +185,7 @@ ans = PROTECT(diag_tC_ptr(L->n, L->p, L->x, + /* is_U = */ FALSE, L->Perm, resultKind)); cholmod_free_factor(&L, &c); @@ -193,9 +194,13 @@ } // using cholmod_spsolve() --> sparse result -SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b) +SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b, SEXP LDL) { - CHM_FR L = internal_chm_factor(a, /*perm*/-1, /*LDL*/-1, /*super*/-1, /*Imult*/0.); + int iLDL = asLogical(LDL); + // When parameter is set to NA in R, let CHOLMOD choose + if(iLDL == NA_LOGICAL) iLDL = -1; + + CHM_FR L = internal_chm_factor(a, /*perm*/-1, iLDL, /*super*/-1, /*Imult*/0.); CHM_SP cx, cb; if(!chm_factor_ok(L)) return R_NilValue;// == "CHOLMOD factorization failed" @@ -211,9 +216,13 @@ } // using cholmod_solve() --> dense result -SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b) +SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP LDL) { - CHM_FR L = internal_chm_factor(a, -1, -1, -1, 0.); + int iLDL = asLogical(LDL); + // When parameter is set to NA in R, let CHOLMOD choose + if(iLDL == NA_LOGICAL) iLDL = -1; + + CHM_FR L = internal_chm_factor(a, /*perm*/-1, iLDL, /*super*/-1, /*Imult*/0.); CHM_DN cx, cb; if(!chm_factor_ok(L)) return R_NilValue;// == "CHOLMOD factorization failed" diff -Nru rmatrix-1.1-3/src/dsCMatrix.h rmatrix-1.1-4/src/dsCMatrix.h --- rmatrix-1.1-3/src/dsCMatrix.h 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/dsCMatrix.h 2014-06-14 16:22:31.000000000 +0000 @@ -12,8 +12,8 @@ SEXP dsCMatrix_Cholesky(SEXP A, SEXP perm, SEXP LDL, SEXP super, SEXP Imult); SEXP dsCMatrix_LDL_D(SEXP Ap, SEXP permP, SEXP resultKind); SEXP dsCMatrix_chol(SEXP x, SEXP pivot); -SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b); -SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b); +SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b, SEXP LDL); +SEXP dsCMatrix_matrix_solve (SEXP a, SEXP b, SEXP LDL); SEXP dsCMatrix_to_dgTMatrix(SEXP x); #endif diff -Nru rmatrix-1.1-3/src/dsyMatrix.c rmatrix-1.1-4/src/dsyMatrix.c --- rmatrix-1.1-3/src/dsyMatrix.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/dsyMatrix.c 2014-06-14 16:22:31.000000000 +0000 @@ -120,12 +120,18 @@ if ((rt && n != adims[0]) || (!rt && m != adims[0])) error(_("Matrices are not conformable for multiplication")); - if (m < 1 || n < 1) { -/* error(_("Matrices with zero extents cannot be multiplied")); */ - } else + if (m >=1 && n >= 1) F77_CALL(dsymm)(rt ? "R" :"L", uplo_P(a), &m, &n, &one, REAL(GET_SLOT(a, Matrix_xSym)), adims, bcp, &m, &zero, vx, &m); + // add dimnames: + if(rt) { // v <- b %*% a : rownames(v) == rownames(b) are already there + SET_VECTOR_ELT(GET_SLOT(val, Matrix_DimNamesSym), 1, + duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); + } else { // v <- a %*% b : colnames(v) == colnames(b) are already there + SET_VECTOR_ELT(GET_SLOT(val, Matrix_DimNamesSym), 0, + duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); + } UNPROTECT(1); return val; } diff -Nru rmatrix-1.1-3/src/init.c rmatrix-1.1-4/src/init.c --- rmatrix-1.1-3/src/init.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/init.c 2014-06-14 16:22:31.000000000 +0000 @@ -73,7 +73,7 @@ CALLDEF(Csparse_vertcat, 2), CALLDEF(pCholesky_validate, 1), CALLDEF(Rsparse_validate, 1), - CALLDEF(diag_tC, 4), + CALLDEF(diag_tC, 2), #ifdef _valid_only_for_old_graph_package CALLDEF(graphNEL_as_dgTMatrix, 2), #endif @@ -152,8 +152,8 @@ CALLDEF(dsCMatrix_Cholesky, 5), CALLDEF(dsCMatrix_LDL_D, 3), CALLDEF(dsCMatrix_chol, 2), - CALLDEF(dsCMatrix_Csparse_solve, 2), - CALLDEF(dsCMatrix_matrix_solve, 2), + CALLDEF(dsCMatrix_Csparse_solve, 3), + CALLDEF(dsCMatrix_matrix_solve, 3), CALLDEF(dsCMatrix_to_dgTMatrix, 1), CALLDEF(dsTMatrix_as_dgTMatrix, 1), CALLDEF(lsTMatrix_as_lgTMatrix, 1), @@ -255,13 +255,13 @@ CALLDEF(CHM_set_common_env, 1), CALLDEF(inv_permutation, 3), - CALLDEF(m_encodeInd, 3), - CALLDEF(m_encodeInd2, 4), + CALLDEF(m_encodeInd, 4), + CALLDEF(m_encodeInd2, 5), CALLDEF(Matrix_rle_i, 2), CALLDEF(Matrix_rle_d, 2), - CALLDEF(R_set_factors, 3), + CALLDEF(R_set_factors, 4), CALLDEF(get_SuiteSparse_version, 0), {NULL, NULL, 0} diff -Nru rmatrix-1.1-3/src/Mutils.c rmatrix-1.1-4/src/Mutils.c --- rmatrix-1.1-3/src/Mutils.c 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/Mutils.c 2014-06-14 16:22:31.000000000 +0000 @@ -163,9 +163,15 @@ } // R interface [for updating the '@ factors' slot of a function *argument* [CARE!] -SEXP R_set_factors(SEXP obj, SEXP val, SEXP name) +SEXP R_set_factors(SEXP obj, SEXP val, SEXP name, SEXP warn) { - return set_factors(obj, val, (char *)CHAR(asChar(name))); + Rboolean do_warn = asLogical(warn); + if(R_has_slot(obj, Matrix_factorSym)) + return set_factors(obj, val, (char *)CHAR(asChar(name))); + else { + if(do_warn) warning(_("Matrix object has no 'factors' slot")); + return val; + } } #if 0 /* unused */ @@ -792,13 +798,14 @@ DUP_MMATRIX_NON_CLASS; - if (isInteger(A) || isLogical(A)) { - A = PROTECT(coerceVector(A, REALSXP)); - nprot++; + if (!isReal(A)) { + if (isInteger(A) || isLogical(A)) { + A = PROTECT(coerceVector(A, REALSXP)); + nprot++; + } else + error(_("invalid class '%s' to dup_mMatrix_as_dgeMatrix"), + class_P(A)); } - if (!isReal(A)) - error(_("invalid class '%s' to dup_mMatrix_as_dgeMatrix"), - class_P(A)); } DUP_MMATRIX_SET_1; @@ -834,41 +841,47 @@ * * @return encoded index; integer if prod(dim) is small; double otherwise */ -SEXP m_encodeInd(SEXP ij, SEXP di, SEXP chk_bnds) +SEXP m_encodeInd(SEXP ij, SEXP di, SEXP orig_1, SEXP chk_bnds) { SEXP ans; - int *ij_di = NULL, n, *Di = INTEGER(di), *IJ, *j_; - Rboolean check_bounds = asLogical(chk_bnds); + int *ij_di = NULL, n, nprot=1; + Rboolean check_bounds = asLogical(chk_bnds), one_ind = asLogical(orig_1); - ij = PROTECT(coerceVector(ij, INTSXP)); + if(TYPEOF(di) != INTSXP) {di = PROTECT(coerceVector(di, INTSXP)); nprot++; } + if(TYPEOF(ij) != INTSXP) {ij = PROTECT(coerceVector(ij, INTSXP)); nprot++; } if(!isMatrix(ij) || (ij_di = INTEGER(getAttrib(ij, R_DimSymbol)))[1] != 2) error(_("Argument ij must be 2-column integer matrix")); n = ij_di[0]; - IJ = INTEGER(ij); - j_ = IJ+n;/* pointer offset! */ + int *Di = INTEGER(di), *IJ = INTEGER(ij), + *j_ = IJ+n;/* pointer offset! */ if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */ ans = PROTECT(allocVector(REALSXP, n)); double *ii = REAL(ans), nr = (double) Di[0]; -#define do_ii_FILL(_i_, _j_) \ +#define do_ii_FILL(_i_, _j_) \ int i; \ if(check_bounds) { \ for(i=0; i < n; i++) { \ if(_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \ ii[i] = NA_INTEGER; \ else { \ - if(_i_[i] < 0 || _i_[i] >= Di[0]) \ + register int i_i, j_i; \ + if(one_ind) { i_i = _i_[i]-1; j_i = _j_[i]-1; } \ + else { i_i = _i_[i] ; j_i = _j_[i] ; } \ + if(i_i < 0 || i_i >= Di[0]) \ error(_("subscript 'i' out of bounds in M[ij]")); \ - if(_j_[i] < 0 || _j_[i] >= Di[1]) \ + if(j_i < 0 || j_i >= Di[1]) \ error(_("subscript 'j' out of bounds in M[ij]")); \ - ii[i] = _i_[i] + _j_[i] * nr; \ + ii[i] = i_i + j_i * nr; \ } \ } \ } else { \ for(i=0; i < n; i++) \ ii[i] = (_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \ - ? NA_INTEGER : _i_[i] + _j_[i] * nr; \ + ? NA_INTEGER \ + : (one_ind ? ((_i_[i]-1) + (_j_[i]-1)*nr) \ + : _i_[i] + _j_[i] *nr); \ } do_ii_FILL(IJ, j_); @@ -878,7 +891,7 @@ do_ii_FILL(IJ, j_); } - UNPROTECT(2); + UNPROTECT(nprot); return ans; } @@ -887,22 +900,24 @@ * * @param i: integer vector * @param j: integer vector of same length as 'i' + * @param orig_1: logical: if TRUE, "1-origin" otherwise "0-origin" * @param di: dim(.), i.e. length 2 integer vector * @param chk_bnds: logical indicating 0 <= ij[,k] < di[k] need to be checked. * * @return encoded index; integer if prod(dim) is small; double otherwise */ -SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP chk_bnds) +SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP orig_1, SEXP chk_bnds) { SEXP ans; - int n = LENGTH(i); - int *Di = INTEGER(di), *i_, *j_; - Rboolean check_bounds = asLogical(chk_bnds); + int n = LENGTH(i), nprot = 1; + Rboolean check_bounds = asLogical(chk_bnds), one_ind = asLogical(orig_1); - if(LENGTH(j) != n || !isInteger(i) || !isInteger(j)) + if(TYPEOF(di)!= INTSXP) {di = PROTECT(coerceVector(di,INTSXP)); nprot++; } + if(TYPEOF(i) != INTSXP) { i = PROTECT(coerceVector(i, INTSXP)); nprot++; } + if(TYPEOF(j) != INTSXP) { j = PROTECT(coerceVector(j, INTSXP)); nprot++; } + if(LENGTH(j) != n) error(_("i and j must be integer vectors of the same length")); - i_ = INTEGER(i); - j_ = INTEGER(j); + int *Di = INTEGER(di), *i_ = INTEGER(i), *j_ = INTEGER(j); if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */ ans = PROTECT(allocVector(REALSXP, n)); @@ -915,7 +930,7 @@ do_ii_FILL(i_, j_); } - UNPROTECT(1); + UNPROTECT(nprot); return ans; } #undef do_ii_FILL diff -Nru rmatrix-1.1-3/src/Mutils.h rmatrix-1.1-4/src/Mutils.h --- rmatrix-1.1-3/src/Mutils.h 2014-03-30 16:34:01.000000000 +0000 +++ rmatrix-1.1-4/src/Mutils.h 2014-06-14 16:22:31.000000000 +0000 @@ -72,7 +72,7 @@ SEXP as_det_obj(double val, int log, int sign); SEXP get_factors(SEXP obj, char *nm); SEXP set_factors(SEXP obj, SEXP val, char *nm); -SEXP R_set_factors(SEXP obj, SEXP val, SEXP name); +SEXP R_set_factors(SEXP obj, SEXP val, SEXP name, SEXP warn); #if 0 SEXP dgCMatrix_set_Dim(SEXP x, int nrow); @@ -287,8 +287,8 @@ SEXP dup_mMatrix_as_geMatrix (SEXP A); SEXP new_dgeMatrix(int nrow, int ncol); -SEXP m_encodeInd (SEXP ij, SEXP di, SEXP chk_bnds); -SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP chk_bnds); +SEXP m_encodeInd (SEXP ij, SEXP di, SEXP orig_1, SEXP chk_bnds); +SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP orig_1, SEXP chk_bnds); SEXP R_all0(SEXP x); SEXP R_any0(SEXP x); diff -Nru rmatrix-1.1-3/tests/dpo-test.R rmatrix-1.1-4/tests/dpo-test.R --- rmatrix-1.1-3/tests/dpo-test.R 2014-01-28 15:34:21.000000000 +0000 +++ rmatrix-1.1-4/tests/dpo-test.R 2014-05-08 10:31:52.000000000 +0000 @@ -6,8 +6,10 @@ stopifnot(c(0,0) == dim(Hilbert(0)), c(9,9) == dim(h9), identical(h9@factors, list())) -str(h9)# no 'factors' -all.equal(c(determinant(h9)$modulus), -96.7369456, tolerance = 2e-8) +str(h9)# no 'factors' 32b: -96.73694669 2.08e-8 +assert.EQ.(c(determinant(h9)$modulus), -96.7369487, tol = 8e-8) +## 64b: -96.73695078 2.15e-8 + ## determinant() now working via chol(): ==> h9 now has factorization stopifnot(names(h9@factors) == "Cholesky", identical(ch9 <- chol(h9), h9@factors$Cholesky)) diff -Nru rmatrix-1.1-3/tests/factorizing.R rmatrix-1.1-4/tests/factorizing.R --- rmatrix-1.1-3/tests/factorizing.R 2014-02-28 22:50:04.000000000 +0000 +++ rmatrix-1.1-4/tests/factorizing.R 2014-04-12 21:37:37.000000000 +0000 @@ -314,20 +314,25 @@ ## the inverse permutation: invP <- solve(P)@perm lDet <- sum(2* log(d))# the "true" value - ldet <- Matrix:::.diag.dsC(Chx = CAp, res.kind = "sumLog") + ldetp <- Matrix:::.diag.dsC(Chx = CAp, res.kind = "sumLog") + ldetp. <- sum(log(Matrix:::.diag.dsC(Chx = CAp, res.kind = "diag") )) ## CA <- Cholesky(A,perm=FALSE) - ldet2 <- Matrix:::.diag.dsC(Chx = CA, res.kind = "sumLog") + ldet <- Matrix:::.diag.dsC(Chx = CA, res.kind = "sumLog") ## not printing CAp : ends up non-integer for n >= 11 mCAp <- as(CAp,"sparseMatrix") print(mCA <- drop0(as(CA, "sparseMatrix"))) stopifnot(identical(A[p,p], as(P %*% A %*% t(P), "symmetricMatrix")), - all.equal(lDet, sum(log(Matrix:::.diag.dsC(Chx= CAp,res.kind="diag")))), relErr(d.^2, Matrix:::.diag.dsC(Chx= CA, res.kind="diag")) < 1e-14, - all.equal(lDet, ldet), - all.equal(lDet, ldet2), relErr(A[p,p], tcrossprod(mCAp)) < 1e-14) + if(FALSE) + rbind(lDet,ldet, ldetp, ldetp.) + ## ==> Empirically, I see lDet = ldet != ldetp == ldetp. + ## if(rr$rcond.A < ...) warning("condition number of A ..." ## <- TODO + cat(1,""); assert.EQ.(lDet, ldet, tol = 1e-14) + cat(2,""); assert.EQ.(ldetp, ldetp., tol = 1e-14) + cat(3,""); assert.EQ.(lDet, ldetp, tol = n^2* 1e-7)# extreme: have seen 0.0011045 !! }## for() mkCholhash <- function(r.all) { @@ -402,7 +407,7 @@ CAinv2 <- as(CAinv2, "symmetricMatrix") stopifnot(identical(CAinv, CAinv2)) -## FINALLY fix this "TODO": +## FINALLY fix "TODO": (not implemented *symbolic* factorization of nMatrix) try( tc <- Cholesky(nT.) ) for(p in c(FALSE,TRUE)) diff -Nru rmatrix-1.1-3/tests/matprod.R rmatrix-1.1-4/tests/matprod.R --- rmatrix-1.1-3/tests/matprod.R 2014-03-12 11:54:33.000000000 +0000 +++ rmatrix-1.1-4/tests/matprod.R 2014-06-14 16:21:38.000000000 +0000 @@ -3,7 +3,9 @@ ### Matrix Products including cross products source(system.file("test-tools.R", package = "Matrix")) -options(warn=1)# show as they happen +doExtras +options(warn=1, # show as they happen + Matrix.verbose = doExtras) ### dimnames -- notably for matrix products ## from ../R/Auxiliaries.R : @@ -110,11 +112,19 @@ (trm %*% trm)[i.lt]) ## crossprod() with numeric vector RHS and LHS -assert.EQ.mat( crossprod(rep(1,5), m5), rbind( colSums(m5))) -assert.EQ.mat( crossprod(rep(1,5), m.), rbind( colSums(m5))) -assert.EQ.mat( crossprod(m5, rep(1,5)), cbind( colSums(m5))) -assert.EQ.mat( crossprod(m., rep(1,5)), cbind( colSums(m5))) +i5 <- rep.int(1, 5) +isValid(S5 <- tcrossprod(tr5), "dpoMatrix")# and inherits from "dsy" +G5 <- as(S5, "generalMatrix")# "dge" +assert.EQ.mat( crossprod(i5, m5), rbind( colSums(m5))) +assert.EQ.mat( crossprod(i5, m.), rbind( colSums(m5))) +assert.EQ.mat( crossprod(m5, i5), cbind( colSums(m5))) +assert.EQ.mat( crossprod(m., i5), cbind( colSums(m5))) +assert.EQ.mat( crossprod(i5, S5), rbind( colSums(S5))) # failed in Matrix 1.1.4 ## tcrossprod() with numeric vector RHS and LHS : +stopifnot(identical(tcrossprod(i5, S5), # <- lost dimnames + tcrossprod(i5, G5) -> m51), + identical(dimnames(m51), list(NULL, LETTERS[1:5])) + ) m51 <- m5[, 1, drop=FALSE] # [6 x 1] m.1 <- m.[, 1, drop=FALSE] ; assert.EQ.mat(m51, m.1) ## The only (M . v) case @@ -578,5 +588,17 @@ } else if(n %% 25 == 0) cat(n, " ") }; cat("\n") +## two with an empty column --- these failed till 2014-06-14 +X <- as(c(1,3,4,5,3), "indMatrix") +Y <- as(c(2,3,4,2,2), "indMatrix") + +## kronecker: +stopifnot(identical(X %x% Y, + as(as.matrix(X) %x% as.matrix(Y), "indMatrix"))) +## crossprod: +(XtY <- crossprod(X, Y))# gave warning in Matrix 1.1-3 +XtY_ok <- as(crossprod(as.matrix(X), as.matrix(Y)), "dgCMatrix") +stopifnot(identical(XtY, XtY_ok)) # not true, previously + cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' diff -Nru rmatrix-1.1-3/tests/Simple.R rmatrix-1.1-4/tests/Simple.R --- rmatrix-1.1-3/tests/Simple.R 2014-03-12 18:44:33.000000000 +0000 +++ rmatrix-1.1-4/tests/Simple.R 2014-06-14 16:22:09.000000000 +0000 @@ -874,9 +874,12 @@ lsUtr <- lst[istri][uniC] (di <- sapply(lsUtr, function(.) dim(get(.)))) ## TODO: use %*%, crossprod(), .. on all those 4 x 4 -- and check "triangular rules" -cat('Time elapsed: ', (.pt <- proc.time()),'\n') # "stats" -## +r <- tryCatch(chol2inv(Diagonal(x=1:10), pi=pi), warning=identity) +stopifnot(grepl("extra argument pi .*chol2inv\\(Diagonal", r$message)) + + +cat('Time elapsed: ', (.pt <- proc.time()),'\n') # "stats" cat("doExtras:",doExtras,"\n") if(doExtras) { cat("checkMatrix() of all: \n---------\n") diff -Nru rmatrix-1.1-3/TODO rmatrix-1.1-4/TODO --- rmatrix-1.1-3/TODO 2014-03-04 07:35:12.000000000 +0000 +++ rmatrix-1.1-4/TODO 2014-05-08 10:31:52.000000000 +0000 @@ -26,7 +26,7 @@ of some samples, i.e., columns of all 0 s. It's mathematically more natural --> typically will be useful. ** DONE R/products.Rout -- [t]crossprod() could/should become more lenient with *vector*s -*** TODO consider analagous changes to base-R [t]crossprod() +*** TODO consider analagous changes to base-R [t]crossprod() ** TODO cor() and cov() at least for y=NULL ("no y"). -> ~/R/MM/Pkg-ex/Matrix/cor_sparse-propos.R <- http://stackoverflow.com/questions/5888287/ ** TODO Investigate the "band changing (and getting) ideas 'band<-' etc, @@ -39,8 +39,8 @@ ** TODO use getOption("Matrix.quiet") in more places [--> less messages/warnings] ** TODO Check for DimNames propagation in coercion and other operations. - Done for (%*%, crossprod, tcrossprod). - + Done for (%*%, crossprod, tcrossprod) -- but not systematically checked (tests/matprod.R) +*** TODO For colSums(), rowSums() ** TODO Report the problem in the Linux ldexp manual page. The second and third calls in the Synopsis should be to ldexpf and ldexpl. @@ -60,6 +60,10 @@ *** DONE .chkName.CHM(name, perm, LDL, super) and .CHM.factor.name() *** TODO use the above ** TODO partly DONE; new arg 'cache=FALSE': allow cache=FALSE to disable the caching +** TODO 0-based vs 1-based indexing: grep -nHE -e '[01]-(orig|ind|base)' *.R + Can I find a *uniform* language '1-based indexing' or '0-origin indexing' ? +*** More systemtic possible via new argumnet 'orig_1' in m_encodeInd(), m_encodeInd2() + -> src/Mutils.c * Generalization of Existing Classes and Methods --------------------------- ** TODO "Math2" , "Math", "Arith": keep triangular and symmetric Matrices when appropriate: