Some basic features of 'MUMPS' (Multifrontal Massively Parallel
sparse direct Solver) are wrapped in a class whose methods can be used
for sequentially solving a sparse linear system (symmetric or not)
with one or many right hand sides (dense or sparse).
There is a possibility to do separately symbolic analysis,
LU (or LDL^t) factorization and system solving.
Third part ordering libraries are included and can be used: 'PORD', 'METIS', 'SCOTCH'.
'MUMPS' method was first described in Amestoy et al. (2001)
MUMPS stands for "a MUltifrontal Massively Parallel sparse direct Solver" see more on their official site http://mumps.enseeiht.fr/. Currently, it is one of the most competitive direct solvers for sparse matrices. On my CPU (Xenon E5-2609 v2 @ 2.50GHz) I have a speedup ranging from 3 to 16 compared to the default solver from Matrix package. In addition, the precision of the solution is equal or even better than with Matrix.
Example of use: library(Matrix) library(rmumps) n=2000 a=Matrix(0, n, n) set.seed(7) ij=sample(1:(nn), 15n) a[ij]=runif(ij) diag(a)=0 diag(a)=-rowSums(a) a[1,1]=a[1,1]-1 am=Rmumps$new(a) b=as.double(a%%(1:n)) # rhs for an exact solution vector 1:n # following time includes symbolic analysis, LU factorization and system solving system.time(x<-am$solve(b)) bb=2b # this second time should be much shorter # as symbolic analysis and LU factorization are already done system.time(xx<-am$solve(bb)) # compare to Matrix corresponding times system.time(xm<-solve(a, b)) system.time(xxm<-solve(a, bb)) # compare to Matrix precision range(x-1:n) # mumps range(xm-1:n) # Matrix
# matrix inversion system.time(aminv <- am$inv()) system.time(ainv <- solve(a)) # the same in Matrix # clean up by hand to avoid possible interference between gc() and # Rcpp object destructor after unloading this namespace rm(am)