[R] writing spdiags function for R

[Ben Bolker]

On 12-04-16 09:32 PM, Moreno I. Coco wrote:
> 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

  A couple of quick comments:

 the semicolons at the end of each line are unnecessary (they're
harmless, but considered bad style -- most common in code of C and
MATLAB coders

  You don't need to source("getband.R") inside the body of spdiags,
probably better to do it outside (so it doesn't get run every time the
function is called)

if (length(which(bd) == T) != 0)

  can be shortened to

  if (any(bd!=0))

  I haven't looked too carefully but can't all the i/j/d/p stuff be
replaced by

  n <- nrow(A)
  p <- 2*n-1
  d <- seq(-(n-1),(n-1))

?


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