Commit 6a5c1447 authored by Reinhard Furrer's avatar Reinhard Furrer

v2.2-0

parent f046e1cd
Package: spam
Type: Package
Title: SPArse Matrix
Version: 2.1-2
Date: 2017-12-21
Version: 2.2-0
Date: 2018-06-19
Authors@R: c(person("Reinhard", "Furrer", role = c("aut", "cre"),
email = "reinhard.furrer@math.uzh.ch"),
person("Florian", "Gerber", role = c("ctb"),
email = "florian.gerber@math.uzh.ch"),
person("Daniel", "Gerber", role = "ctb",
email = "daniel_gerber_2222@hotmail.com"),
email = "florian.gerber@math.uzh.ch"),
person("Roman", "Flury", role = c("ctb"),
email = "roman.flury@math.uzh.ch"),
person("Daniel", "Gerber", role = "ctb",
email = "daniel_gerber_2222@hotmail.com"),
person("Kaspar", "Moesinger", role = "ctb",
email = "kaspar.moesinger@gmail.com"),
email = "kaspar.moesinger@gmail.com"),
person("Youcef", "Saad", role = "ctb",
comment = "SPARSEKIT http://www-users.cs.umn.edu/~saad/software/SPARSKIT/"),
comment = "SPARSEKIT http://www-users.cs.umn.edu/~saad/software/SPARSKIT/"),
person(c("Esmond", "G."), "Ng", role = "ctb",
comment = "Fortran Cholesky routines"),
comment = "Fortran Cholesky routines"),
person(c("Barry", "W."), "Peyton", role = "ctb",
comment = "Fortran Cholesky routines"),
comment = "Fortran Cholesky routines"),
person(c("Joseph", "W.H."), "Liu", role = "ctb",
comment = "Fortran Cholesky routines"),
comment = "Fortran Cholesky routines"),
person(c("Alan", "D."), "George", role = "ctb",
comment = "Fortran Cholesky routines")
)
comment = "Fortran Cholesky routines"),
person(c("Lehoucq", "B."), "Rich", role = "ctb",
comment = "ARPACK"),
person(c("Maschhoff"), "Kristi", role = "ctb",
comment = "ARPACK"),
person(c("Sorensen", "C."), "Danny", role = "ctb",
comment = "ARPACK"),
person(c("Yang"), "Chao", role = "ctb",
comment = "ARPACK"))
Depends:
R (>= 3.1),
dotCall64,
......@@ -42,7 +51,6 @@ Description: Set of functions for sparse matrix algebra.
(2) based on transparent and simple structure(s),
(3) tailored for MCMC calculations within G(M)RF.
(4) and it is fast and scalable (with the extension package spam64).
LazyLoad: Yes
LazyData: Yes
License: LGPL-2
LazyData: true
License: LGPL-2 | BSD_3_clause + file LICENSE
URL: https://git.math.uzh.ch/reinhard.furrer/spam
......@@ -106,6 +106,7 @@ export(
"determinant.spam",
"determinant.spam.chol.NgPeyton",
"det",# << "identical" as base - but with correct determinant()
"chol.spam",
"solve.spam",
"forwardsolve.spam",
......@@ -174,7 +175,8 @@ export(
"bdiag.spam",
"var.spam",
"eigen.spam",
"eigen.spam",
"eigen_approx",
"bandwidth",
......
......@@ -369,7 +369,7 @@ setMethod("t","spam",t.spam)
stop("argument 'x' should be of mode 'numeric' (or 'spam')")
as.spam.spam <- function(x, eps = getOption("spam.eps")) {
force64 <- getOption("spam.force64")
force64 <- getOption("spam.force64")
if(force64)
SS <- .format64()
......@@ -396,7 +396,8 @@ as.spam.spam <- function(x, eps = getOption("spam.eps")) {
nz <- z$rowpointers[dimx[1]+1]-1
if(nz==0) {
return(.newSpam(
dimension=dimx
dimension=dimx,
force64 = force64
))
}
......@@ -466,7 +467,7 @@ as.spam.matrix <- function(x, eps = getOption("spam.eps")) {
"as.spam.numeric" <- function(x, eps = getOption("spam.eps")) {
if (eps<.Machine$double.eps) stop("'eps' should not be smaller than machine precision",call.=FALSE)
force64 = getOption("spam.force64")
force64 <- getOption("spam.force64")
poselements <- (abs(x)>eps)
if (any(!is.finite(x))) {
......@@ -479,7 +480,7 @@ as.spam.matrix <- function(x, eps = getOption("spam.eps")) {
if (identical(nz,0)){
return(
.newSpam(
rowpointers = c(1,rep_len64(2,lx)),
# rowpointers = c(1,rep_len64(2,lx)),
dimension = c(lx,1),
force64 = force64
)
......@@ -506,7 +507,9 @@ as.spam.dist <- function(x, eps = getOption("spam.eps")) {
dimx <- attr(x,"Size")
size <- dimx*(dimx-1)/2
if( getOption("spam.force64") || dimx^2 > 2147483647)
force64 <- getOption("spam.force64")
if(force64 || dimx^2 > 2147483647)
SS <- .format64()
else
SS <- .format32
......@@ -532,15 +535,18 @@ as.spam.dist <- function(x, eps = getOption("spam.eps")) {
PACKAGE = SS$package
)
nz <- z$rowpointers[dimx+1]-1
if(nz==0) return( .newSpam( rowpointers = c( 1, rep( 2, dimx)),
dimension = c( dimx, dimx)))
if(nz==0) return( .newSpam(
# rowpointers = c( 1, rep( 2, dimx)),
dimension = c( dimx, dimx),
force64 = force64))
return(.newSpam(
entries=z$entries[1:nz],
colindices=z$colindices[1:nz],
rowpointers=z$rowpointers[1:(dimx+1)],
dimension=c(dimx,dimx)
dimension=c(dimx,dimx),
force64 = force64
))
}
......@@ -622,7 +628,8 @@ spam.numeric <- function(x, nrow = 1, ncol = 1, eps = getOption("spam.eps")) {
return( .newSpam( entries = entries,
colindices = colindices,
rowpointers = rowpointers,
dimension = dimx))
dimension = dimx,
force64 = force64))
}
setOldClass(c("dist", "numeric"))
......@@ -1141,7 +1148,7 @@ function (x, rw, cl,value)
yrow,
ycol,
y,
y = vector_dc("double",nrow*ycol),
x@entries,
x@colindices,
......@@ -1241,51 +1248,46 @@ function (x, rw, cl,value)
if(nz==0){#trap zero matrix
z$entries <- 0
z$colindices <- 1
z$rowpointers <- c(1,rep(2,xn))
} else {
entries <- z$entries
colindices <- z$colindices
rowpointers <- z$rowpointers
z <- NULL
length(entries) <- nz
length(colindices) <- nz
z <- .C64("sortrows",
SIGNATURE=c(SS$signature, "double", SS$signature, SS$signature),
xn,
entries=entries,
colindices=colindices,
rowpointers=rowpointers,
INTENT=c("r", "rw", "rw", "rw"),
NAOK = getOption("spam.NAOK"),
PACKAGE = SS$package)
return(.newSpam(
dimension=c(xn,yl)
))
}
entries <- z$entries
colindices <- z$colindices
rowpointers <- z$rowpointers
z <- NULL
length(colindices) <- nz
entries <- z$entries
colindices <- z$colindices
rowpointers <- z$rowpointers
z <- NULL
length(entries) <- nz
length(colindices) <- nz
z <- .C64("sortrows",
SIGNATURE=c(SS$signature, "double", SS$signature, SS$signature),
xn,
entries=entries,
colindices=colindices,
rowpointers=rowpointers,
INTENT=c("r", "rw", "rw", "rw"),
NAOK = getOption("spam.NAOK"),
PACKAGE = SS$package)
return(.newSpam(
entries=entries,
colindices=colindices,
rowpointers=rowpointers,
entries=z$entries,
colindices=z$colindices,
rowpointers=z$rowpointers,
dimension=c(xn,yl)
))
}
.spam.matmul.vector <- function(x,y)
{
# --- CHANGED ---
# Refactored the function .spam.matmul into the functions
# .spam.matmul and .spam.matmul.vector
# if we have x*Y, we calculate t(t(Y)*x)
# Use "is.spam(y)" instead of "is.vector(x)" because the later might be
# misleading in case that x has any attributes attached.
......@@ -1297,26 +1299,26 @@ function (x, rw, cl,value)
b <- y
}
SS <- .format.spam(A)
nrow <- A@dimension[1]
ncol <- A@dimension[2]
if(length(b) != ncol) stop("not conformable for multiplication")
z <- .C64("amux",
NAOK = getOption("spam.NAOK"),
SIGNATURE=c(SS$signature, "double", "double", "double",
SS$signature, SS$signature),
nrow,
b,
y = vector_dc("double",nrow),
A@entries,
A@colindices,
A@rowpointers,
INTENT=c("r", "r", "w", "r", "r", "r"),
PACKAGE = SS$package)$y
if(is.spam(y))
dim(z) <- c(1,nrow)
else
......
......@@ -117,9 +117,9 @@
## print("2")
return(
.newSpam(
entries=x@entries,
colindices=x@colindices,
rowpointers=c(1,rep_len64(2,value[1])),
# entries=x@entries,
# colindices=x@colindices,
# rowpointers=c(1,rep_len64(2,value[1])),
dimension=value,
force64=force64
)
......
......@@ -66,6 +66,8 @@ nearest.dist <- function( x, y=NULL, method = "euclidean",
if (is.na(method)) stop("invalid distance method")
force64 <- getOption("spam.force64")
if (method == 4) {
if (is.null(R))
p <- ifelse( miles,3963.34,6378.388)
......@@ -115,7 +117,7 @@ nearest.dist <- function( x, y=NULL, method = "euclidean",
}
# EXPLICIT-STORAGE-FORMAT
if(max(n1,n2,nnz) > 2147483647 - 1 | getOption("spam.force64"))
if(max(n1,n2,nnz) > 2147483647 - 1 || force64)
SS <- .format64()
else
SS <- .format32
......@@ -161,7 +163,7 @@ nearest.dist <- function( x, y=NULL, method = "euclidean",
nnz <- nnz*getOption("spam.nearestdistincreasefactor")*n1/(d$iflag-1)
# EXPLICIT-STORAGE-FORMAT
if(max(n1,n2,nnz) > 2147483647 - 1 | getOption("spam.force64"))
if(max(n1,n2,nnz) > 2147483647 - 1 || force64)
SS <- .format64()
else
SS <- .format32
......@@ -180,11 +182,19 @@ nearest.dist <- function( x, y=NULL, method = "euclidean",
length(d$entries) <- d$nnz
length(d$colindices) <- d$nnz
if(d$nnz == 0) {
return(.newSpam(
dimension=c(n1,n2),
force64 = force64
))
}
return(.newSpam(
entries=d$entries,
colindices=d$colindices,
rowpointers=d$rowpointers,
dimension=c(n1,n2)
dimension=c(n1,n2),
force64 = force64
))
}
......
# HEADER ####################################################
# This is file spam/R/eigen.R. #
# It is part of the R package spam, #
# --> https://CRAN.R-project.org/package=spam #
# --> https://CRAN.R-project.org/package=spam64 #
# --> https://git.math.uzh.ch/reinhard.furrer/spam #
# by Reinhard Furrer [aut, cre], Florian Gerber [ctb], #
# Daniel Gerber [ctb], Kaspar Moesinger [ctb], #
# Youcef Saad [ctb] (SPARSEKIT), #
# Esmond G. Ng [ctb] (Fortran Cholesky routines), #
# Barry W. Peyton [ctb] (Fortran Cholesky routines), #
# Joseph W.H. Liu [ctb] (Fortran Cholesky routines), #
# Alan D. George [ctb] (Fortran Cholesky routines), #
# Esmond G. Ng [ctb] (Fortran Cholesky routines), #
# Barry W. Peyton [ctb] (Fortran Cholesky routines), #
# Joseph W.H. Liu [ctb] (Fortran Cholesky routines), #
# Alan D. George [ctb] (Fortran Cholesky routines) #
# HEADER END ################################################
setMode <- function(sMode, symmetric, silent = FALSE){
if (symmetric) {
modes <- c("LM", "SM", "LA", "SA") #,"BE")
if (!(sMode %in% modes)) {
sMode <- 'LM'
if (!silent) { warning(paste("the control option \"mode\" is set to", sMode), call. = FALSE) } }
} else {
modes <- c("LR", "SR", "LI", "SI", "SM", "LM")
if (!(sMode %in% modes)) {
sMode <- 'LR'
if(!silent) { warning(paste("the control option \"mode\" is set to", sMode), call. = FALSE) } }
}
if (sMode == 'LM') { imode <- as.integer(1) }
if (sMode == 'SM') { imode <- as.integer(2) }
if (sMode == 'LR') { imode <- as.integer(3) }
if (sMode == 'SR') { imode <- as.integer(4) }
if (sMode == 'LI') { imode <- as.integer(5) }
if (sMode == 'SI') { imode <- as.integer(6) }
if (sMode == 'LA') { imode <- as.integer(7) }
if (sMode == 'SA') { imode <- as.integer(8) }
if (sMode == 'BE') { imode <- as.integer(9) }
return(imode)
}
mk_cmplxentries <- function(z)
{
if (dim(z)[2] != 2)
stop("wrong format from fortran return: dn_eigen_f")
cmplx <- NULL
cmplx <- sapply(1:length(z[ , 1]), function(x) { complex(real = z[x, 1], imaginary = z[x, 2]) })
if (is.null(cmplx))
stop("error while aggregating fortran return from two vectors to a complex")
return(cmplx)
}
getEigenval <- function(values, mode, dim, nEig, symmetric){
if (symmetric) {
result <- rep(NA, dim)
orderedInd <- order(values[1:nEig], decreasing = TRUE)
values <- values[orderedInd]
# sort works also for complex values
if (identical(mode %% 2, 0)) { # even modes are SA, SM, SR, SI
result[(dim - nEig + 1):dim] <- values
} else { # odd modes are LA, LA, LR, LI
result[1:nEig] <- values }
} else {
result <- matrix(NA, nrow = dim, ncol = 2)
orderedInd <- order(values[1:nEig, 1], decreasing = TRUE)
values <- values[orderedInd, ]
if (identical(mode %% 2, 0)) {
result[(dim - nEig + 1):dim, ] <- values
} else {
result[1:nEig, ] <- values }
result <- mk_cmplxentries(result)
}
return(list("eigenvalues" = result , "order" = orderedInd))
}
getEigenvec <- function(v, sym, dimen, nEig, orderind, eigenvalues, cmplxeps){
if (is.null(v))
stop("fortran returned NULL eigenvectors")
if (sym) {
v <- matrix(v, nrow = dimen, ncol = nEig, byrow = FALSE)
v <- v[ , orderind, drop = FALSE]
} else {
v <- matrix(v[1:(dimen*nEig*2)], nrow = dimen, ncol = nEig*2, byrow = FALSE)
v_real <- v[ , seq(1, nEig*2, by = 2)]
v_imag <- v[ , seq(2, nEig*2, by = 2)]
v_real <- v_real[ , orderind, drop = FALSE]
v_imag <- v_imag[ , orderind, drop = FALSE]
rm(v)
eigenvalues <- stats::na.omit(eigenvalues)
v <- matrix(NA, nrow = dimen, ncol = nEig)
v <- sapply(1:nEig, function(x) {
if (abs(Im(eigenvalues[x])) > cmplxeps) {
mk_cmplxentries(cbind(v_real[ , x], v_imag[ , x]))
} else {
mk_cmplxentries(cbind(v_real[ , x], base::rep.int(0, times = dimen)))
}
})
}
return(v)
}
eigen_approx <- function(x,
nev,
ncv,
nitr,
mode,
only.values = FALSE,
verbose = FALSE,
f_routine){
devMode <- FALSE
# check & parse arguments
if (x@dimension[1] <= nev)
stop("nev: the number of eigenvalues to calculate must be smaller than the matrix dimensions", call. = devMode)
if (f_routine != "ds_eigen_f" && f_routine != "dn_eigen_f")
stop("non valid fortran routine is specified", call. = devMode)
f_mode <- setMode(sMode = mode, symmetric = ifelse(identical(f_routine, "ds_eigen_f"), TRUE, FALSE))
fortran_object <- result <- list(NULL)
if(getOption("spam.force64") || .format.spam(x)$package != "spam") {
SS <- .format64()
} else {
SS <- .format32 }
# Fortran call: symmetric matrices
if (f_routine == "ds_eigen_f") {
fortran_object <- .C64 (f_routine,
SIGNATURE = c(SS$signature, SS$signature, SS$signature, SS$signature,
SS$signature, SS$signature, "double",
SS$signature, SS$signature, "double",
"double", SS$signature, SS$signature),
maxnev = nev,
ncv = ncv,
maxitr = nitr,
n = x@dimension[1],
iwhich = f_mode,
na = x@dimension[1],
a = x@entries,
ja = x@colindices,
ia = x@rowpointers,
v = vector_dc("double", x@dimension[1]*ncv),
d = vector_dc("double", nev),
vf = ifelse(devMode, 1L, 0L),
iparam = integer_dc(11),
INTENT = c("r", "r", "r", "r", "r", "r", "r", "r", "r",
"rw", "rw", "r", "rw"),
NAOK = getOption("spam.NAOK"),
PACKAGE = SS$package)
if (is.null(fortran_object)) {
stop("error while calling fortran routine, check (control) arguments", call. = devMode) }
result <- list ("nEigenVal" = nev,
"Mode" = f_mode,
"Eigenvectors" = if (!only.values) { fortran_object$v[1:(x@dimension[1]*nev)] } else { NULL },
"Eigenvalues" = fortran_object$d,
"nconv" = fortran_object$iparam[5])
}
# Fortran call: nonsymmetric matrices
if (f_routine == "dn_eigen_f") {
fortran_object <- .C64 (f_routine,
SIGNATURE = c(SS$signature, SS$signature, SS$signature, SS$signature,
SS$signature, SS$signature, "double",
SS$signature, SS$signature, "double",
"double", "double", SS$signature,SS$signature),
maxnev = nev,
ncv = ncv,
maxitr = nitr,
n = x@dimension[1],
iwhich = f_mode,
na = x@dimension[1],
a = x@entries,
ja = x@colindices,
ia = x@rowpointers,
v = vector_dc("double", x@dimension[1]*ncv),
dr = vector_dc("double", nev),
di = vector_dc("double", nev),
vf = ifelse(devMode, 1L, 0L),
iparam = integer_dc(11),
INTENT = c("r", "r", "r", "r", "r", "r", "r", "r", "r",
"rw", "rw", "rw", "r", "rw"),
NAOK = getOption("spam.NAOK"),
PACKAGE = SS$package)
if (is.null(fortran_object)) {
stop("error while calling fortran routine, check (control) arguments", call. = devMode) }
result <- list ("nEigenVal" = nev,
"Mode" = f_mode,
"Eigenvectors" = if (!only.values) { fortran_object$v[1:(x@dimension[1]*nev*2)] } else { NULL },
"Eigenvalues" = cbind(fortran_object$dr, fortran_object$di),
"nconv" = fortran_object$iparam[5])
}
if (verbose) {
cat("\n used options/arguments:")
if (identical(f_routine, "dn_eigen_f") ) { issym <- "non symmetric" } else { issym <- "symmetric" }
cat(paste("\n FORTRAN routine:", issym, "matrices"))
cat(paste("\n nitr:", nitr))
cat(paste("\n ncv:", ncv, "\n"))
cat("\n summary FORTRAN return:")
cat(paste("\n ", nev, "eigenvalues requested and", result$nconv, "converged\n"))
}
if (nev != result$nconv)
warning(paste("\n only", result$nconv, "instead of", nev ,"eigenvalues converged, try to increase 'control = list(nitr = .., ncv = ..)'"), call. = devMode)
if (is.null(result)) {
stop("unpredicted error, check control options of the eigen.spam function.", call. = devMode)
}
return(result)
}
eigen.spam <- function (x, nev = 10, symmetric, only.values = FALSE, control = list()){
con <- list(mode = 'LM',
nitr = NULL,
ncv = NULL,
spamflag = FALSE,
verbose = FALSE,
cmplxeps = .Machine$double.eps)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC]))
warning("unknown names in control: ", paste(noNms, collapse = ", "), call. = TRUE)
ifelse(!con$verbose, vFlag <- FALSE, vFlag <- TRUE)
# arpack routines cant handle 'small' matrices
minDimARPACK <- 100
resContainer <- list(NULL)
if (!is.spam(x)) {
try(x <- as.spam(x)) }
if (missing(symmetric)){
symmetric <- isSymmetric.spam(x)
} else {
if (!identical(symmetric, isSymmetric.spam(x))) {