Hi Ben,
thank you soo much:) your code worked, and we
managed to halven computing time.
I am copying below the final version of the function,
if anybody is willing to use/improve it.
thanks again,
Moreno
###########################################################
spdiags = function(A){
require(Matrix)
source("getband.R")
# Find all nonzero diagonals
i = rep(seq(1, nrow(A),1),nrow(A));
j = sort(i);
d = sort(j-i);
d = unique(d);
p = length(d); ##the nr. col of the new matrix
m = nrow(A); n = ncol(A);
B = Matrix(FALSE, nrow = min(c(m,n)), ncol = p, sparse = TRUE);
count = 0
for (k in 1:p){
if (m >= n){
i = max(1, 1+d[k]):min(n, m+d[k]);
} else { i = max(1, 1-d[k]):min(m,n-d[k]); }
if (length(i) != 0){
bd = getband(A, d[k]);
if (length(which(bd) == T) != 0){
count = count + 1;
# print( paste( "non-zero diagonal", count) )
B[i,count] = bd
} else {
# print (paste ("ncol", ncol(B), "diag", k ) )
B = B[, -ncol(B)] }
}
}
return (list( B = B, d = d) )
}
############################################################
############################################################
Quoting Ben Bolker <bbol...@gmail.com> on Fri, 13 Apr 2012 20:09:52 +0000:
Ben Bolker <bbolker <at> gmail.com> writes:
I'm not quite sure how to do it, but I think you should look
at the ?band function in Matrix. In combination with diag() of a
suitably truncated matrix, you should be able to extract bands
of sparse matrices efficiently ...
getband <- function(A,k) {
n <- nrow(A)
if (abs(k)>(n-1)) stop("bad band requested")
if (k>0) {
v <- seq(n-k) ## -seq((n-k+1),n)
w <- seq(k+1,n) ## -seq(n-k-1)
} else if (k<0) {
v <- seq(-k+1,n)
w <- seq(n+k)
} else return(diag(A))
diag(band(A,k,k)[v,w,drop=FALSE])
}
PS: I think this should extract the k^th off-diagonal
band in a way that should (?) work reasonably efficiently
with sparse matrices. I have not tested it carefully,
nor benchmarked it.
Ben Bolker
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
