Solving a matrix equation in PETSc

General
1

Hi,

I wonder what is the best approach to solve a matrix equation of the form AX=B in PETSc, where A is full rank. I have read about the matSolve method from petsc4py.PETSc, but I cannot find a simple usage example. Any help or workaround would really be appreciated.

Thanks.

2

How do you want to solve this? That each row of B is a new RHS?

3

Each column of B is a new RHS, as you would do with np.linalg.solve

4

Then I would first do an LU decomposition once, by calling PC.setUp() after setting the the ksp type preonly and preconditioner to LU.
Otherwise there is no benefit, as an iterative solver would be dependent on the RHS

1 Like
5

Hi @dokken,

Thanks for your answer.

My best solution is the following:

ksp = PETSc.KSP().create(A.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
for j in range(B.size[1):
    b = B.getDenseColumnVec(j)
    vec = X.getDenseColumnVec(j)
    ksp.solve(b, vec)
    X.restoreDenseColumnVec(j)

but this requires B to be dense… Can you think of a better solution for the case in which B is sparse, and hence I cannot use getDenseColumnVec?

6

Looking at the PETSc documentation:
https://petsc.org/release/manualpages/Mat/MatMatSolve/
you can use: PETSc.Mat.matSolve: src/binding/petsc4py/src/petsc4py/PETSc/Mat.pyx · main · PETSc / petsc · GitLab