R/kronecker.R

Defines functions .kroneckerDimnames

## METHODS FOR GENERIC: kronecker
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.kroneckerDimnames <- function(dnX, dX = lengths(dnX, FALSE),
                               dnY, dY = lengths(dnX, FALSE),
                               sep = ":") {
    dnr <- list(NULL, NULL)
    if(identical(dnX, dnr) && identical(dnY, dnr))
        return(NULL)
    for(i in 1:2) {
        if(have.sX <- (nY <- dY[i]) && !is.null(sX <- dnX[[i]]))
            sX <- rep    (sX,  each = nY)
        if(have.sY <- (nX <- dX[i]) && !is.null(sY <- dnY[[i]]))
            sY <- rep.int(sY, times = nX)
        dnr[[i]] <-
            if(have.sX && have.sY)
                paste0(sX, sep, sY)
            else if(have.sX)
                paste0(sX, sep    )
            else if(have.sY)
                paste0(    sep, sY)
            else if(nX && nY)
                rep.int(sep, nX * nY)
    }
    dnr
}

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              r <- new("ddiMatrix")
              r@Dim <- dX * dY
              if((uX <- X@diag != "N") & (uY <- Y@diag != "N"))
                  r@diag <- "U"
              else if(uX)
                  r@x <- rep.int(as.double(Y@x), dX[1L])
              else if(uY)
                  r@x <- rep(as.double(X@x), each = dY[1L])
              else r@x <- rep(X@x, each = dY[1L]) * Y@x
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY)))
                  r@Dimnames <- dnr
              r
	  })

if(FALSE) { # --NOT YET--
setMethod("kronecker", signature(X = "diagonalMatrix", Y = "denseMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- X@diag != "N"
              uY <- FALSE
              shape <- .M.shape(Y)
              r <- new(`substr<-`("d.CMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  uplo <- r@uplo <- Y@uplo
                  if(shape == "t")
                      uY <- Y@diag != "N"
              }
              if(!all(dr)) {
                  r@p <- integer(dr[2L] + 1)
                  if(uX && uY)
                      r@diag <- "U"
              } else {
                  m <- dY[1L]
                  nX <- dX[1L]
                  nY <- length(y <- Y@x)
                  if(shape != "g" && nY > 1L && nY == prod(dY)) {
                      Y <- pack(Y)
                      nY <- length(y <- Y@x)
                  }
                  if(as.double(nX) * nY > .Machine$integer.max)
                      stop("number of nonzero entries cannot exceed 2^31-1")
                  if(!uX && uY) {
                      diag(Y) <- TRUE
                      nY <- length(y <- Y@x)
                  }
                  if(is.logical(y) && .M.kind(Y) == "n")
                      y <- y | is.na(y)
                  if(shape == "g") {
                      r@p <- seq.int(0L, by = m, length.out = dr[2L] + 1)
                      r@i <-
                          rep(seq.int(0L, by = m, length.out = nX),
                              each = nY) +
                          0:(m-1L)
                  } else if(uplo == "U") {
                      r@p <- c(0L, cumsum(rep.int(1:m, nX)))
                      r@i <-
                          rep(seq.int(0L, by = m, length.out = nX),
                              each = nY) +
                          sequence.default(nvec = 1:m, from = 0L)
                  } else {
                      r@p <- c(0L, cumsum(rep.int(m:1, nX)))
                      r@i <-
                          rep(seq.int(0L, by = m, length.out = nX),
                              each = nY) +
                          sequence.default(nvec = m:1, from = 0:(m-1L))
                  }
                  r@x <-
                      if(uX)
                          rep.int(as.double(y), nX)
                      else as.double(y) * rep(X@x, each = nY)
                  if(uX && uY)
                      r <- ..diagN2U(r, sparse = TRUE)
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })

setMethod("kronecker", signature(X = "denseMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              shape <- .M.shape(X)
              uX <- FALSE
              uY <- Y@diag != "N"
              r <- new(`substr<-`("d.CMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  uplo <- r@uplo <- X@uplo
                  if(shape == "t")
                      uX <- X@diag != "N"
              }
              if(!all(dr)) {
                  r@p <- integer(dr[2L] + 1)
                  if(uX && uY)
                      r@diag <- "U"
              } else {
                  m <- dX[1L]
                  n <- dX[2L]
                  nX <- length(x <- X@x)
                  nY <- dY[1L]
                  if(shape != "g" && nX > 1L && nX == prod(dX)) {
                      X <- pack(X)
                      nX <- length(x <- X@x)
                  }
                  if(as.double(nX) * nY > .Machine$integer.max)
                      stop("number of nonzero entries cannot exceed 2^31-1")
                  if(uX && !uY) {
                      diag(X) <- TRUE
                      nX <- length(x <- X@x)
                  }
                  if(is.logical(x) && .M.kind(X) == "n")
                      x <- x | is.na(x)
                  if(shape == "g") {
                      x. <- function() {
                          j <- rep(1:n, each = nY)
                          as.vector(`dim<-`(as.double(x), dX)[, j])
                      }
                      r@p <- seq.int(0L, by = m, length.out = dr[2L] + 1)
                      r@i <- rep.int(sequence.default(nvec = rep.int(m, nY),
                                                      from = 0:(nY-1L),
                                                      by = nY),
                                     n)
                      r@x <-
                          if(uY)
                              x.()
                          else x.() * rep(Y@x, each = m)
                  } else if(uplo == "U") {
                      rep.1.n <- rep(1:n, each = nY)
                      s <- sequence.default(
                          nvec = rep.1.n,
                          from = rep(1L + cumsum(0:(n-1L)), each = nY),
                          by = 1L)

                      r@p <- c(0L, cumsum(rep.1.n))
                      r@i <- sequence.default(nvec = rep.1.n,
                                              from = rep.int(0:(nY-1L), n),
                                              by = nY)
                      r@x <-
                          if(uY)
                              as.double(x)[s]
                          else
                              as.double(x)[s] *
                                  rep.int(rep.int(Y@x, n), rep.1.n)
                  } else {
                      rep.n.1 <- rep(n:1, each = nY)
                      s <- sequence.default(
                          nvec = rep.n.1,
                          from = rep(1L + cumsum(c(0L, if(n > 1L) n:2)),
                                     each = nY),
                          by = 1L)
                      r@p <- c(0L, cumsum(rep.n.1))
                      r@i <- sequence.default(nvec = rep.n.1,
                                              from = 0:(n*nY-1L),
                                              by = nY)
                      r@x <-
                          if(uY)
                              as.double(x)[s]
                          else
                              as.double(x)[s] *
                                  rep.int(rep.int(Y@x, n), rep.n.1)
                  }
                  if(uX && uY)
                      r <- ..diagN2U(r, sparse = TRUE)
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })
} # --NOT YET--

setMethod("kronecker", signature(X = "denseMatrix", Y = "denseMatrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              shape <- switch(.M.shape(X),
                              g = "g",
                              t = if(.M.shape(Y) == "t" && X@uplo == Y@uplo)
                                      "t"
                                  else "g",
                              s = if(.M.shape(Y) == "s")
                                      "s"
                                  else "g")
              r <- new(switch(shape,
                              g = "dgeMatrix",
                              t = "dtrMatrix",
                              s = "dsyMatrix"))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  r@uplo <- X@uplo
                  if(shape == "t" && X@diag != "N" && Y@diag != "N")
                      r@diag <- "U"
              }
              if(all(dr)) {
                  X <- .dense2m(X)
                  Y <- .dense2m(Y)
                  storage.mode(X) <- "double"
                  storage.mode(Y) <- "double"
                  r@x <-
                      as.vector(X[rep(seq_len(dX[1L]), each = dY[1L]),
                                  rep(seq_len(dX[2L]), each = dY[2L])]) *
                      as.vector(Y[rep.int(seq_len(dY[1L]), dX[1L]), ])
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
          })

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "CsparseMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- X@diag != "N"
              uY <- FALSE
              shape <- .M.shape(Y)
              r <- new(`substr<-`("d.CMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  r@uplo <- Y@uplo
                  if(shape == "t" && (uX & (uY <- Y@diag != "N")))
                      r@diag <- "U"
              }
              if(!all(dr))
                  r@p <- integer(dr[2L] + 1)
              else {
                  if(!uX && uY)
                      Y <- ..diagU2N(Y)
                  if((nY <- (p <- Y@p)[length(p)]) == 0L)
                      r@p <- integer(dr[2L] + 1)
                  else if(as.double(nX <- dX[1L]) * nY > .Machine$integer.max)
                      stop("number of nonzero entries cannot exceed 2^31-1")
                  else {
                      head. <-
                          if(length(Y@i) > nY)
                              function(x) x[seq_len(nY)]
                          else identity
                      r@p <- c(0L, cumsum(rep.int(p[-1L] - p[-length(p)], nX)))
                      r@i <-
                          rep(seq.int(0L, by = dY[1L], length.out = nX),
                              each = nY) +
                          head.(Y@i)
                      r@x <-
                          if(uX) {
                              if(.M.kind(Y) != "n")
                                  rep.int(as.double(head.(Y@x)), nX)
                              else rep.int(1, nX * nY)
                          } else {
                              if(.M.kind(Y) != "n")
                                  rep(as.double(X@x), each = nY) *
                                      as.double(head.(Y@x))
                              else rep(as.double(X@x), each = nY)
                          }
                  }
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })

setMethod("kronecker", signature(X = "CsparseMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- FALSE
              uY <- Y@diag != "N"
              shape <- .M.shape(X)
              r <- new(`substr<-`("d.CMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  r@uplo <- X@uplo
                  if(shape == "t" && (uX <- X@diag != "N") && uY)
                      r@diag <- "U"
              }
              if(!all(dr))
                  r@p <- integer(dr[2L] + 1)
              else {
                  if(uX && !uY)
                      X <- ..diagU2N(X)
                  if((nX <- (p <- X@p)[length(p)]) == 0L)
                      r@p <- integer(dr[2L] + 1)
                  else if(as.double(nY <- dY[1L]) * nX > .Machine$integer.max)
                      stop("number of nonzero entries cannot exceed 2^31-1")
                  else {
                      dp <- p[-1L] - p[-length(p)]
                      j. <- which(dp > 0L)
                      nj. <- length(j.)

                      rep.dp <- rep(dp[j.], each = nY)
                      s. <- sequence.default(
                          nvec = rep.dp,
                          from = rep(1L + p[j.], each = nY),
                          by = 1L)

                      r@p <- c(0L, cumsum(rep(dp, each = nY)))
                      r@i <- (nY * X@i)[s.] +
                          rep.int(rep.int(0:(nY-1L), nj.), rep.dp)
                      r@x <-
                          if(uY) {
                              if(.M.kind(X) != "n")
                                 as.double(X@x)[s.]
                              else
                                 rep.int(1, nX * nY)
                          } else {
                              if(.M.kind(X) != "n")
                                 rep.int(rep.int(as.double(Y@x), nj.), rep.dp) *
                                     as.double(X@x)[s.]
                              else
                                 rep.int(rep.int(as.double(Y@x), nj.), rep.dp)
                          }
                  }
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })

setMethod("kronecker", signature(X = "CsparseMatrix", Y = "CsparseMatrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- uY <- FALSE
              if((sX <- .M.shape(X)) == "t")
                  uX <- X@diag != "N"
              if((sY <- .M.shape(Y)) == "t")
                  uY <- Y@diag != "N"
              shape <-
                  switch(sX,
                         g = "g",
                         t = if(sY == "t" && X@uplo == Y@uplo) "t" else "g",
                         s = if(sY == "s") "s" else "g")
              cl <- "d.CMatrix"
              if(!all(dr <- dX * dY)) {
                  substr(cl, 2L, 2L) <- shape
                  r <- new(cl)
                  r@Dim <- dr
                  r@p <- integer(dr[2L] + 1)
                  if(shape != "g") {
                      r@uplo <- X@uplo
                      if(shape == "t" && uX && uY)
                          r@diag <- "U"
                  }
              } else {
                  substr(cl, 2L, 2L) <- if(shape == "s") "g" else shape
                  r <- new(cl)
                  r@Dim <- dr
                  if(uX)
                      X <- ..diagU2N(X)
                  if(uY)
                      Y <- ..diagU2N(Y)
                  if(sY == "s")
                      Y <- .sparse2g(Y)
                  else if(sX == "s")
                      X <- .sparse2g(X)
                  if((nX <- (pX <- X@p)[length(pX)]) == 0L ||
                     (nY <- (pY <- Y@p)[length(pY)]) == 0L)
                      r@p <- integer(dr[2L] + 1)
                  else if(as.double(nX) * nY > .Machine$integer.max)
                      stop("number of nonzero entries cannot exceed 2^31-1")
                  else {
                      dpX <- pX[-1L] - (pX. <- pX[-length(pX)])
                      dpY <- pY[-1L] - (pY. <- pY[-length(pY)])

                      rep.dpX <- rep(dpX, each = dY[2L])
                      rep.dpY <- rep.int(dpY, dX[2L])

                      t1 <- rep.int(rep.dpY, rep.dpX)

                      s1 <- sequence.default(
                          nvec = rep.dpX,
                          from = rep(1L + pX., each = dY[2L]),
                          by = 1L)

                      s2 <- sequence.default(
                          nvec = t1,
                          from = rep.int(rep.int(1L + pY., dX[2L]), rep.dpX),
                          by = 1L)

                      r@p <- c(0L, cumsum(rep.dpX * dpY))
                      r@i <- rep.int((dY[1L] * X@i)[s1], t1) + Y@i[s2]
                      r@x <-
                          if(.M.kind(X) != "n") {
                              if(.M.kind(Y) != "n")
                                  rep.int(as.double(X@x)[s1], t1) *
                                      as.double(Y@x)[s2]
                              else
                                  rep.int(as.double(X@x)[s1], t1)
                          } else {
                              if(.M.kind(Y) != "n")
                                  as.double(Y@x)[s2]
                              else
                                  rep.int(1, nX * nY)
                          }
                  }
                  if(shape == "t") {
                      r@uplo <- X@uplo
                      if(uX && uY)
                          r <- ..diagN2U(r, sparse = TRUE)
                  } else if(shape == "s")
                      r <- .Call(R_sparse_force_symmetric, r, X@uplo)
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
          })

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "RsparseMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              .tCR2RC(kronecker(t(X), .tCR2RC(Y), FUN, make.dimnames, ...)))

setMethod("kronecker", signature(X = "RsparseMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              .tCR2RC(kronecker(.tCR2RC(X), t(Y), FUN, make.dimnames, ...)))

setMethod("kronecker", signature(X = "RsparseMatrix", Y = "RsparseMatrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              .tCR2RC(kronecker(.tCR2RC(X), .tCR2RC(Y), FUN, make.dimnames, ...)))

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "TsparseMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- X@diag != "N"
              uY <- FALSE
              shape <- .M.shape(Y)
              r <- new(`substr<-`("d.TMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  r@uplo <- Y@uplo
                  if(shape == "t" && (uX & (uY <- Y@diag != "N")))
                      r@diag <- "U"
              }
              if(all(dr)) {
                  if(any((kind <- .M.kind(Y)) == c("n", "l")))
                      Y <- .Call(Tsparse_aggregate, Y)
                  if(!uX && uY)
                      Y <- ..diagU2N(Y)
                  nX <- dX[1L]
                  nY <- length(Y@i)
                  r@i <-
                      rep(seq.int(0L, by = dY[1L], length.out = nX),
                          each = nY) +
                      Y@i
                  r@j <-
                      rep(seq.int(0L, by = dY[2L], length.out = nX),
                          each = nY) +
                      Y@j
                  r@x <-
                      if(uX) {
                          if(kind != "n")
                              rep.int(as.double(Y@x), nX)
                          else rep.int(1, as.double(nX) * nY)
                      } else {
                          if(kind != "n")
                              rep(as.double(X@x), each = nY) * as.double(Y@x)
                          else rep(as.double(X@x), each = nY)
                      }
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })

setMethod("kronecker", signature(X = "TsparseMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- FALSE
              uY <- Y@diag != "N"
              shape <- .M.shape(X)
              r <- new(`substr<-`("d.TMatrix", 2L, 2L, shape))
              r@Dim <- dr <- dX * dY
              if(shape != "g") {
                  r@uplo <- X@uplo
                  if(shape == "t" && (uX <- X@diag != "N") && uY)
                      r@diag <- "U"
              }
              if(all(dr)) {
                  if(any((kind <- .M.kind(X)) == c("n", "l")))
                      X <- .Call(Tsparse_aggregate, X)
                  if(uX && !uY)
                      X <- ..diagU2N(X)
                  nX <- length(X@i)
                  nY <- dY[1L]
                  r@i <- rep(nY * X@i, each = nY) + 0:(nY-1L)
                  r@j <- rep(nY * X@j, each = nY) + 0:(nY-1L)
                  r@x <-
                      if(uY) {
                          if(kind != "n")
                              rep.int(as.double(Y@x), nY)
                          else rep.int(1, as.double(nX) * nY)
                      } else {
                          if(kind != "n")
                              rep(as.double(X@x), each = nY) * as.double(Y@x)
                          else rep(as.double(X@x), each = nY)
                      }
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
	  })

setMethod("kronecker", signature(X = "TsparseMatrix", Y = "TsparseMatrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              uX <- uY <- FALSE
              if((sX <- .M.shape(X)) == "t")
                  uX <- X@diag != "N"
              if((sY <- .M.shape(Y)) == "t")
                  uY <- Y@diag != "N"
              shape <-
                  switch(sX,
                         g = "g",
                         t = if(sY == "t" && X@uplo == Y@uplo) "t" else "g",
                         s = if(sY == "s") "s" else "g")
              cl <- "d.TMatrix"
              if(!all(dr <- dX * dY)) {
                  substr(cl, 2L, 2L) <- shape
                  r <- new(cl)
                  r@Dim <- dr
                  if(shape != "g") {
                      r@uplo <- X@uplo
                      if(shape == "t" && uX && uY)
                          r@diag <- "U"
                  }
              } else {
                  substr(cl, 2L, 2L) <- if(shape == "s") "g" else shape
                  r <- new(cl)
                  r@Dim <- dr
                  if(uX)
                      X <- ..diagU2N(X)
                  if(uY)
                      Y <- ..diagU2N(Y)
                  if(sY == "s")
                      Y <- .sparse2g(Y)
                  else if(sX == "s")
                      X <- .sparse2g(X)
                  nY <- length(Y@i)
                  r@i <- i. <- rep(dY[1L] * X@i, each = nY) + Y@i
                  r@j <-       rep(dY[2L] * X@j, each = nY) + Y@j
                  r@x <-
                      if(.M.kind(X) != "n") {
                          if(.M.kind(Y) != "n")
                              rep(as.double(X@x), each = nY) * as.double(Y@x)
                          else
                              rep(as.double(X@x), each = nY)
                      } else {
                          if(.M.kind(Y) != "n")
                              rep.int(as.double(X@x), nY)
                          else
                              rep.int(1, length(i.))
                      }

                  if(shape == "t") {
                      r@uplo <- X@uplo
                      if(uX && uY)
                          r <- ..diagN2U(r, sparse = TRUE)
                  } else if(shape == "s")
                      r <- .Call(R_sparse_force_symmetric, r, X@uplo)
              }
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY))) {
                  if(shape == "s" && !isSymmetricDN(dnr))
                      r <- .sparse2g(r)
                  r@Dimnames <- dnr
              }
              r
          })

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "indMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, as(Y, "nsparseMatrix"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "indMatrix", Y = "diagonalMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(as(X, "nsparseMatrix"), Y, FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "indMatrix", Y = "indMatrix"),
	  function(X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(!(missing(FUN) || identical(FUN, "*")))
                  stop("method for kronecker() must use default FUN=\"*\"")
              if(any(as.double(dX <- X@Dim) * (dY <- Y@Dim) >
                     .Machine$integer.max))
		  stop("dimensions cannot exceed 2^31-1")
              r <- new("indMatrix")
              r@Dim <- dX * dY
              r@perm <- dY[2L] * rep(X@perm - 1L, each = dY[1L]) +
                  rep.int(Y@perm, dX[1L])
              if(make.dimnames &&
                 !is.null(dnr <- .kroneckerDimnames(dimnames(X), dX,
                                                    dimnames(Y), dY)))
                  r@Dimnames <- dnr
              r
	  })

## Catch everything else with these:

setMethod("kronecker", signature(X = "Matrix", Y = "matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, .m2ge(Y, "d"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "Matrix", Y = "vector"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, .m2ge(Y, "d"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "matrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(.m2ge(X, "d"), Y, FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "vector", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(.m2ge(X, "d"), Y, FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "denseMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(as(X, "CsparseMatrix"), Y, FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "CsparseMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, as(Y, "CsparseMatrix"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "RsparseMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, as(Y, "RsparseMatrix"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "TsparseMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, as(Y, "TsparseMatrix"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "diagonalMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(X, as(Y, "CsparseMatrix"), FUN, make.dimnames, ...))

setMethod("kronecker", signature(X = "indMatrix", Y = "Matrix"),
          function(X, Y, FUN = "*", make.dimnames = FALSE, ...)
              kronecker(as(X, "CsparseMatrix"), Y, FUN, make.dimnames, ...))


## MJ: no longer needed ... replacement above
if(FALSE) {
setMethod("kronecker", signature(X = "Matrix", Y = "ANY",
				 FUN = "ANY", make.dimnames = "ANY"),
	  function(X, Y, FUN, make.dimnames, ...) {
	      if(is(X, "sparseMatrix"))
		  warning("using slow kronecker() method")
	      X <- as(X, "matrix") ; Matrix(callGeneric()) })

setMethod("kronecker", signature(X = "ANY", Y = "Matrix",
				 FUN = "ANY", make.dimnames = "ANY"),
	  function(X, Y, FUN, make.dimnames, ...) {
	      if(is(Y, "sparseMatrix"))
		  warning("using slow kronecker() method")
	      Y <- as(Y, "matrix") ; Matrix(callGeneric()) })

tmp <- function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
    kronecker(as(X, "TsparseMatrix"), Y,
	      FUN = FUN, make.dimnames = make.dimnames, ...)
}
setMethod("kronecker", signature(X="diagonalMatrix", Y="ANY"          ), tmp)
setMethod("kronecker", signature(X="diagonalMatrix", Y="Matrix"       ), tmp)
setMethod("kronecker", signature(X="ANY",            Y="sparseMatrix" ), tmp)
## the above could recurse infinitely :
setMethod("kronecker", signature(X="sparseMatrix",   Y="TsparseMatrix"), tmp)

tmp <- function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
    kronecker(X, as(Y, "TsparseMatrix"),
	      FUN = FUN, make.dimnames = make.dimnames, ...)
}
setMethod("kronecker", signature(X="ANY",           Y="diagonalMatrix"), tmp)
setMethod("kronecker", signature(X="Matrix",        Y="diagonalMatrix"), tmp)
setMethod("kronecker", signature(X="sparseMatrix",  Y="ANY"           ), tmp)
setMethod("kronecker", signature(X="TsparseMatrix", Y="sparseMatrix"  ), tmp)
rm(tmp)

## from ./dgTMatrix.R :
setMethod("kronecker", signature(X = "dgTMatrix", Y = "dgTMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...)
      {
	  if (FUN != "*") stop("kronecker method must use default 'FUN'")
	  ## otherwise we don't know that many results will be zero
	  ydim <- Y@Dim
	  xi <- X@i
	  xnnz <- length(xi)
	  yi <- Y@i
	  ynnz <- length(yi)
	  new("dgTMatrix", Dim = X@Dim * ydim,
	      i = rep.int(yi, xnnz) + ydim[1] * rep.int(xi, rep.int(ynnz, xnnz)),
	      j = rep.int(Y@j, xnnz) + ydim[2] * rep.int(X@j, rep.int(ynnz, xnnz)),
	      ## faster than x = as.vector(outer(Y@x, X@x, FUN = FUN)
	      x = as.vector(Y@x %*% t(X@x)))
      })

## triangularity -- should be preserved "when obvious":
setMethod("kronecker", signature(X = "dtTMatrix", Y = "dtTMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...)
      {
	  if (FUN != "*") stop("kronecker method must use default 'FUN'")
	  ## otherwise we don't know that many results will be zero
	  if(X@uplo != Y@uplo) { ## result not triangular
	      X <- .sparse2g(X)
	      Y <- .sparse2g(Y)
	      return(callGeneric())
	  }
	  ## else: both 'uplo' are the same -- result *is* triangular
	  ## d.U <- (dX <- X@diag == "U") && (dY <- Y@diag == "U")
	  if(Y@diag == "U")
	      Y <- .diagU2N(Y, "dtTMatrix")
	  ydim <- Y@Dim
	  if(X@diag != "U") {
	      xi <- X@i
	      xj <- X@j
	      xx <- X@x
	  } else { ## X@diag == "U"
	      nx <- X@Dim[1] # triangular matrices are square
	      ii <- seq_len(nx) - 1L
	      xi <- c(X@i, ii)
	      xj <- c(X@j, ii)
	      xx <- c(X@x, rep.int(1, nx))
	  }
	  xnnz <- length(xi)
	  yi <- Y@i
	  ynnz <- length(yi)
	  new("dtTMatrix", Dim = X@Dim * ydim,
	      i = rep.int(yi,  xnnz) + ydim[1] * rep.int(xi, rep.int(ynnz, xnnz)),
	      j = rep.int(Y@j, xnnz) + ydim[2] * rep.int(xj, rep.int(ynnz, xnnz)),
	      ## faster than x = as.vector(outer(Y@x, X@x, FUN = FUN)
	      x = as.vector(Y@x %*% t(xx)),
	      uplo = X@uplo,
	      diag = "N" # if(d.U) { "U" , but drop the entries}  else "N"
	      )
      })

setMethod("kronecker", signature(X = "dtTMatrix", Y = "dgTMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(it <- isTriangular(Y))
                  ## improve: also test for unit diagonal
                  Y <- if(attr(it, "kind") == "U") triu(Y) else tril(Y)
	      else
                  X <- .sparse2g(X)
	      callGeneric() #-> dtT o dtT   or	 dgT o dgT
	  })

setMethod("kronecker", signature(X = "dgTMatrix", Y = "dtTMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(it <- isTriangular(X))
                  ## improve: also test for unit diagonal
                  X <- if(attr(it, "kind") == "U") triu(X) else tril(X)
	      else
		  Y <- .sparse2g(Y)
	      callGeneric() #-> dtT o dtT   or	 dgT o dgT
	  })

setMethod("kronecker", signature(X = "TsparseMatrix", Y = "TsparseMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(.hasSlot(X, "uplo") && !.hasSlot(X, "diag"))
                  X <- .sparse2g(X)
              if(.hasSlot(Y, "uplo") && !.hasSlot(Y, "diag"))
                  Y <- .sparse2g(Y)
              X <- ..sparse2d(X)
              Y <- ..sparse2d(Y)
              callGeneric()
	  })

setMethod("kronecker", signature(X = "dsparseMatrix", Y = "dsparseMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
              if(.hasSlot(X, "uplo") && !.hasSlot(X, "diag"))
                  X <- .sparse2g(X)
              if(.hasSlot(Y, "uplo") && !.hasSlot(Y, "diag"))
                  Y <- .sparse2g(Y)
              if(.hasSlot(X, "p"))
                  X <- .CR2T(X)
              if(.hasSlot(Y, "p"))
                  Y <- .CR2T(Y)
              callGeneric()
	  })
} ## MJ

Try the Matrix package in your browser

Any scripts or data that you put into this service are public.

Matrix documentation built on Nov. 11, 2022, 9:06 a.m.