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 <[email protected]>:
>
> On Tue, Oct 1, 2019 at 4:07 AM Jan Grießer <[email protected]
> <mailto:[email protected]>> 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
> <[email protected] <mailto:[email protected]>>:
> 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
> <[email protected] <mailto:[email protected]>> 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 <[email protected]
> > <mailto:[email protected]>> 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.
> > <[email protected] <mailto:[email protected]>>:
> >
> > 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 <[email protected]
> > > <mailto:[email protected]>> 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
> > > <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
> > > <[email protected] <mailto:[email protected]>>:
> > > 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
> > > <[email protected] <mailto:[email protected]>> 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
> > > > <[email protected] <mailto:[email protected]>> 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/>
