How to compute matrix determinant using MUMPS?


93
views
0
12 weeks ago by
Dear all:

I need to compute the determinant of my stiffness matrix. Using the following code I am able to set up MUMPS solver through PETSc:
A = dolfin.PETScMatrix()

solver = dolfin.PETScKrylovSolver()
solver.ksp().setType("preonly")
solver.ksp().getPC().setType("lu")
solver.ksp().getPC().setFactorSolverPackage("mumps")
solver.ksp().setOperators(A=A.mat())
However, I can't seem to be able to access the Factor Matrix to apply the control commands as in the PETSc example. If I believe this commit by Jan Blechta, the controls commands are there, but in petsc4py 3.8.1 there is no setFactorSetUpSolverType within PC, such that solver.ksp().getPC().getFactorMatrix() returns error code 58. Any help would be super appreciated. Thanks in advance.

Martin
Community: FEniCS Project

1 Answer


3
12 weeks ago by
It looks like the underlying problem is that the factorization has not been computed yet.  The following runs without error for me in the 2017.2 Docker container:

from dolfin import *

mesh = dolfin.UnitSquareMesh(10,10)
V = dolfin.FunctionSpace(mesh,"Lagrange",1)
u = dolfin.TrialFunction(V)
v = dolfin.TestFunction(V)
A = assemble(u*v*dx)

solver = dolfin.PETScKrylovSolver()
solver.ksp().setType("preonly")
solver.ksp().getPC().setType("lu")
solver.ksp().getPC().setFactorSolverPackage("mumps")
solver.ksp().setOperators(A=as_backend_type(A).mat())

# Also works:
#
#B = assemble(v*dx)
#u = Function(V)
#solver.solve(A,u.vector(),B)

x = as_backend_type(Function(V).vector()).vec()
y = as_backend_type(Function(V).vector()).vec()
solver.ksp().getPC().apply(x,y)

solver.ksp().getPC().getFactorMatrix()
​
There may also be some less-involved way to pre-compute the factorization, but nothing stood out to me while browsing the petsc4py API:

http://www.mcs.anl.gov/petsc/petsc4py-current/docs/apiref/index.html
Thanks David! I am now able to get the factorization matrix, which apparently is allocated during factorization and not at solver instantiation. Anyway, I still cannot get the correct determinant:
import dolfin
import numpy
import petsc4py

mesh = dolfin.UnitSquareMesh(1,1)
dx = dolfin.dx(mesh)

V = dolfin.FunctionSpace(mesh, "Lagrange", 1)
u = dolfin.TrialFunction(V)
v = dolfin.TestFunction(V)

A = dolfin.PETScMatrix()
B = dolfin.PETScVector()

dolfin.assemble_system(
    u*v*dx,
    v*dx,
    A_tensor=A,
    b_tensor=B)

print A.array()
print numpy.linalg.det(A.array()) # -> 9.6450617284e-05

solver = dolfin.PETScKrylovSolver()

options = dolfin.PETScOptions()
options.set("ksp_type", "preonly");
options.set("pc_type", "lu");
options.set("pc_factor_mat_solver_package", "mumps");
options.set("mat_mumps_icntl_33", 1);
solver.ksp().setFromOptions()

solver.ksp().setOperators(A=A.mat())

u = dolfin.Function(V)
solver.solve(
    u.vector(),
    B)

print solver.ksp().getPC().getFactorMatrix().getMumpsRinfo(33) # -> 0
print solver.ksp().getPC().getFactorMatrix().getMumpsRinfog(33) # -> 0​
Any idea? Thanks again.

Martin
written 12 weeks ago by Martin Genet  
The following reproduces numpy's determinant (just based on imitating the PETSc example), although I'm not entirely sure what all of the MUMPS flags mean:

fm = solver.ksp().getPC().getFactorMatrix()
rinfo12 = fm.getMumpsRinfog(12)
rinfo13 = fm.getMumpsRinfog(13)
infog34 = fm.getMumpsInfog(34)
print(rinfo12*(2.0**infog34))​
written 12 weeks ago by David Kamensky  
Thanks so much David! Martin
written 12 weeks ago by Martin Genet  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »