...
[iramuteq] / Rscripts / mysvd.R
1 #http://www.mail-archive.com/rcpp-devel@lists.r-forge.r-project.org/msg01513.html
2
3 write.sparse <- function (m, to) {
4   ## Writes in a format that SVDLIBC can read
5   stopifnot(inherits(m, "dgCMatrix"))
6   fh <- file(to, open="w")
7   
8   wl <- function(...) cat(..., "\n", file=fh)
9   
10   ## header
11   wl(dim(m), length(m@x))
12   
13   globalCount <- 1
14   nper <- diff(m@p)
15   for(i in 1:ncol(m)) {
16     wl(nper[i])  ## Number of entries in this column
17     if (nper[i]==0) next
18     for(j in 1:nper[i]) {
19       wl(m@i[globalCount], m@x[m@p[i]+j])
20       globalCount <- globalCount+1
21     }
22   }
23 }
24
25 write.sparse2 <- function(m, to) {
26     stopifnot(inherits(m, "dgCMatrix"))
27     fh <- file(to, open="w")
28   
29     wl <- function(...) cat(..., "\n", file=fh)
30     wl(dim(m), length(m@x))   
31
32 }
33
34 print('FIXME : svdpath hardcoded !!')
35 my.svd <- function(x, nu, nv) {
36   stopifnot(nu==nv)
37   rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
38   if (rc != 0)
39     stop("Couldn't run external svd code")
40   d <- scan("/tmp/sout-S", skip=1)
41 #FIXME : sometimes, svdlibc doesn't find solution with 2 dimensions, but does with 3 
42   if (length(d)==1) {
43       nu <- nu + 1
44       rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
45       d <- scan("/tmp/sout-S", skip=1)
46   }
47   ut <- matrix(scan('/tmp/sout-Ut',skip=1),nrow=nu,byrow=TRUE)
48   if (nrow(ut)==3) {
49       ut <- ut[-3,]
50   }
51   vt <- NULL
52   list(d=d, u=-t(ut), v=vt)
53 }