iramuteq.org Git - iramuteq/blob - 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 my.svd <- function(x, nu, nv) {
  25   stopifnot(nu==nv)
  26   rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
  27   if (rc != 0)
  28     stop("Couldn't run external svd code")
  29   d <- scan("/tmp/sout-S", skip=1)
  30 #FIXME : sometimes, libsvdc doesn't find solution with 2 dimensions, but does with 3
  31   if (length(d)==1) {
  32       nu <- nu + 1
  33       rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
  34       d <- scan("/tmp/sout-S", skip=1)
  35   }
  36   ut <- matrix(scan('/tmp/sout-Ut',skip=1),nrow=nu,byrow=TRUE)
  37   if (nrow(ut)==3) {
  38       ut <- ut[-3,]
  39   }
  40   vt <- NULL
  41   list(d=d, u=-t(ut), v=vt)
  42 }