R/ddmatrix_qr.r

In pbdDMAT: 'pbdR' Distributed Matrix Methods

#' QR Decomposition Methods
#' 
#' \code{qr()} takes the QR decomposition.
#' 
#' \code{qr.Q()} recovers Q from the output of \code{qr()}.
#' 
#' \code{qr.R()} recovers R from the output of \code{qr()}.
#' 
#' \code{qr.qy()} multiplies \code{y} by Q.
#' 
#' \code{qr.qty()} multiplies \code{y} by the transpose of Q.
#' 
#' Functions for forming a QR decomposition and for using the outputs of these
#' numerical QR routines.
#' 
#' @param x,y 
#' numeric distributed matrices for \code{qr()}. Otherwise, \code{x}
#' is a list, namely the return from \code{qr()}.
#' @param tol 
#' logical value, determines whether or not columns are zero
#' centered.
#' @param complete 
#' logical expression of length 1.  Indicates whether an
#' arbitrary orthogonal completion of the Q or X matrices is to be made, or
#' whether the R matrix is to be completed by binding zero-value rows beneath
#' the square upper triangle.
#' @param Dvec 
#' Not implemented for objects of class \code{ddmatrix}.  vector
#' (not matrix) of diagonal values.  Each column of the returned Q will be
#' multiplied by the corresponding diagonal value.  Defaults to all 1's.
#' 
#' @return 
#' \code{qr()} returns a list consisting of: \code{qr} - \code{rank} -
#' calculated numerical rank, \code{tau} - \code{pivot} - \code{"class"} -
#' attribute "qr".
#' 
#' @keywords Methods Linear Algebra
#' @aliases qr.Q qr.R qr.qty qr.qy
#' @name qr
#' @rdname qr
NULL



setGeneric(name="qr", 
  function(x, ...)
    standardGeneric("qr"),
  package="pbdDMAT"
)

setGeneric(name="qr.Q", 
  function(x, ...)
    standardGeneric("qr.Q"),
  package="pbdDMAT"
)

setGeneric(name="qr.R", 
  function(x, ...)
    standardGeneric("qr.R"),
  package="pbdDMAT"
)

setGeneric(name="qr.qy", 
  function(x, ...)
    standardGeneric("qr.qy"),
  package="pbdDMAT"
)

setGeneric(name="qr.qty", 
  function(x, ...)
    standardGeneric("qr.qty"),
  package="pbdDMAT"
)



#' @rdname qr
#' @export
setMethod("qr", signature(x="ddmatrix"), 
  function (x, tol = 1e-07)
  {
    descx <- base.descinit(x@dim, x@bldim, x@ldim, ICTXT=x@ICTXT)
    
    m <- descx[3L]
    n <- descx[4L]
    
    ICTXT = get.comm.from.ICTXT(x@ICTXT)
    out <- base.rpdgeqpf(tol=tol, m=m, n=n, x=x@Data, descx=descx, comm=ICTXT)
    
    if (base.ownany(dim=x@dim, bldim=x@bldim, ICTXT=x@ICTXT))
      x@Data <- out$qr
    else
      x@Data <- matrix(0)
    
    ret <- list(qr=x, rank=out$rank, tau=out$tau, pivot=out$pivot)
    attr(ret, "class") <- "qr"
    
    ret
  }
)



#' @rdname qr
#' @export
setMethod("qr.Q", signature(x="ANY"), 
  function(x, complete = FALSE,  Dvec)
    {
      # x is of class qr
      
      if (is.ddmatrix(x$qr))
      {
        # complete/Dvec options
        qr <- x$qr
        
        if (qr@dim[1L] < qr@dim[2L])
          qr <- qr[, 1L:x$rank]
        
        # Matrix descriptors
        descqr <- base.descinit(qr@dim, qr@bldim, qr@ldim, ICTXT=qr@ICTXT)
        
        m <- descqr[3]
        n <- descqr[4]
        
        k <- x$rank
        
        ret <- base.rpdorgqr(m=m, n=n, k=k, qr=qr@Data, descqr=descqr, tau=x$tau)
        
        if (base.ownany(dim=qr@dim, bldim=qr@bldim, ICTXT=qr@ICTXT)){
          qr@Data <- ret
        }
        
        return( qr )
        
      }
      else
      {
        dqr <- dim(x$qr)
        cmplx <- mode(x$qr) == "complex"
        ret <- base::qr.Q(qr=x, complete=complete, Dvec=Dvec)
      }
      
      ret
  }
)



# qr.R
dmat.qr.R <- function(qr, complete=FALSE)
{
  ret <- qr$qr
  
  if (!complete){
    if (min(ret@dim)!=ret@dim[1L])
      ret <- ret[1L:min(ret@dim), ]
  }
  
  descx <- base.descinit(dim=ret@dim, bldim=ret@bldim, ldim=ret@ldim, ICTXT=ret@ICTXT)
  
  ret@Data <- base.tri2zero(x=ret@Data, descx=descx, uplo='L', diag='N')
  
  # not particularly efficient, but I don't expect this to get any real use...
  rank <- qr$rank
  n <- ret@dim[1L]
  p <- ret@dim[2L]
  mn <- min(ret@dim)
  
  if (rank < p){
    if (n>p)
      for (i in (rank+1L):mn)
        ret[i,i] <- 0
  }
  
  ret
}



#' @rdname qr
#' @export
setMethod("qr.R", signature(x="ANY"), 
  function(x, complete = FALSE) 
  {
    qr <- x
    
    if (is.ddmatrix(qr$qr)){
      ret <- dmat.qr.R(qr=qr, complete=complete)
    } else {
      ret <- base::qr.R(qr=qr, complete=complete)
    }
    
    ret
  }
)



#' @rdname qr
#' @export
setMethod("qr.qy", signature(x="ANY"), 
  function(x, y)
  {
    if (is.ddmatrix(x$qr)){
      
      qr <- x$qr
      
      # Matrix descriptors
      descqr <- base.descinit(qr@dim, qr@bldim, qr@ldim, ICTXT=qr@ICTXT)
      descc <- base.descinit(y@dim, y@bldim, y@ldim, ICTXT=y@ICTXT)
      
      m <- descqr[3L]
      n <- y@dim[2L]
      k <- x$rank
      
      out <- base.rpdormqr(side='L', trans='N', m=m, n=n, k=k, qr=qr@Data, descqr=descqr, tau=x$tau, c=y@Data, descc=descc)
      
      y@Data <- out
      
      return(y)
      
    } else {
      ret <- base::qr.qy(qr=x, y=y)
    }
    
    ret
  }
)



#' @rdname qr
#' @export
setMethod("qr.qty", signature(x="ANY"), 
  function(x, y)
  {
    if (is.ddmatrix(x$qr)){
      
      qr <- x$qr
      
      # Matrix descriptors
      descqr <- base.descinit(qr@dim, qr@bldim, qr@ldim, ICTXT=qr@ICTXT)
      descc <- base.descinit(y@dim, y@bldim, y@ldim, ICTXT=y@ICTXT)
      
      m <- descqr[3L]
      n <- y@dim[2L]
      k <- x$rank
      
      out <- base.rpdormqr(side='L', trans='T', m=m, n=n, k=k, qr=qr@Data, descqr=descqr, tau=x$tau, c=y@Data, descc=descc)
      
      y@Data <- out
      
      return(y)
      
    } else {
      ret <- base::qr.qty(qr=x, y=y)
    }
    
    ret
  }
)



#' @rdname qr
#' @export
lq <- function (x)
{
  # Matrix descriptors
  descx <- base.descinit(x@dim, x@bldim, x@ldim, ICTXT=x@ICTXT)
  
  m <- descx[3L]
  n <- descx[4L]
  
  # qr
  out <- base.rpdgelqf(m=m, n=n, x=x@Data, descx=descx)
  
  if (base.ownany(dim=x@dim, bldim=x@bldim, ICTXT=x@ICTXT))
    x@Data <- out$lq
  
  ret <- list(lq=x, tau=out$tau)
  
  attr(ret, "class") <- "lq"
  
  ret
}


#' @rdname qr
#' @export
lq.Q <- function(x, complete = FALSE)
{
  if (!is(x, "lq"))
    comm.stop("'x' must be an 'lq' object")
  if (!is.ddmatrix(x$lq))
    comm.stop("badly formed input; lq not a 'ddmatrix'")
  
  lq <- x$lq
  
  # Matrix descriptors
  desc <- base.descinit(lq@dim, lq@bldim, lq@ldim, ICTXT=lq@ICTXT)
  
  m <- desc[3]
  n <- desc[4]
  
  k <- min(m, n)
  
  ret <- base.rpdorglq(m=m, n=n, k=k, lq=lq@Data, desc=desc, tau=x$tau)
  
  if (base.ownany(dim=lq@dim, bldim=lq@bldim, ICTXT=lq@ICTXT))
    lq@Data <- ret
  
  lq
}