Package {rmumps}

Rcpp port of MUMPS library for LU or LDL^t factorization of sparse matrices

Description

Creates a MUMPS compatible object storing a sparse matrix. Gives a possibility to do separately symbolic analysis, factorization and system solving.

Details

Create a new Rmumps object with A <- Rmumps$new(asparse) then solve a linear system with one or many right hand sides x <- solve(A, b). Cf. Rmumps

Author(s)

Serguei Sokol, INRA

Maintainer: Serguei Sokol (sokol at insa-toulouse.fr)

References

MUMPS official site https://mumps-solver.org

Sokol S. rmumps: Wrapper for 'MUMPS' Library. https://CRAN.R-project.org/package=rmumps

Examples

  # 3 x 3 system given as triplets (0-based indices): A = diag(4, 5, 6)
  A <- Rmumps$new(0:2, 0:2, c(4, 5, 6), 3L)
  x <- solve(A, c(4, 10, 18))
  x  # 1 2 3

Exported Constants

Description

Integer constants defining permutation types and exported from rmumps are following:

• RMUMPS_PERM_AMD

• RMUMPS_PERM_AMF

• RMUMPS_PERM_SCOTCH

• RMUMPS_PERM_PORD

• RMUMPS_PERM_METIS

• RMUMPS_PERM_QAMD

• RMUMPS_PERM_AUTO

They are all regrouped in a named vector RMUMPS_PERM where names are items above and values are corresponding constants.

Examples

if (requireNamespace("slam", quietly=TRUE)) {
    am=rmumps::Rmumps$new(slam::as.simple_triplet_matrix(diag(1:3)))
    am$set_permutation(RMUMPS_PERM_AMD)
    am$solve(1:3)
}

Rcpp Exported Class Wrapping MUMPS library

Description

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.

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

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).

Author(s)

Serguei Sokol, INRA

References

MUMPS official site https://mumps-solver.org

Sokol S. rmumps: Wrapper for 'MUMPS' Library. https://CRAN.R-project.org/package=rmumps

Examples

 if (requireNamespace("Matrix", quietly = TRUE)) {
  # prepare random sparse matrix
  library(Matrix)
  library(rmumps)
  n=300
  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()
 }

Delete via Pointer

Description

This is a C wrapper to Rmumps::~Rmumps() destructor. Available in R too. In C++ code can be used as rmumps::Rmumps__del_ptr(pm)

Usage

Rmumps__del_ptr(pm)

Arguments

Value

No return value, called for its side effect: the object and its MUMPS instance are freed; pm must not be used afterwards.

Get Permutation Parameter

Description

This is a C wrapper to Rmumps::get_permutation() method. Available in R too. In C++ code can be used as rmumps::Rmumps__get_permutation(pm)

Usage

Rmumps__get_permutation(pm)

Arguments

Value

integer defining permutation method used before matrix decomposition.

Pointer Interface to Rmumps Methods

Description

C wrappers to methods of the Rmumps class, for objects held as pointers (e.g. created with Rmumps__ptr_ijv or obtained with Rmumps__xptr). Available in R too. In C++ code they can be used as rmumps::Rmumps__symbolic(pm) etc. (LinkingTo: rmumps, Rcpp).

Usage

Rmumps__xptr(obj)

Rmumps__symbolic(pm)

Rmumps__numeric(pm)

Rmumps__solvetptr(pm, pb, lrhs, nrhs)

Rmumps__set_icntl(pm, iv, ii)

Rmumps__get_icntl(pm)

Rmumps__set_cntl(pm, v, iv)

Rmumps__get_cntl(pm)

Rmumps__get_infos(pm)

Rmumps__det(pm)

Rmumps__log2det(pm)

Rmumps__get_sym_perm(pm)

Rmumps__get_uns_perm(pm)

Rmumps__set_perm_in(pm, p_in)

Rmumps__mumps_version(pm)

Arguments

Value

Rmumps__xptr: a pointer of type XPtr<Rmumps> to the object inside obj, without a finalizer (obj keeps owning it and must stay alive). Rmumps__get_icntl, Rmumps__get_cntl, Rmumps__get_infos, Rmumps__det, Rmumps__log2det, Rmumps__get_sym_perm, Rmumps__get_uns_perm, Rmumps__mumps_version: the same values as the corresponding methods of Rmumps. The other functions are called for their side effects and return NULL invisibly.

Examples

a <- Rmumps$new(c(0L, 1L), c(0L, 1L), c(2, 4), 2L)
pm <- Rmumps__xptr(a)
Rmumps__numeric(pm)
Rmumps__det(pm)

Construct via Triplet Pointers

Description

This is a C wrapper to Rmumps::Rmumps(i, j, v, n, nz, sym) constructor. Available in R too. In C++ code can be used as rmumps::Rmumps__ptr_ijv(pi, pj, pa, n, nz, sym)

Usage

Rmumps__ptr_ijv(pi, pj, pa, n, nz, sym)

Arguments

Value

pointer of type XPtr<Rmumps> pointing to newly created object. To avoid memory leakage, it is user's responsibility to call Rmumps__del_ptr(pm) in a due moment (where pm is the returned pointer).

Set Matrix via Pointer

Description

This is a C wrapper to Rmumps::set_mat_ptr(a) method. Available in R too. In C++ code can be used as rmumps::Rmumps__set_mat_ptr(pm). Using this method invalidates previous numeric decomposition (but not symbolic one).

Usage

Rmumps__set_mat_ptr(pm, pa)

Arguments

Value

No return value, called for its side effect: the object now holds the new matrix values (the next solve refactorizes).

Set Permutation Parameter

Description

This is a C wrapper to Rmumps::set_permutation(permutation) method. Available in R too. In C++ code can be used as rmumps::Rmumps__set_permutation(pm, permutation)

Usage

Rmumps__set_permutation(pm, permutation)

Arguments

Value

No return value, called for its side effect: the ordering method (ICNTL(7)) of the object is set.

Solve via Pointer

Description

This is a C wrapper to Rmumps::solveptr() method. Available in R too. In C++ code can be used as rmumps::Rmumps__solveptr(pobj, pb, lrhs, nrhs)

Usage

Rmumps__solveptr(pobj, pb, lrhs, nrhs)

Arguments

Value

No return value, called for its side effect: the right hand sides pointed to by pb are overwritten by the solutions.

Explore via Triplet

Description

This is a C wrapper to Rmumps::triplet() method. Available in R too. In C++ code can be used as rmumps::Rmumps__triplet(pm)

Usage

Rmumps__triplet(pm)

Arguments

Value

a list with sparse triplet described with fields i, j, v

Solve a Transposed System

Description

Solves A^T x = b for an Rmumps object holding A.

Usage

solvet(a, b, ...)

Arguments

Value

The solution: a vector or a matrix, like solve().

Examples

a <- Rmumps$new(c(0L, 1L, 0L), c(0L, 1L, 1L), c(2, 4, 1), 2L)
solvet(a, c(2, 9))