Givens-based QR decomposition in R or C++: poor performance implementation (or lack thereof) for rectangular matrices

My own Givens-based QR decomposition function in R (pseudocode from Turin's University linear algebra department) is the following:

>Givens.fn<-function(V) {

   tol<-1.0*10^-14
   m<-dim(V)[1]
   n<-dim(V)[2]

   spId<-bandSparse(m,m,0,list(rep(1, m+1))) # computed just once 
   Q<-spId                                   
   R<-V
   sapply(1:n, function(j){                      
      if (j<m) {    
         vapply((j+1):m, function(i){          # vectorized internal loop  
            if(abs(R[i,j])>tol) {                             
               G<-spId                  
               x<-R[j,j]
               y<-R[i,j]
               norm<-sqrt(x^2+y^2)
               c<-x/norm
               s<-y/norm
               G[j,j]=c
               G[i,i]=c
               G[j,i]=s
               G[i,j]=-s                         
               Q<<-G%*%Q
               R<<-G%*%R
            } 
            return(1)                   # --> saves 15% execution time!
         },FUN.VALUE=1.0)
      }
   })

  Q<-t(Q)

return(list(Q,R))
}

It works fine, but the givens() function provided in R by pracma is 134x faster on an total elapsed time around 10" (on my intel i5-3570 4core, 8GB desktop):

>m<-20; n<-20
>set.seed(1)
>X <- as.matrix(replicate(n, runif(m)))
>library(rbenchmark)
>library(pracma)
>benchmark(Givens.fn(X),
      givens(X),
      order = "elapsed",
      replications = 10)

          test replications elapsed relative user.self sys.self user.child sys.child
2    givens(X)           10    0.08    1.000      0.08     0.00         NA        NA
1 Givens.fn(X)           10   10.77  134.625     10.69     0.08         NA        NA

Unfortunately, the pracma version accepts just square matrices, which is not my case. Any further code improvements? Any other Givens-based QR functions (even in C++)? Thanks.