[petsc-users] Computing part of the inverse of a large matrix
Zhang, Hong
hzhang at mcs.anl.gov
Tue Oct 8 09:54:24 CDT 2019
Jan:
Mumps manual says:
--3 MUMPS was called with an invalid value for JOB. This may happen if the analysis (JOB=1)was not performed (or failed) before the factorization (JOB=2), or the factorization was not performed (or failed) before the solve (JOB=3), or the initialization phase (JOB=-1) was performed a second time on an instance not freed (JOB=-2). See description of JOB in Section 4. This error also occurs if JOB does not contain the same value on all processes on entry to MUMPS. INFO(2) is then set to the local value of JOB.
However, your code seems identical to petsc/src/ksp/ksp/examples/tests/ex27.c which works fine. Can you test your matrices using ex27.c?
Hong
Hey,
I tried your proposed way. The factorization seems to work but the I get an error when solving A*X=B. When I run the appended code
I get the following error message:
Traceback (most recent call last):
File "inverse_matrix.py", line 145, in <module>
compute_inverse_mat(dynamical_matrix_nn, args.dynamical_matrix_dimension, args.m_cols_rows)
File "inverse_matrix.py", line 130, in compute_inverse_mat
greens_function = fac_dyn_matrix.matSolve(B,x)
File "PETSc/Mat.pyx", line 1509, in petsc4py.PETSc.Mat.matSolve
petsc4py.PETSc.Error: error code 76
[1] MatMatSolve() line 3380 in /Users/griesserj/Libaries/petsc/src/mat/interface/matrix.c
[1] MatMatSolve_MUMPS() line 1177 in /Users/griesserj/Libaries/petsc/src/mat/impls/aij/mpi/mumps/mumps.c
[1] Error in external library
[1] Error reported by MUMPS in solve phase: INFOG(1)=-3
Do you have an idea where this error can come from or did I miss to set some necessary options?
Thank you very much in advance!
Code:
def compute_inverse_mat(dynamical_matrix_nn, dyn_matrix_dim, m_cols_rows):
# Set up the matrix B
B = PETSc.Mat().create(comm=comm)
B.setSizes((dyn_matrix_dim,m_cols_rows))
B.setType("dense")
B.setFromOptions()
B.setUp()
Rstart, Rend = B.getOwnershipRange()
print("rank, size, start_frame, end_frame \n", rank, " / ", size, " / ", Rstart, " / ", Rend)
# Fill the matrix B
if (Rstart < m_cols_rows) and (Rend <= m_cols_rows):
for i in range(Rstart, Rend):
B[i,i] = 1
if (Rstart < m_cols_rows) and (Rend >= m_cols_rows):
for i in range(Rstart, m_cols_rows):
B[i,i] = 1
# Assemble B
B.assemble()
# Set up matrix x
x = PETSc.Mat().create(comm=comm)
x.setSizes((dyn_matrix_dim,m_cols_rows))
x.setType("dense")
x.setFromOptions()
x.setUp()
x.assemble()
# Create the linear solver and set various options
ksp = PETSc.KSP().create(comm=comm)
# This implements a stub method that applies ONLY the preconditioner.
ksp.setType("preonly")
# Set the matrix associated with the linear system
ksp.setOperators(dynamical_matrix_nn, dynamical_matrix_nn)
# Define the preconditioner object
pc = ksp.getPC()
# Set the preconditioner to LU-factorization
pc.setType("lu")
# Use MUMPS
pc.setFactorSolverType("mumps")
# Prepares for the use of the preconditioner
pc.setFromOptions()
pc.setUp()
# Sets up the inernal data structures for the later use of an iterative solver
ksp.setFromOptions()
ksp.setUp()
# Get the factorized dynamical matrix
fac_dyn_matrix = pc.getFactorMatrix()
# Compute part of the inverse by solving
# A*X=B
greens_function = fac_dyn_matrix.matSolve(B,x)
Am 01.10.2019 um 10:09 schrieb Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>:
On Tue, Oct 1, 2019 at 4:07 AM Jan Grießer <griesser.jan at googlemail.com<mailto:griesser.jan at googlemail.com>> wrote:
Hey Matt,
Can you elaborate a little bit on your idea for calculating the inverse matrix ?
It is exactly what you were doing before, except you use KSP with
-ksp_type preonly -pc_type lu -pc_mat_solver_package mumps
and then MatMatSolve on the identity matrix.
Thanks,
Matt
Greetings Jan
Am Mo., 30. Sept. 2019 um 17:50 Uhr schrieb Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>:
I think the easier way to do it is to use a KSP which is configured to do preonly and LU. That will do the right thing in parallel.
Matt
On Mon, Sep 30, 2019 at 11:47 AM Smith, Barry F. via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
The Python wrapper for PETSc may be missing some functionality; there is a manual process involved in creating new ones. You could poke around the petsc4py source and see how easy it would be to add more functionality that you need.
> On Sep 30, 2019, at 10:13 AM, Jan Grießer <griesser.jan at googlemail.com<mailto:griesser.jan at googlemail.com>> wrote:
>
> I configured PETSc with MUMPS and tested it already for the spectrum slicing method in Slepc4py but i have problems in setting up the LU factorization in the PETSc4py. Since i do not find the corresponding methods and commands in the source code. Thats why is was wondering if this is even possible in the python version.
>
> Am Mo., 30. Sept. 2019 um 16:57 Uhr schrieb Smith, Barry F. <bsmith at mcs.anl.gov<mailto:bsmith at mcs.anl.gov>>:
>
> If you want a parallal LU (and hence the ability to build the inverse in parallel) you need to configure PETSc with --download-mumps --download-scalapack
>
> Barry
>
>
> > On Sep 30, 2019, at 9:44 AM, Jan Grießer <griesser.jan at googlemail.com<mailto:griesser.jan at googlemail.com>> wrote:
> >
> > Is the MatMumpsGetInverse also wrapped to the python version in PETSc4py ? If yes is there any example for using it ?
> > My other question is related to the LU factoriation (https://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix).
> > Is the LU factorization only possible for sequential Aij matrices ? I read in the docs that this is the case for ordering.
> > After setting up my matrix A, B and x i tried:
> > r, c = dynamical_matrix_nn.getOrdering("nd")
> > fac_dyn_matrix = dynamical_matrix_nn.factorLU(r,c)
> >
> > resulting in an error:
> > [0] No support for this operation for this object type
> > [0] Mat type mpiaij
> >
> > Am Fr., 27. Sept. 2019 um 16:26 Uhr schrieb Zhang, Hong <hzhang at mcs.anl.gov<mailto:hzhang at mcs.anl.gov>>:
> > See ~petsc/src/mat/examples/tests/ex214.c on how to compute selected entries of inv(A) using mumps.
> > Hong
> >
> > On Fri, Sep 27, 2019 at 8:04 AM Smith, Barry F. via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
> >
> > MatMumpsGetInverse() maybe useful. Also simply using MatMatSolve() with the first 1000 columns of the identity and "throwing away" the part you don't need may be most effective.
> >
> > Barry
> >
> >
> >
> > > On Sep 27, 2019, at 3:34 AM, Jan Grießer via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
> > >
> > > Hi all,
> > > i am using petsc4py. I am dealing with rather large sparse matrices up to 600kx600k and i am interested in calculating a part of the inverse of the matrix(I know it will be a dense matrix). Due to the nature of my problem, I am only interested in approximately the first 1000 rows and 1000 columns (i.e. a large block in the upper left ofthe matrix). Before I start to play around now, I wanted to ask if there is a clever way to tackle this kind of problem in PETSc in principle. For any input I would be very grateful!
> > > Greetings Jan
> >
>
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20191008/2ac28815/attachment-0001.html>
More information about the petsc-users
mailing list