From 91a1b380d7065ea1bbfd87ad5a61cf31a5a9df03 Mon Sep 17 00:00:00 2001 From: Frederic Bertrand Date: Mon, 19 Apr 2021 15:11:04 +0200 Subject: [PATCH] CRAN 1.0.0 --- .Rbuildignore | 1 + .gitignore | 1 + DESCRIPTION | 6 +- NAMESPACE | 34 ++- R/bigalgebra.R | 420 ++++++++++++++++++++++++++++++++----- R/utility.R | 2 +- R/wrappers.R | 91 +++++++- R/zzz.R | 4 +- README.Rmd | 4 +- README.md | 4 +- docs/authors.html | 4 +- docs/index.html | 8 +- docs/pkgdown.yml | 2 +- docs/reference/daxpy.html | 67 ++---- docs/reference/dcopy.html | 219 +++++++++++++++++++ docs/reference/dgeev.html | 342 ++++++++++++++++++++++++++++++ docs/reference/dgemm.html | 278 ++++++++++++++++++++++++ docs/reference/dgeqrf.html | 243 +++++++++++++++++++++ docs/reference/dgesdd.html | 324 ++++++++++++++++++++++++++++ docs/reference/dpotrf.html | 241 +++++++++++++++++++++ docs/reference/dscal.html | 201 ++++++++++++++++++ docs/reference/index.html | 42 ++++ man/daxpy.Rd | 72 +++---- man/dcopy.Rd | 47 +++++ man/dgeev.Rd | 115 ++++++++++ man/dgemm.Rd | 85 ++++++++ man/dgeqrf.Rd | 54 +++++ man/dgesdd.Rd | 129 ++++++++++++ man/dpotrf.Rd | 59 ++++++ man/dscal.Rd | 34 +++ src/Makevars | 1 - src/R_init_bigalgebra.c | 6 + src/bigalgebra.cpp | 158 ++++++++++++-- src/bigalgebra.h | 18 +- 34 files changed, 3132 insertions(+), 184 deletions(-) create mode 100644 docs/reference/dcopy.html create mode 100644 docs/reference/dgeev.html create mode 100644 docs/reference/dgemm.html create mode 100644 docs/reference/dgeqrf.html create mode 100644 docs/reference/dgesdd.html create mode 100644 docs/reference/dpotrf.html create mode 100644 docs/reference/dscal.html create mode 100644 man/dcopy.Rd create mode 100644 man/dgeev.Rd create mode 100644 man/dgemm.Rd create mode 100644 man/dgeqrf.Rd create mode 100644 man/dgesdd.Rd create mode 100644 man/dpotrf.Rd create mode 100644 man/dscal.Rd diff --git a/.Rbuildignore b/.Rbuildignore index 3fe0eb1..04ecfc5 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -48,3 +48,4 @@ genlogo.R ^vignettes/.*_cache$ ^vignettes/.*_files$ +cran_submit_text.txt diff --git a/.gitignore b/.gitignore index 37afac0..de7de0d 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ src/*.o src/*.so src/*.dll +cran_submit_text.txt diff --git a/DESCRIPTION b/DESCRIPTION index 99f3a34..ed28478 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: bigalgebra Type: Package -Title: 'BLAS' Routines for Native R Matrices and 'big.matrix' Objects +Title: 'BLAS' and 'LAPACK' Routines for Native R Matrices and 'big.matrix' Objects Version: 1.0.0 -Date: 2021-04-10 +Date: 2021-04-13 Depends: bigmemory (>= 4.0.0) Imports: methods LinkingTo: bigmemory, BH, Rcpp @@ -13,7 +13,7 @@ Authors@R: c( person(given = "John W.", family= "Emerson", role = c("aut"))) Author: Frederic Bertrand [cre, ctb] (), Michael J. Kane [aut], Bryan Lewis [aut], John W. Emerson [aut] Maintainer: Frederic Bertrand -Description: Provides arithmetic functions for R matrix and 'big.matrix' objects. +Description: Provides arithmetic functions for R matrix and 'big.matrix' objects as well as functions for QR factorization, Cholesky factorization, General eigenvalue, and Singular value decomposition (SVD). A method matrix multiplication and an arithmetic method -for matrix addition, matrix difference- allows for mixed type operation -a matrix class object and a big.matrix class object- and pure type operation for two big.matrix class objects. License: LGPL-3 | Apache License 2.0 Encoding: UTF-8 Copyright: (C) 2014 Michael J. Kane, Bryan Lewis, and John W. Emerson diff --git a/NAMESPACE b/NAMESPACE index acae8c7..3393c93 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,14 +1,42 @@ import(methods) importFrom("bigmemory", "typeof") importFrom("bigmemory", "filebacked.big.matrix") +importFrom("utils", "capture.output") + #S4 methods exportMethods("%*%") exportMethods("Arith") + #Functions +# Copy a matrix +# Y := X +export("dcopy") +# Multiply by a scalar +# Y := ALPHA * Y +export("dscal") +# Add two matrices. +# Y := ALPHA * X + Y export("daxpy") - +# Matrix Multiply +# C := ALPHA * op(A) * op(B) + BETA * C +export("dgemm") +# QR factorization +# return 0 if successful, -i if ith argument has illegal value +export("dgeqrf") +# Cholesky factorization +# return 0 if successful, <0 if -i-th argument is invalid, > 0 if leading minor +# is not positive definite +export("dpotrf") +# General eigenvalue +# return 0 if successful, <0 i-th argument has illegal value, >0 QR +# algorithm failed. +# for now, VL and VR have to be matrices but they could be NULL +export("dgeev") +# Singular value decomposition (SVD) +# Returns: = 0 if successful +# < 0 if INFO = -i had an illegal value +# > 0 if DBDSDC did not converge +export("dgesdd") # #' @useDynLib bigalgebra, .registration = TRUE useDynLib(bigalgebra, .registration = TRUE) - -importFrom("utils", "capture.output") diff --git a/R/bigalgebra.R b/R/bigalgebra.R index 8bce92e..d1eea05 100644 --- a/R/bigalgebra.R +++ b/R/bigalgebra.R @@ -24,6 +24,37 @@ check_matrix = function(A, classes=c('big.matrix', 'matrix'), # Do I need a function to add a scalar to each element of a matrix? +#' @title Copy a vector. +#' +#' @description Copy double precision DX to double precision DY. +#' For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY), +#' where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is +#' defined in a similar way using INCY. +#' @param N number of elements in input vector(s) +#' @param X double precision vector with N elements +#' @param INCX storage spacing between elements of DX +#' @param Y double precision vector with N elements +#' @param INCY storage spacing between elements of DY +#' +#' @return DY copy of vector DX (unchanged if N .LE. 0) +#' @references C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. Krogh, Basic linear algebra subprograms for Fortran usage, Algorithm No. 539, Transactions on Mathematical Software 5, 3 (September 1979), pp. 308-323. +#' @export +#' +#' @examples +#' \donttest{ +#' set.seed(4669) +#' A = big.matrix(3, 2, type="double", init=1, dimnames=list(NULL, +#' c("alpha", "beta")), shared=FALSE) +#' B = big.matrix(3, 2, type="double", init=0, dimnames=list(NULL, +#' c("alpha", "beta")), shared=FALSE) +#' +#' dcopy(X=A,Y=B) +#' A[,]-B[,] +#' +#' # The big.matrix file backings will be deleted when garbage collected. +#' rm(A,B) +#' gc() +#' } # Copy a matrix # Y := X dcopy = function(N=NULL, X, INCX=1, Y, INCY=1) @@ -34,11 +65,31 @@ dcopy = function(N=NULL, X, INCX=1, Y, INCY=1) { N = as.double(nrow(X))*as.double(ncol(X)) } - .Call('dcopy_wrapper', N, X, as.double(INCX), Y, as.double(INCY), + .Call(`_dcopy_wrapper`, N, X, as.double(INCX), Y, as.double(INCY), X.is.bm, Y.is.bm) return(0) } +#' @title Scales a vector by a constant. +#' +#' @param N an integer. Number of elements in input vector(s) +#' @param ALPHA a real number. The scalar alpha +#' @param Y a big matrix to scale by ALPHA +#' @param INCY an integer. Storage spacing between elements of Y. +#' +#' @return Update Y. +#' @export +#' +#' @examples +#' set.seed(4669) +#' A = big.matrix(3, 2, type="double", init=1, dimnames=list(NULL, +#' c("alpha", "beta")), shared=FALSE) +#' dscal(ALPHA=2,Y=A) +#' A[,] +#' +#' # The big.matrix file backings will be deleted when garbage collected. +#' rm(A) +#' gc() # Multiply by a scalar # Y := ALPHA * Y dscal = function(N=NULL, ALPHA, Y, INCY=1) @@ -48,61 +99,91 @@ dscal = function(N=NULL, ALPHA, Y, INCY=1) { N = as.double(nrow(Y))*as.double(ncol(Y)) } - .Call('dscal_wrapper', N, as.double(ALPHA), Y, as.double(INCY), Y.is.bm) + .Call(`_dscal_wrapper`, N, as.double(ALPHA), Y, as.double(INCY), Y.is.bm) return(0) } -# Add two matrices. -# Y := ALPHA * X + Y -daxpy = function(N=NULL, ALPHA=1, X, INCX=1, Y, INCY=1) -{ - X.is.bm = check_matrix(X) - Y.is.bm = check_matrix(Y) - if (is.null(N)) - { - N = as.double(nrow(X))*as.double(ncol(X)) - } - .Call('daxpy_wrapper', as.double(N), as.double(ALPHA), X, as.double(INCX), - Y, as.double(INCY), X.is.bm, Y.is.bm) - return(0) -} +# # Add two matrices. +# # Y := ALPHA * X + Y +# daxpy = function(N=NULL, ALPHA=1, X, INCX=1, Y, INCY=1) +# { +# X.is.bm = check_matrix(X) +# Y.is.bm = check_matrix(Y) +# if (is.null(N)) +# { +# N = as.double(nrow(X))*as.double(ncol(X)) +# } +# .Call("_daxpy_wrapper", as.double(N), as.double(ALPHA), X, as.double(INCX), +# Y, as.double(INCY), X.is.bm, Y.is.bm) +# return(0) +# } -# Matrix Multiply -# C := ALPHA * op(A) * op(B) + BETA * C -# This is function provides dgemm functionality. The matrix types -# are handled on the C++ side. -dgemm = function(TRANSA='n', TRANSB='n', M=NULL, N=NULL, K=NULL, - ALPHA=1, A, LDA=NULL, B, LDB=NULL, BETA=0, C, LDC=NULL, COFF=0) -{ - A.is.bm = check_matrix(A) - B.is.bm = check_matrix(B) - C.is.bm = check_matrix(C) - # The matrices look OK. Now, if they haven't been specified, let's - # specify some reasonable dimension information. - if ( is.null(M) ) - { - M = ifelse ( is_transposed(TRANSA), ncol(A), nrow(A) ) - } - if ( is.null(N) ) - { - N = ifelse ( is_transposed(TRANSB), nrow(B), ncol(B) ) - } - if ( is.null(K) ) - { - K = ifelse ( is_transposed(TRANSA), nrow(A), ncol(A) ) - } - if ( is.null(LDA) ) LDA = nrow (A) - if ( is.null(LDB) ) LDB = nrow (B) - if ( is.null(LDC) ) LDC = nrow (C) - - .Call('dgemm_wrapper', as.character(TRANSA), as.character(TRANSB), - as.double(M), as.double(N), as.double(K), as.double(ALPHA), A, - as.double(LDA), B, as.double(LDB), - as.double(BETA), C, as.double(LDC), as.logical(A.is.bm), - as.logical(B.is.bm), as.logical(C.is.bm), COFF) - return(invisible(C)) -} +# # Matrix Multiply +# # C := ALPHA * op(A) * op(B) + BETA * C +# # This is function provides dgemm functionality. The matrix types +# # are handled on the C++ side. +# dgemm = function(TRANSA='n', TRANSB='n', M=NULL, N=NULL, K=NULL, +# ALPHA=1, A, LDA=NULL, B, LDB=NULL, BETA=0, C, LDC=NULL, COFF=0) +# { +# A.is.bm = check_matrix(A) +# B.is.bm = check_matrix(B) +# C.is.bm = check_matrix(C) +# # The matrices look OK. Now, if they haven't been specified, let's +# # specify some reasonable dimension information. +# if ( is.null(M) ) +# { +# M = ifelse ( is_transposed(TRANSA), ncol(A), nrow(A) ) +# } +# if ( is.null(N) ) +# { +# N = ifelse ( is_transposed(TRANSB), nrow(B), ncol(B) ) +# } +# if ( is.null(K) ) +# { +# K = ifelse ( is_transposed(TRANSA), nrow(A), ncol(A) ) +# } +# if ( is.null(LDA) ) LDA = nrow (A) +# if ( is.null(LDB) ) LDB = nrow (B) +# if ( is.null(LDC) ) LDC = nrow (C) +# +# .Call('dgemm_wrapper', as.character(TRANSA), as.character(TRANSB), +# as.double(M), as.double(N), as.double(K), as.double(ALPHA), A, +# as.double(LDA), B, as.double(LDB), +# as.double(BETA), C, as.double(LDC), as.logical(A.is.bm), +# as.logical(B.is.bm), as.logical(C.is.bm), COFF) +# return(invisible(C)) +# } +#' @title QR factorization +#' +#' @description DGEQRF computes a QR factorization of a real M-by-N matrix A: A = Q * R. +#' @param M an integer. The number of rows of the matrix A. M >= 0. +#' @param N an integer. The number of columns of the matrix A. N >= 0. +#' @param A the M-by-N big matrix A. +#' @param LDA an integer. The leading dimension of the array A. LDA >= max(1,M). +#' @param TAU a min(M,N) matrix. The scalar factors of the elementary reflectors. +#' @param WORK a (MAX(1,LWORK)) matrix. On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +#' @param LWORK an integer. The dimension of th array WORK. +#' +#' @return M-by-N big matrix A. The elements on and above the diagonal of the +#' array contain the min(M,N)-by-N upper trapezoidal matrix R (R is upper +#' triangular if m >= n); the elements below the diagonal, with the array TAU, +#' represent the orthogonal matrix Q as a product of min(m,n) elementary +#' reflectors. +#' @export +#' +#' @examples +#' hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } +#' h9 <- hilbert(9); h9 +#' qr(h9)$rank #--> only 7 +#' qrh9 <- qr(h9, tol = 1e-10) +#' qrh9$rank +#' C <- as.big.matrix(h9) +#' dgeqrf(A=C) +#' +#' # The big.matrix file backings will be deleted when garbage collected. +#' rm(C) +#' gc() # QR factorization # return 0 if successful, -i if ith argument has illegal value dgeqrf=function(M=NULL, N=NULL, A, LDA=NULL, TAU=NULL, WORK=NULL, @@ -136,13 +217,61 @@ dgeqrf=function(M=NULL, N=NULL, A, LDA=NULL, TAU=NULL, WORK=NULL, TAU.is.bm = check_matrix(TAU) WORK.is.bm = check_matrix(WORK) INFO = 0 - .Call('dgeqrf_wrapper', as.double(M), as.double(N), A, as.double(LDA), + .Call(`_dgeqrf_wrapper`, as.double(M), as.double(N), A, as.double(LDA), TAU, WORK, as.double(LWORK), as.double(INFO), A.is.bm, TAU.is.bm, WORK.is.bm) return(INFO) } -# Cholesky factorization + + + +#' @title Cholesky factorization +#' +#' @description DPOTRF computes the Cholesky factorization of a real symmetric positive definite matrix A. +#' +#' The factorization has the form +#' \describe{ +#' \item{}{A = U**T * U, if UPLO = 'U', or} +#' \item{}{A = L * L**T, if UPLO = 'L',} +#' } +#' where U is an upper triangular matrix and L is lower triangular. +#' +#' This is the block version of the algorithm, calling Level 3 BLAS. +#' @param UPLO a character. +#' \describe{ +#' \item{}{'U': Upper triangle of A is stored;} +#' \item{}{'L': Lower triangle of A is stored.} +#' } +#' @param N an integer. The order of the matrix A. N >= 0. +#' @param A a big.matrix, dimension (LDA,N). +#' @param LDA an integer. Dimension of the array A. LDA >= max(1,N). +#' +#' @return updates the big matrix A with the result, INFO is an integer +#' \describe{ +#' \item{}{= 0: successful exit} +#' \item{}{< 0: if INFO = -i, the i-th argument had an illegal value} +#' \item{}{> 0: if INFO = i, the leading minor of order i is not positive definite, and the factorization could not be completed.} +#' } +#' Terms laying out of the computed triangle should be discarded. +#' @export +#' +#' @examples +#' set.seed(4669) +#' A = matrix(rnorm(16),4) +#' B = as.big.matrix(A %*% t(A)) +#' C = A %*% t(A) +#' chol(C) +#' dpotrf(UPLO='U', N=4, A=B, LDA=4) +#' D <- B[,] +#' D[lower.tri(D)]<-0 +#' D +#' D-chol(C) +#' t(D)%*%D-C +#' +#' #' # The big.matrix file backings will be deleted when garbage collected. +#' rm(A,B,C,D) +#' gc() # return 0 if successful, <0 if -i-th argument is invalid, > 0 if leading minor # is not positive definite dpotrf=function(UPLO='U', N=NULL, A, LDA=NULL) @@ -157,11 +286,92 @@ dpotrf=function(UPLO='U', N=NULL, A, LDA=NULL) } A.is.bm = check_matrix(A) INFO = 0 - .Call('dpotrf_wrapper', as.character(UPLO), as.double(N), A, as.double(LDA), + .Call(`_dpotrf_wrapper`, as.character(UPLO), as.double(N), A, as.double(LDA), as.double(INFO), A.is.bm) return(INFO) } +#' @title DGEEV computes eigenvalues and eigenvectors. +#' +#' @description DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices. +#' +#' DGEEV computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. +#' The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. +#' The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate-transpose of u(j). +#' +#' The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. +#' +#' @param JOBVL a character. +#' \describe{ +#' \item{}{= 'N': left eigenvectors of A are not computed;} +#' \item{}{= 'V': left eigenvectors of A are computed.} +#' } +#' @param JOBVR a character. +#' \describe{ +#' \item{}{= 'N': right eigenvectors of A are not computed;} +#' \item{}{= 'V': right eigenvectors of A are computed.} +#' } +#' @param N an integer. The order of the matrix A. N >= 0. +#' @param A a matrix of dimension (LDA,N), the N-by-N matrix A. +#' @param LDA an integer. The leading dimension of the matrix A. LDA >= max(1,N). +#' @param WR a vector of dimension (N). WR contain the real part of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. +#' @param WI a vector of dimension (N). WI contain the imaginary part of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. +#' @param VL a matrx of dimension (LDVL,N) +#' \describe{ +#' \item{}{If JOBVL = 'V', the left eigenvectors u(j) are stored one +#' after another in the columns of VL, in the same order +#' as their eigenvalues.} +#' \item{}{If JOBVL = 'N', VL is not referenced.} +#' \item{}{If the j-th eigenvalue is real, then u(j) = VL(:,j), +#' the j-th column of VL.} +#' \item{}{If the j-th and (j+1)-st eigenvalues form a complex +#' conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and +#' u(j+1) = VL(:,j) - i*VL(:,j+1).} +#' } +#' @param LDVL an integer. The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. +#' @param VR a matrix of dimension (LDVR,N). +#' \describe{ +#' \item{}{If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues.} +#' \item{}{If JOBVR = 'N', VR is not referenced.} +#' \item{}{If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR.} +#' \item{}{If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and v(j+1) = VR(:,j) - i*VR(:,j+1).} +#' } +#' @param LDVR an integer. The leading dimension of the array VR. LDVR >= 1; if JOBVR = 'V', LDVR >= N. +#' @param WORK a matrix of dimension (MAX(1,LWORK)) +#' @param LWORK an integer. The dimension of the array WORK.LWORK >= max(1,3*N), and if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good performance, LWORK must generally be larger. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. +#' +#' @return WR, WI, VR, VL and Work. On exit, A has been overwritten. +#' @export +#' +#' @examples +#' \donttest{ +#' set.seed(4669) +#' A = matrix(rnorm(16),4) +#' WR= matrix(0,nrow=4,ncol=1) +#' WI= matrix(0,nrow=4,ncol=1) +#' VL = matrix(0,ncol=4,nrow=4) +#' eigen(A) +#' dgeev(A=A,WR=WR,WI=WI,VL=VL) +#' VL +#' WR +#' WI +#' +#' rm(A,WR,WI,VL) +#' +#' A = as.big.matrix(matrix(rnorm(16),4)) +#' WR= matrix(0,nrow=4,ncol=1) +#' WI= matrix(0,nrow=4,ncol=1) +#' VL = as.big.matrix(matrix(0,ncol=4,nrow=4)) +#' eigen(A[,]) +#' dgeev(A=A,WR=WR,WI=WI,VL=VL) +#' VL[,] +#' WR[,] +#' WI[,] +#' +# The big.matrix file backings will be deleted when garbage collected. +#' rm(A,WR,WI,VL) +#' gc() +#' } # General eigenvalue # return 0 if successful, <0 i-th argument has illegal value, >0 QR # algorithm failed. @@ -211,12 +421,110 @@ dgeev=function(JOBVL='V', JOBVR='V', N=NULL, A, LDA=NULL, WR, WI, VL, VR.is.bm = check_matrix(VR) WORK.is.bm = check_matrix(WORK) INFO=0 - .Call('dgeev_wrapper', as.character(JOBVL), as.character(JOBVR), as.double(N), A, as.double(LDA), WR, WI, VL, as.double(LDVL), VR, as.double(LDVR), + .Call(`_dgeev_wrapper`, as.character(JOBVL), as.character(JOBVR), as.double(N), A, as.double(LDA), WR, WI, VL, as.double(LDVL), VR, as.double(LDVR), WORK, as.double(LWORK), as.double(INFO), A.is.bm, WR.is.bm, WI.is.bm, VL.is.bm, VR.is.bm, WORK.is.bm) return(INFO) } + +#' @title DGESDD computes the singular value decomposition (SVD) of a real matrix. +#' +#' @description DGESDD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and right singular vectors. If singular vectors are desired, it uses a divide-and-conquer algorithm. +#' +#' The SVD is written +#' +#' A = U * SIGMA * transpose(V) +#' +#' where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. +#' +#' Note that the routine returns VT = V**T, not V. +#' +#' @param JOBZ a character. Specifies options for computing all or part of the matrix U: +#' \describe{ +#' \item{}{= 'A': all M columns of U and all N rows of V**T are returned in the arrays U and VT;} +#' \item{}{= 'S': the first min(M,N) columns of U and the first min(M,N) rows of V**T are returned in the arrays U and VT;} +#' \item{}{= 'O': If M >= N, the first N columns of U are overwritten on the array A and all rows of V**T are returned in the array VT; otherwise, all columns of U are returned in the array U and the first M rows of V**T are overwritten in the array A;} +#' \item{}{= 'N': no columns of U or rows of V**T are computed.} +#' } +#' @param M an integer. The number of rows of the input matrix A. M >= 0. +#' @param N an integer. The number of columns of the input matrix A. N >= 0. +#' @param A the M-by-N matrix A. +#' @param LDA an integer. The leading dimension of the matrix A. LDA >= max(1,M). +#' @param S a matrix of dimension (min(M,N)). The singular values of A, sorted so that S(i) >= S(i+1). +#' @param U U is a matrx of dimension (LDU,UCOL) +#' \describe{ +#' \item{}{UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; UCOL = min(M,N) if JOBZ = 'S'.} +#' \item{}{If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M orthogonal matrix U;} +#' \item{}{if JOBZ = 'S', U contains the first min(M,N) columns of U (the left singular vectors, stored columnwise);} +#' \item{}{if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.} +#' } +#' @param LDU an integer. The leading dimension of the matrix U. LDU >= 1; if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. +#' @param VT VT is matrix of dimension (LDVT,N) +#' \describe{ +#' \item{}{If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the N-by-N orthogonal matrix V**T;} +#' \item{}{if JOBZ = 'S', VT contains the first min(M,N) rows of V**T (the right singular vectors, stored rowwise);} +#' \item{}{if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.} +#' } +#' @param LDVT an integer. The leading dimension of the matrix VT. LDVT >= 1; if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; if JOBZ = 'S', LDVT >= min(M,N). +#' @param WORK a matrix of dimension (MAX(1,LWORK)) +#' @param LWORK an integer. The dimension of the array WORK. LWORK >= 1. +#' If LWORK = -1, a workspace query is assumed. The optimal +#' size for the WORK array is calculated and stored in WORK(1), +#' and no other work except argument checking is performed. +#' +#' Let mx = max(M,N) and mn = min(M,N). +#' \describe{ +#' \item{}{If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).} +#' \item{}{If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).} +#' \item{}{If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.} +#' \item{}{If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.} +#' } +#' These are not tight minimums in all cases; see comments inside code. +#' For good performance, LWORK should generally be larger; +#' a query is recommended. +#' +#' @return IWORK an integer matrix dimension of (8*min(M,N)) +#' A is updated. +#' \describe{ +#' \item{}{if JOBZ = 'O', A is overwritten with the first N columns of U (the left singular vectors, stored columnwise) if M >= N; A is overwritten with the first M rows of V**T (the right singular vectors, stored rowwise) otherwise.} +#' \item{}{if JOBZ .ne. 'O', the contents of A are destroyed.} +#' } +#' INFO an integer +#' \describe{ +#' \item{}{= 0: successful exit.} +#' \item{}{< 0: if INFO = -i, the i-th argument had an illegal value.} +#' \item{}{> 0: DBDSDC did not converge, updating process failed.} +#' } +#' @export +#' +#' @examples +#' \donttest{ +#' set.seed(4669) +#' A = matrix(rnorm(12),4,3) +#' S = matrix(0,nrow=3,ncol=1) +#' U = matrix(0,nrow=4,ncol=4) +#' VT = matrix(0,ncol=3,nrow=3) +#' dgesdd(A=A,S=S,U=U,VT=VT) +#' S +#' U +#' VT +#' +#' rm(A,S,U,VT) +#' +#' A = as.big.matrix(matrix(rnorm(12),4,3)) +#' S = as.big.matrix(matrix(0,nrow=3,ncol=1)) +#' U = as.big.matrix(matrix(0,nrow=4,ncol=4)) +#' VT = as.big.matrix(matrix(0,ncol=3,nrow=3)) +#' dgesdd(A=A,S=S,U=U,VT=VT) +#' S[,] +#' U[,] +#' VT[,] +#' +#' rm(A,S,U,VT) +#' gc() +#' } +# Singular value decomposition (SVD) # Returns: = 0 if successful # < 0 if INFO = -i had an illegal value # > 0 if DBDSDC did not converge @@ -279,7 +587,7 @@ dgesdd = function( JOBZ='A', M=NULL, N=NULL, A, LDA=NULL, S, U, LDU=NULL, } WORK.is.bm = check_matrix(WORK) INFO = 0 - .Call('dgesdd_wrapper', as.character(JOBZ), as.double(M), as.double(N), A, + .Call(`_dgesdd_wrapper`, as.character(JOBZ), as.double(M), as.double(N), A, as.double(LDA), S, U, as.double(LDU), VT, as.double(LDVT), WORK, as.double(LWORK), as.double(INFO), A.is.bm, S.is.bm, U.is.bm, VT.is.bm, WORK.is.bm) diff --git a/R/utility.R b/R/utility.R index 3980d8a..cfdbe88 100644 --- a/R/utility.R +++ b/R/utility.R @@ -42,7 +42,7 @@ anon_matrix = function(m, n, type, val=NULL) descriptorfile=d) address = capture.output(print(ans@address)) path = paste(p,f,sep="//") - if(bigdebug()) warning(paste("Creating anonymous maitrx ",path)) + if(bigdebug()) warning(paste("Creating anonymous matrix ",path)) assign(address, path, envir=.bigalgebra_env) reg.finalizer(ans@address, finalize_anon_matrix, onexit=TRUE) ans diff --git a/R/wrappers.R b/R/wrappers.R index a5d967a..6296a59 100644 --- a/R/wrappers.R +++ b/R/wrappers.R @@ -12,6 +12,54 @@ # +#' Matrix Multiply +#' +#' @description This is function provides dgemm functionality, which DGEMM +#' performs one of the matrix-matrix operations. +#' C := ALPHA * op(A) * op(B) + BETA * C. +#' +#' @param TRANSA a character. TRANSA specifies the form of op( A ) to be used in the matrix multiplication as follows: +#' \describe{ +#' \item{}{TRANSA = 'N' or 'n', op( A ) = A.} +#' \item{}{TRANSA = 'T' or 't', op( A ) = A**T.} +#' \item{}{TRANSA = 'C' or 'c', op( A ) = A**T.} +#' } +#' @param TRANSB a character. TRANSB specifies the form of op( B ) to be used in the matrix multiplication as follows: +#' #' \describe{ +#' \item{}{TRANSA = 'N' or 'n', op( B ) = B.} +#' \item{}{TRANSA = 'T' or 't', op( B ) = B**T.} +#' \item{}{TRANSA = 'C' or 'c', op( B ) = B**T.} +#' } +#' @param M an integer. M specifies the number of rows of the matrix op( A ) and of the matrix C. M must be at least zero. +#' @param N an integer. N specifies the number of columns of the matrix op( B ) and of the matrix C. N must be at least zero. +#' @param K an integer. K specifies the number of columns of the matrix op( A ) and the number of rows of the matrix op( B ). K must be at least zero. +#' @param ALPHA a real number. Specifies the scalar alpha. +#' @param A a matrix of dimension (LDA, ka), where ka is k when TRANSA = 'N' or 'n', and is m otherwise. Before entry with TRANSA = 'N' or 'n', the leading m by k part of the array A must contain the matrix A, otherwise the leading k by m part of the array A must contain the matrix A. +#' @param LDA an integer. +#' @param B a matrix of dimension ( LDB, kb ), where kb is n when TRANSB = 'N' or 'n', and is k otherwise. Before entry with TRANSB = 'N' or 'n', the leading k by n part of the array B must contain the matrix B, otherwise the leading n by k part of the array B must contain the matrix B. +#' @param LDB an integer. +#' @param BETA a real number. Specifies the scalar beta +#' @param C a matrix of dimension ( LDC, N ). Before entry, the leading m by n part of the array C must contain the matrix C, except when beta is zero, in which case C need not be set on entry. On exit, the array C is overwritten by the m by n matrix ( alpha*op( A )*op( B ) + beta*C ). +#' @param LDC an integer. +#' @param COFF offset for C. +#' +#' @return Update C with the result. +#' @export +#' +#' @examples +#' require(bigmemory) +#' A = as.big.matrix(matrix(1, nrow=3, ncol=2)) +#' B <- big.matrix(2, 3, type="double", init=-1, +#' dimnames=list(NULL, c("alpha", "beta")), shared=FALSE) +#' C = big.matrix(3, 3, type="double", init=1, +#' dimnames=list(NULL, c("alpha", "beta", "gamma")), shared=FALSE) +#' 2*A[,]%*%B[,]+0.5*C[,] +#' E = dgemm(ALPHA=2.0, A=A, B=B, BETA=0.5, C=C) +#' E[,] # Same result +#' +#' # The big.matrix file backings will be deleted when garbage collected. +#' rm(A,B,C,E) +#' gc() # Matrix Multiply # C := ALPHA * op(A) * op(B) + BETA * C # This is function provides dgemm functionality. @@ -49,12 +97,53 @@ dgemm = function(TRANSA='N', TRANSB='N', M=NULL, N=NULL, K=NULL, as.logical(B.is.bm), as.logical(C.is.bm), COFF) } +#' @title BLAS daxpy functionality +#' +#' @description This function implements the function Y := A * X + Y where X and Y may be either native double-precision valued R matrices or numeric vectors, or double-precision valued \code{\link[bigmemory]{big.matrix}} objects, and A is a scalar. +#' @param A Optional numeric scalar value to scale the matrix \code{X} by, with a default value of 1. +#' @param X Requried to be either a native \R \code{\link{matrix}} or numeric vector, or a \code{\link[bigmemory]{big.matrix}} object +#' @param Y Optional native \R \code{\link{matrix}} or numeric vector, or a \code{\link[bigmemory]{big.matrix}} object +#' +#' @details At least one of either \code{X} or \code{Y} must be a \code{big.matrix}. All values must be of type \code{double} (the only type presently supported by the bigalgebra package). +#' +#' This function is rarely necessary to use directly since the bigalgebra package defines standard arithmetic operations and scalar multiplication. It is more efficient to use \code{daxpy} directly when both scaling and matrix addition are required, in which case both operations are performed in one step. +#' +#' @return The output value depends on the classes of input values \code{X} and \code{Y} and on the value of the global option \code{bigalgebra.mixed_arithmetic_returns_R_matrix}. +#' +#' If \code{X} and \code{Y} are both big matrices, or \code{Y} is missing, \code{options("bigalgebra.mixed_arithmetic_returns_R_matrix")} is \code{FALSE}, then a \code{big.matrix} is returned. The returned \code{big.matrix} is backed by a temporary file mapping that will be deleted when the returned result is garbage collected by R (see the examples). +#' +#' Otherwise, a standard R matrix is returned. The dimensional shape of the output is taken from \code{X}. If input \code{X} is dimensionless (that is, lacks a dimension attribute), then the output is a column vector. +#' +#' @references \url{https://www.netlib.org/blas/daxpy.f} +#' @author Michael J. Kane +#' @seealso \code{\link[bigmemory]{bigmemory}} +#' @export +#' +#' @examples +#' require(bigmemory) +#'A = matrix(1, nrow=3, ncol=2) +#'B <- big.matrix(3, 2, type="double", init=0, +#' dimnames=list(NULL, c("alpha", "beta")), shared=FALSE) +#'C = B + B # C is a new big matrix +#'D = A + B # D defaults to a regular R matrix, to change this, set the option: +#'# options(bigalgebra.mixed_arithmetic_returns_R_matrix=FALSE) +#'E = daxpy(A=1.0, X=B, Y=B) # Same kind of result as C +#'print(C[]) +#'print(D) +#'print(E[]) +#' +#'# The C and E big.matrix file backings will be deleted when garbage collected: +#'# (We enable debugging to see this explicitly) +#'options(bigalgebra.DEBUG=TRUE) +#'rm(C,E) +#'gc() +#' # Vector addition and scaling # Y := A * X + Y # A is a scalar double # X is either a big.matrix or regular R matrix. # Y is an optional matrix or vector of the same length as X. -# Returns a new matrix or big matrx with the same dimensions as X. If +# Returns a new matrix or big matrix with the same dimensions as X. If # X is a dimension-less R vector, returns a column. Returned value type # depends on the arguments and the value of the option # options("bigalgebra.mixed_airthmetic_returns_R_matrix")[[1]]. diff --git a/R/zzz.R b/R/zzz.R index b5eb6bf..22fde8f 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -3,7 +3,7 @@ # Packgage library and global options .onLoad <- function(libname, pkgname) { - library.dynam("bigalgebra", pkgname, libname); +# library.dynam("bigalgebra", pkgname, libname); options(bigalgebra.temp_pattern="matrix_") options(bigalgebra.mixed_arithmetic_returns_R_matrix=TRUE) options(bigalgebra.tempdir=tempdir) @@ -11,7 +11,7 @@ } .onUnload <- function(libpath) { - library.dynam.unload("bigalgebra", libpath); +# library.dynam.unload("bigalgebra", libpath); options(bigalgebra.temp_pattern=c()) options(bigalgebra.mixed_arithmetic_returns_R_matrix=c()) options(bigalgebra.tempdir=c()) diff --git a/README.Rmd b/README.Rmd index 683d601..640f262 100644 --- a/README.Rmd +++ b/README.Rmd @@ -28,8 +28,8 @@ knitr::opts_chunk$set( [![DOI](https://zenodo.org/badge/136206211.svg)](https://zenodo.org/badge/latestdoi/136206211) -This package provides arithmetic functions for native `R` matrices and -`bigmemory::big.matrix` objects. +This package provides arithmetic functions for native `R` matrices and `bigmemory::big.matrix` objects as well as functions for QR factorization, Cholesky factorization, General eigenvalue, and Singular value decomposition (SVD). A method matrix multiplication and an arithmetic method -for matrix addition, matrix difference- allows for mixed type operation -a matrix class object and a big.matrix class object- and pure type operation for two big.matrix class objects. + The package defines a number of global options that begin with `bigalgebra`. diff --git a/README.md b/README.md index fad4b4e..bf8bc5a 100644 --- a/README.md +++ b/README.md @@ -19,8 +19,8 @@ [![DOI](https://zenodo.org/badge/136206211.svg)](https://zenodo.org/badge/latestdoi/136206211) -This package provides arithmetic functions for native `R` matrices and -`bigmemory::big.matrix` objects. +This package provides arithmetic functions for native `R` matrices and `bigmemory::big.matrix` objects as well as functions for QR factorization, Cholesky factorization, General eigenvalue, and Singular value decomposition (SVD). A method matrix multiplication and an arithmetic method -for matrix addition, matrix difference- allows for mixed type operation -a matrix class object and a big.matrix class object- and pure type operation for two big.matrix class objects. + The package defines a number of global options that begin with `bigalgebra`. diff --git a/docs/authors.html b/docs/authors.html index 7488045..a289968 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -122,9 +122,9 @@

Citation

Source: inst/CITATION -

Frederic Bertrand, Michael J. Kane, John Emerson and Stephen Weston(2021). 'BLAS' Routines for Native R Matrices and 'big.matrix' Objects, R package version 1.0.0.

+

Frederic Bertrand, Michael J. Kane, John Emerson and Stephen Weston(2021). 'BLAS' and 'LAPACK' Routines for Native R Matrices and 'big.matrix' Objects, R package version 1.0.0.

@Manual{,
-  title = {'BLAS' Routines for Native R Matrices and 'big.matrix' Objects},
+  title = {'BLAS' and 'LAPACK' Routines for Native R Matrices and 'big.matrix' Objects},
   author = {F. Bertrand and Michael J. Kane and John Emerson and Stephen Weston},
   publisher = {manual},
   year = {2021},
diff --git a/docs/index.html b/docs/index.html
index 4c10f68..2ada50f 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -5,7 +5,7 @@
 
 
 
-BLAS Routines for Native R Matrices and big.matrix Objects • bigalgebra
+BLAS and LAPACK Routines for Native R Matrices and big.matrix Objects • bigalgebra
 
 
 
@@ -17,8 +17,8 @@
 
 
 
-
-
+
+
 
 
 
-
This package provides arithmetic functions for native R matrices and bigmemory::big.matrix objects.

+
This package provides arithmetic functions for native R matrices and bigmemory::big.matrix objects as well as functions for QR factorization, Cholesky factorization, General eigenvalue, and Singular value decomposition (SVD). A method matrix multiplication and an arithmetic method -for matrix addition, matrix difference- allows for mixed type operation -a matrix class object and a big.matrix class object- and pure type operation for two big.matrix class objects.

 
The package defines a number of global options that begin with bigalgebra.

 
They include:

 
Option Default value

diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index 1cbcc26..b74dd0c 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -2,5 +2,5 @@ pandoc: 2.11.2
 pkgdown: 1.6.1
 pkgdown_sha: ~
 articles: {}
-last_built: 2021-04-10T13:56Z
+last_built: 2021-04-16T10:42Z
 
diff --git a/docs/reference/daxpy.html b/docs/reference/daxpy.html
index c8da6ca..9bc5ff2 100644
--- a/docs/reference/daxpy.html
+++ b/docs/reference/daxpy.html
@@ -47,10 +47,7 @@
 
 
 
-
+
 
 
 
@@ -123,81 +120,60 @@
   

     

     
BLAS daxpy functionality

-    
+    Source: R/wrappers.R
     
     

 
     

-    
This function implements the function Y := A * X + Y where X and Y may be
-either native double-precision valued R matrices or numeric vectors, or
-double-precision valued big.matrix objects, and A is a
-scalar.

+    
This function implements the function Y := A * X + Y where X and Y may be either native double-precision valued R matrices or numeric vectors, or double-precision valued big.matrix objects, and A is a scalar.

     

 
-    
daxpy(A=1, X, Y)

+    
daxpy(A = 1, X, Y)

 
     
Arguments

     
-      
+      
-      
+      
-      
+      

     
     

       A
Optional numeric scalar value to scale the matrix X by,
-           with a default value of 1.
Optional numeric scalar value to scale the matrix X by, with a default value of 1.

     

     

       X
Requried to be either a native R matrix
-           or numeric vector,
-           or a big.matrix object
Requried to be either a native R matrix or numeric vector, or a big.matrix object

     

     

       Y
Optional native R matrix or numeric vector,
-           or a big.matrix object
Optional native R matrix or numeric vector, or a big.matrix object

     

     

 
-    
Details

-
-    
At least one of either X or Y must be a big.matrix.
-All values must be of type double (the only type presently supported
-by the bigalgebra package).

-
This function is rarely necessary to use directly since the bigalgebra
-package defines standard arithmetic operations and scalar multiplication.
-It is more efficient to use daxpy directly when both scaling and
-matrix addition are required, in which case both operations are performed
-in one step.

     
Value

 
-    
The output value depends on the classes of input values X and Y
-and on the value of the global option
-bigalgebra.mixed_arithmetic_returns_R_matrix.

-
If X and Y are both big matrices, or Y is missing,
-options("bigalgebra.mixed_arithmetic_returns_R_matrix") is FALSE,
-then a big.matrix is returned. The returned big.matrix is backed
-by a temporary file mapping that will be deleted when the returned result is
-garbage collected by R (see the examples).

-
Otherwise, a standard R matrix is returned. The dimensional shape of the output
-is taken from X. If input X is dimensionless (that is, lacks a
-dimension attribute), then the output is a column vector.

+    
The output value depends on the classes of input values X and Y and on the value of the global option bigalgebra.mixed_arithmetic_returns_R_matrix.

+
If X and Y are both big matrices, or Y is missing, options("bigalgebra.mixed_arithmetic_returns_R_matrix") is FALSE, then a big.matrix is returned. The returned big.matrix is backed by a temporary file mapping that will be deleted when the returned result is garbage collected by R (see the examples).

+
Otherwise, a standard R matrix is returned. The dimensional shape of the output is taken from X. If input X is dimensionless (that is, lacks a dimension attribute), then the output is a column vector.

+    
Details

+
+    
At least one of either X or Y must be a big.matrix. All values must be of type double (the only type presently supported by the bigalgebra package).

+
This function is rarely necessary to use directly since the bigalgebra package defines standard arithmetic operations and scalar multiplication. It is more efficient to use daxpy directly when both scaling and matrix addition are required, in which case both operations are performed in one step.

     
References

 
     
https://www.netlib.org/blas/daxpy.f

-    
Author

-
-    
Michael J. Kane

     
See also

 
     
bigmemory

+    
Author

+
+    
Michael J. Kane

 
     
Examples

     
require(bigmemory)
 A = matrix(1, nrow=3, ncol=2)
 B <- big.matrix(3, 2, type="double", init=0,
-                dimnames=list(NULL, c("alpha", "beta")), shared=FALSE)
+               dimnames=list(NULL, c("alpha", "beta")), shared=FALSE)
 C = B + B   # C is a new big matrix
 D = A + B   # D defaults to a regular R matrix, to change this, set the option:
-            # optons(bigalgebra.mixed_arithmetic_returns_R_matrix=FALSE)
+# options(bigalgebra.mixed_arithmetic_returns_R_matrix=FALSE)
 E = daxpy(A=1.0, X=B, Y=B)  # Same kind of result as C
 print(C[])
 
#>      [,1] [,2]
@@ -218,8 +194,9 @@ 
Examp
 rm(C,E)
 gc()
 
#>           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
-#> Ncells  761688 40.7    1430568 76.5         NA  1430568 76.5
-#> Vcells 1421587 10.9    8388608 64.0      65536  2561285 19.6

+#> Ncells  765022 40.9    1434878 76.7         NA  1434878 76.7
+#> Vcells 1427548 10.9    8388608 64.0      65536  2523363 19.3

+