Rmumps function

Rcpp Exported Class Wrapping MUMPS library

This class can be used for storing sparse matrix and solving corresponding linear system with one or many right hand sides. There is a possibility to do separately symbolic analysis, LU factorization and system solving. 1.1

class

Examples

## Not run:

  # prepare random sparse matrix
  library(Matrix)
  library(rmumps)
  n=2000
  a=Matrix(0, n, n)
  set.seed(7)
  ij=sample(1:(n*n), 15*n)
  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<-solve(am, b))
  bb=2*b
  # this second time should be much shorter
  # as symbolic analysis and LU factorization are already done
  system.time(xx<-solve(am, 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 <- solve(am))
  system.time(ainv <- solve(a)) # the same in Matrix
  
  # symmetric matrix
  asy=as(a+t(a), "symmetricMatrix")
  bs=as.double(asy%*%(1:n)) # rhs for 1:n solution
  au=asy
  # Here, we keep only diagonal and upper values of asy matrix.
  # It could be also diagonal and lower values.
  au[row(au)>col(au)]=0
  ams=Rmumps$new(au, sym=1)
  system.time(xs<-solve(ams, bs)) # rmumps
  system.time(xsm<-solve(asy, bs))# Matrix
  # compare to Matrix precision
  range(xs-1:n)  # mumps
  range(xsm-1:n) # Matrix

  # clean up by hand to avoid possible interference between gc() and
  # Rcpp object destructor after unloading this namespace
  rm(am, ams)
  gc()
## End(Not run)

References

MUMPS official site http://mumps.enseeiht.fr

Sokol S (2020). Rmumps: Rcpp port of MUMPS. rmumps package version 5.2.1-X, <URL: http://CRAN.R-project.org/package=rmumps>.

Author(s)

Serguei Sokol, INRA

Note

When creating a symmetric matrix (sym=1 or sym=2), the upper (or lower) mart of the input matrix must be zeroed.

For meaning of entries in MUMPS vectors cntl, icntl, info, rinfo, infog and rinfog cf. original documentation of MUMPS project.

No need to call symbolic() and numeric() methods before a solve() call.

If in constructor, a parameter copy is set to FALSE, no rhs neither matrix copying is done. The solution is written "in place" thus overwriting rhs (watch out side effects)

For a detailed error diagnostic (e.g. when factorizing a singular matrix), use method get_infos() and cf. MUMPS documentation on the official MUMPS site).

Fields

• sym:: integer (read only), 0=non symmetric matrix, 1=symmetric with pivots on diagonal or 2=general symmetric
• copy:: logical, copy or not rhs and matrix values
• mrhs:: numeric matrix, multiple rhs (always overwritten with solution)
• rhs:: numeric vector, single rhs (always overwritten with solution)

Methods

• new(asp, sym=0, copy=TRUE):: constructor from Matrix::dgTMatrix class (or from convertible to it) and slam::simple_triplet_matrix class
• new(i, j, x, n, copy=TRUE):: constructor from triade rows, cols, vals
• symbolic():: do symbolic analysis (stored internally)
• numeric():: do LU or LDL^t factorization (stored internally)
• solve(b):: solve single rhs (if b is a vector) or multiple rhs if b is a matrix (can be dense or sparse). Return the solution(s).
• solvet(b):: same as solve() but solves with transposed matrix
• det():: Return determinant of the matrix
• inv():: Return inverse of the matrix)
• set_mat_data(x):: updates matrix entries (x must be in the same order as in previous calls
• set_icntl(iv, ii):: set ICNTL parameter vector
• get_icntl():: get ICNTL parameter vector
• set_cntl(v, iv):: set CNTL parameter vector
• get_cntl():: get CNTL parameter vector
• get_infos():: get a named list of information vectors: info, rinfo, infog and rinfog
• dim():: Return a dimension vector of the matrix
• nrow():: Return a row number of the matrix
• ncol():: Return a column number of the matrix
• print():: Print summary information on the matrix
• show():: Print summary information on the matrix
• set_keep():: Set KEEP array elements (undocumented feature of MUMPS)
• get_keep():: Get a copy of KEEP array elements (length=500)
• set_permutation(perm):: Set permutation type which can impact storage and factorization performances. Parameter perm can take one of the following predefined integer values RMUMPS_PERM_AMD, RMUMPS_PERM_AMF, RMUMPS_PERM_SCOTCH, RMUMPS_PERM_PORD, RMUMPS_PERM_METIS, RMUMPS_PERM_QAMD. This method should be called once and before symbolic analysis of the matrix. If it is called afterward, a new symbolic and numeric factorization will be performed when one of other methods (e.g. solve()) will request them. In other words, previous symbolic and numeric factorizations are canceled by this method.
• get_permutation():: get permutation type currently set in the object
• mumps_version():: Return a string with MUMPS version used in rmumps