[R] writing spdiags function for R

[Moreno I. Coco]

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