Type:	Package
Title:	'BLAS' and 'LAPACK' Routines for Native R Matrices and 'big.matrix' Objects
Version:	3.0.0
Date:	2025-10-05
Depends:	bigmemory (≥ 4.0.0)
Imports:	methods
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0)
LinkingTo:	bigmemory, BH, Rcpp
Author:	Frederic Bertrand [cre, aut], Michael J. Kane [aut], Bryan Lewis [aut], John W. Emerson [aut]
Maintainer:	Frederic Bertrand <frederic.bertrand@lecnam.net>
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
LazyLoad:	yes
NeedsCompilation:	yes
RoxygenNote:	7.3.3
URL:	https://fbertran.github.io/bigalgebra/, https://github.com/fbertran/bigalgebra/
BugReports:	https://github.com/fbertran/bigalgebra/issues/
VignetteBuilder:	knitr
Config/testthat/edition:	3
Packaged:	2025-10-05 14:30:23 UTC; bertran7
Repository:	CRAN
Date/Publication:	2025-10-05 15:10:02 UTC

Arithmetic routines for native R matrices and big.matrix objects.

Description

This package provides arithmetic functions for native R matrices and big.matrix objects.

Details

This package provides arithmetic functions for native R matrices and big.matrix objects.

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

The bigalgebra.tempdir option must be a function that returns a temporary directory path used to big matrix results of BLAS and LAPACK operations. The deault value is simply the default R tempdir function.

The bigalgebra.temp_pattern is a name prefix for file names of generated big matrix objects output as a result of BLAS and LAPACK operations.

The bigalgebra.mixed_arithmetic_returns_R_matrix option determines whether arithmetic operations involving an R matrix or vector and a big.matrix matrix or vector return a big matrix (when the option is FALSE), or return a normal R matrix (TRUE).

The package is built, by default, with R's native BLAS libraries, which use 32-bit signed integer indexing. The default build is limited to vectors of at most 2**31 - 1 entries and matrices with at most 2**31 - 1 rows and 2**31 - 1 columns (note that standard R matrices are limtied to 2**31 - 1 total entries).

The package includes a reference BLAS implementation that supports 64-bit integer indexing, relaxing the limitation on vector lengths and matrix row and column limits. Installation of this package with the 64-bit reference BLAS implementation may be performed from the command-line install:

REFBLAS=1 R CMD INSTALL bigalgebra

where "bigalgebra" is the source package (for example, bigalgebra_0.8.4.tar.gz).

The package may also be build with user-supplied external BLAS and LAPACK libraries, in either 32- or 64-bit varieties. This is an advanced topic that requires additional Makevars modification, and may include adjustment of the low-level calling syntax depending on the library used.

Feel free to contact us for help installing and running the package.

Author(s)

Frédéric Bertrand, Michael J. Kane, Bryan Lewis, John W. Emerson

Maintainer: Frédéric Bertrand <frederic.bertrand@lecnam.net>

References

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

https://www.netlib.org/lapack/

See Also

bigmemory, big.matrix

Examples

# Testing the development of the user-friendly operators:
# if you have any problems, please email us!  - Jay & Mike 4/29/2010

library("bigmemory")
A <- big.matrix(5,4, type="double", init=0,
                dimnames=list(NULL, c("alpha", "beta")))
B <- big.matrix(4,4, type="double", init=0,
                dimnames=list(NULL, c("alpha", "beta")))

C <- A 
D <- A[] 

print(C - D)       # Compare the results (subtraction of an R matrix from a
                   # big.matrix)

# The next example illustrates mixing R and big.matrix objects. It returns by
# default (see # options("bigalgebra.mixed_arithmetic_returns_R_matrix")
D <- matrix(rnorm(16),4)
E <- A 

Class "big.matrix" arithmetic methods

Description

Arithmetic operations for big.matrices

Methods

%*%

signature{x="big.matrix", y="big.matrix"}: ...

%*%

signature{x="matrix", y="big.matrix"}: ...

%*%

signature{x="big.matrix", y="matrix"}: ...

Arith

signature{x="big.matrix", y="big.matrix"}: ...

Arith

signature{x="big.matrix", y="matrix"}: ...

Arith

signature{x="matrix", y="big.matrix"}: ...

Arith

signature{x="big.matrix", y="numeric"}: ...

Arith

signature{x="numeric", y="big.matrix"}: ...

Notes

Miscellaneous arithmetic methods for matrices and big.matrices. See also options("bigalgebra.mixed_arithmetic_returns_R_matrix").

Author(s)

B. W. Lewis <blewis@illposed.net>

Add two double-precision vectors.

Description

Compute double precision DY := DX + DY for N elements, applying the specified storage increments. This routine mirrors the BLAS DAXPY operation with a unit scaling factor.

Usage

dadd(N = NULL, X, INCX = 1L, Y, INCY = 1L)

Arguments

Details

The implementation delegates to the BLAS DAXPY routine with a unit scaling factor, making it equivalent to daxpy(1.0, X, Y) while exposing an interface consistent with other low-level wrappers such as dcopy and dscal.

Value

The modified object Y containing the element-wise sums.

Examples

set.seed(4669)
X <- matrix(runif(6), 3, 2)
Y <- matrix(runif(6), 3, 2)
dadd(X = X, Y = Y)
all.equal(Y, X + Y)

Sum of absolute values

Description

Sum of absolute values

Usage

dasum(N = NULL, X, INCX = 1L)

Arguments

Value

Numeric scalar.

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 big.matrix objects, and A is a scalar.

Usage

daxpy(A = 1, X, Y)

Arguments

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.

Author(s)

Michael J. Kane

References

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

See Also

bigmemory

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()

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.

Usage

dcopy(N = NULL, X, INCX = 1, Y, INCY = 1)

Arguments

Value

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.

Examples

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()

Dot product of two vectors

Description

Dot product of two vectors

Usage

ddot(N = NULL, X, INCX = 1L, Y, INCY = 1L)

Arguments

Value

Numeric scalar containing the dot product.

Examples

ddot(X = 1:3, Y = c(2, 4, 6))

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.

Usage

dgeev(
  JOBVL = NULL,
  JOBVR = NULL,
  N = NULL,
  A,
  LDA = NULL,
  WR,
  WI,
  VL = NULL,
  LDVL = NULL,
  VR = NULL,
  LDVR = NULL,
  WORK = NULL,
  LWORK = NULL
)

Arguments

Value

WR, WI, VR, VL and Work. On exit, A has been overwritten.

Examples

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[,]

rm(A,WR,WI,VL)
gc()

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.

Usage

dgemm(
  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
)

Arguments

Value

Update C with the result.

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()

QR factorization

Description

DGEQRF computes a QR factorization of a real M-by-N matrix A: A = Q * R.

Usage

dgeqrf(
  M = NULL,
  N = NULL,
  A,
  LDA = NULL,
  TAU = NULL,
  WORK = NULL,
  LWORK = NULL
)

Arguments

Value

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.

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()

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.

Usage

dgesdd(
  JOBZ = "A",
  M = NULL,
  N = NULL,
  A,
  LDA = NULL,
  S,
  U,
  LDU = NULL,
  VT,
  LDVT = NULL,
  WORK = NULL,
  LWORK = NULL
)

Arguments

Value

IWORK an integer matrix dimension of (8*min(M,N)) A is updated.

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.

if

JOBZ .ne. 'O', the contents of A are destroyed.

INFO an integer

= 0:

successful exit.

< 0:

if INFO = -i, the i-th argument had an illegal value.

> 0:

DBDSDC did not converge, updating process failed.

Examples

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()

Element-wise (Hadamard) product

Description

Computes Z := X \circ Y. When 'Z' is missing it is allocated automatically.

Usage

dhprod(N = NULL, X, INCX = 1L, Y, INCY = 1L, Z, INCZ = 1L)

Arguments

Value

The updated object 'Z'.

Examples

dhprod(X = 1:4, Y = rep(2, 4))

Euclidean norm (2-norm)

Description

Euclidean norm (2-norm)

Usage

dnrm2(N = NULL, X, INCX = 1L)

Arguments

Value

Numeric scalar containing the Euclidean norm.

Cholesky factorization

Description

DPOTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.

The factorization has the form

A =

U**T * U, if UPLO = 'U', or

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.

Usage

dpotrf(UPLO = "U", N = NULL, A, LDA = NULL)

Arguments

Value

updates the big matrix A with the result, INFO is an integer

= 0:

successful exit

< 0:

if INFO = -i, the i-th argument had an illegal value

> 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.

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()

Product of vector elements

Description

Product of vector elements

Usage

dprdct(N = NULL, X, INCX = 1L)

Arguments

Value

Numeric scalar equal to the product of elements of 'X'.

High-accuracy dot product

Description

Forms the dot product using long-double accumulation to mitigate rounding error.

Usage

dqddot(N = NULL, X, INCX = 1L, Y, INCY = 1L)

Arguments

Value

Numeric scalar equal to the dot product of 'X' and 'Y'.

Examples

dqddot(X = 1:3, Y = c(2, 4, 6))

Scales a vector by a constant.

Description

Scales a vector by a constant.

Usage

dscal(N = NULL, ALPHA, Y, INCY = 1L)

Arguments

Value

Update Y.

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()

Fill a vector or matrix with a constant value

Description

Fill a vector or matrix with a constant value

Usage

dset(N = NULL, ALPHA, X, INCX = 1L)

Arguments

Value

Invisibly returns 'X' after modification.

Examples

x <- matrix(0, 2, 3)
dset(ALPHA = 5, X = x)
x

Element-wise square root

Description

Applies the square root to each entry of 'X' in place, supporting both base R and ['bigmemory::big.matrix'] inputs.

Usage

dsqrt(N = NULL, X, INCX = 1L)

Arguments

Value

Invisibly returns the modified object 'X'.

Examples

vals <- matrix(c(1, 4, 9, 16), 2)
dsqrt(X = vals)
vals

In-place vector subtraction

Description

Updates 'Y' by subtracting 'X', i.e. Y := Y - X.

Usage

dsub(N = NULL, X, INCX = 1L, Y, INCY = 1L)

Arguments

Value

Invisibly returns the modified 'Y'.

Examples

x <- 1:4
y <- rep(10, 4)
dsub(X = x, Y = y)
y

Sum of elements

Description

Sum of elements

Usage

dsum(N = NULL, X, INCX = 1L)

Arguments

Value

Numeric scalar giving the sum of elements of 'X'.

Swap two double-precision vectors.

Description

Exchange the elements of two double precision vectors in place. For I = 0 to N-1, swap DX(LX + I * INCX) with DY(LY + I * INCY) where LX and LY depend on the increment signs. When an optimized BLAS is available the implementation dispatches to DSWAP; otherwise a portable C loop performs the exchange.

Usage

dswap(N = NULL, X, INCX = 1L, Y, INCY = 1L)

Arguments

Details

When an optimized BLAS is available the implementation delegates to the Fortran DSWAP routine. Otherwise a portable C fallback performs the exchange directly while respecting the supplied vector increments.

Value

Invisibly returns NULL; both X and Y are modified in place.

See Also

[dcopy()], [daxpy()] and [dscal()].

Examples

set.seed(4670)
X <- matrix(runif(6), 3, 2)
Y <- matrix(runif(6), 3, 2)
X_original <- X
Y_original <- Y
dswap(X = X, Y = Y)
all.equal(X, Y_original)
all.equal(Y, X_original)

Symmetric matrix-matrix multiplication

Description

Computes C := \alpha \operatorname{op}(A) B + \beta C when 'A' is symmetric.

Usage

dsymm(
  SIDE = "L",
  UPLO = "U",
  M = NULL,
  N = NULL,
  ALPHA = 1,
  A,
  LDA = NULL,
  B,
  LDB = NULL,
  BETA = 0,
  C,
  LDC = NULL
)

Arguments

Value

Invisibly returns 'C'.

Examples

A <- matrix(c(2, 1, 1, 3), 2, 2)
B <- diag(2)
C <- matrix(0, 2, 2)
dsymm(A = A, B = B, C = C)
C

Generalised AXPBY operation

Description

Computes the linear combination Y := \alpha X + \beta Y in place.

Usage

dvcal(N = NULL, ALPHA = 1, X, INCX = 1L, BETA = 1, Y, INCY = 1L)

Arguments

Value

Invisibly returns the modified 'Y'.

Examples

x <- 1:5
y <- rep(2, 5)
dvcal(ALPHA = 2, X = x, BETA = -1, Y = y)
y

Three-dimensional cross product

Description

Three-dimensional cross product

Usage

dxyz(X, Y, Z)

Arguments

Value

The updated object 'Z' containing the cross product.

Examples

dxyz(X = c(1, 0, 0), Y = c(0, 1, 0))

Index of the maximum absolute value

Description

Index of the maximum absolute value

Usage

idamax(N = NULL, X, INCX = 1L)

Arguments

Value

Integer index (1-based).

Index of the minimum absolute value

Description

Index of the minimum absolute value

Usage

idamin(N = NULL, X, INCX = 1L)

Arguments

Value

Integer index (1-based).

Index of the maximum element

Description

Index of the maximum element

Usage

idmax(N = NULL, X, INCX = 1L)

Arguments

Value

Integer index (1-based).

Index of the minimum element

Description

Index of the minimum element

Usage

idmin(N = NULL, X, INCX = 1L)

Arguments

Value

Integer index (1-based) of the smallest entry in 'X'.