SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Diff of /pkg/Matrix/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2197, Mon Jun 2 14:34:42 2008 UTC revision 2206, Wed Jun 25 14:57:02 2008 UTC
# Line 52  Line 52 
52  ### Bert's code built on a post by Andy Liaw who most probably was influenced  ### Bert's code built on a post by Andy Liaw who most probably was influenced
53  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
54  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
55    if(FALSE)##--- currently unused:
56  bdiag <- function(...) {  .bdiag <- function(lst) {
57      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
58      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
59      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
60      ## make sure we had all matrices:      ## make sure we had all matrices:
61      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
62          stop("some arguments are not matrices")          stop("some arguments are not matrices")
63      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
64                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
65      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
66        r@Dim <- as.integer(csdim[nrow(csdim),])
67      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
68      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
69          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
70          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
71              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
72          else ## square matrix          else ## square matrix
73              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
74        }
75        r
76    }
77    ### Doug Bates needed something like bdiag() for lower-triangular
78    ### (Tsparse) Matrices and provided a much more efficient implementation:
79    .bdiag <- function(lst) {
80        ### block-diagonal matrix [a dgTMatrix] from list of matrices
81        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
82    
83        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
84                       as, "TsparseMatrix")
85    
86        if(nl == 1) return(Tlst[[1]])
87        ## else
88        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
89        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
90        new("dgTMatrix", Dim = c(i_off[nl+1], j_off[nl + 1]),
91            i = unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k])),
92            j = unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])),
93            x = unlist(lapply(Tlst, slot, "x")))
94      }      }
95      ## slightly debatable if we really should return Csparse.. :  
96      as(ret, "CsparseMatrix")  bdiag <- function(...) {
97        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
98        if(nA == 1 && !is.list(...))
99            return(as(..., "CsparseMatrix"))
100        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
101        if(length(alis) == 1)
102            return(as(alis[[1]], "CsparseMatrix"))
103    
104        ## else : two or more arguments
105        as(.bdiag(alis), "CsparseMatrix")
106  }  }
107    
108    
# Line 619  Line 648 
648  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
649  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
650            function(e1,e2) {            function(e1,e2) {
651                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
652                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
653                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
654                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 666 
666    
667  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
668            function(e1,e2) {            function(e1,e2) {
669                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
670                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
671                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
672                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 685 
685  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
686  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
687            function(e1,e2) {            function(e1,e2) {
688                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
689                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
690                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
691                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 704 
704    
705  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
706            function(e1,e2) {            function(e1,e2) {
707                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
708                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
709                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
710                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 724 
724  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
725  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
726            function(e1,e2) {            function(e1,e2) {
727                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
728                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
729                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
730                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 743 
743    
744  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
745            function(e1,e2) {            function(e1,e2) {
746                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
747                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
748                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
749                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 763 
763  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
764  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
765            function(e1,e2) {            function(e1,e2) {
766                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
767                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
768                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
769                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 781 
781    
782  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
783            function(e1,e2) {            function(e1,e2) {
784                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
785                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
786                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
787                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 851  Line 880 
880            function(object) {            function(object) {
881                d <- dim(object)                d <- dim(object)
882                cl <- class(object)                cl <- class(object)
883                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
884                            d[1], d[2], cl))                            d[1], d[2], cl))
885                  if(d[1] < 50) {
886                      cat("\n")
887                prDiag(object)                prDiag(object)
888                  } else {
889                      cat(", with diagonal entries\n")
890                      show(diag(object))
891                      invisible(object)
892                  }
893            })            })

Legend:
Removed from v.2197  
changed lines
  Added in v.2206

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge