<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">Hey,<div class="">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 </div><div class="">I get the following error message:</div><div class=""><br class=""></div><div class=""><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">Traceback (most recent call last):</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> File "inverse_matrix.py", line 145, in <module></span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> compute_inverse_mat(dynamical_matrix_nn, args.dynamical_matrix_dimension, args.m_cols_rows)</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> File "inverse_matrix.py", line 130, in compute_inverse_mat</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> greens_function = fac_dyn_matrix.matSolve(B,x)</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""> File "PETSc/Mat.pyx", line 1509, in petsc4py.PETSc.Mat.matSolve</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">petsc4py.PETSc.Error: error code 76</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">[1] MatMatSolve() line 3380 in /Users/griesserj/Libaries/petsc/src/mat/interface/matrix.c</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">[1] MatMatSolve_MUMPS() line 1177 in /Users/griesserj/Libaries/petsc/src/mat/impls/aij/mpi/mumps/mumps.c</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">[1] Error in external library</span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">[1] Error reported by MUMPS in solve phase: INFOG(1)=-3</span></div></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""><br class=""></span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class="">Do you have an idea where this error can come from or did I miss to set some necessary options?</div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><br class=""></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class="">Thank you very much in advance!</div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class=""><br class=""></span></div><div style="margin: 0px; font-stretch: normal; font-size: 11px; line-height: normal; font-family: Menlo;" class=""><span style="font-variant-ligatures: no-common-ligatures" class="">Code:</span></div><div class=""><div class="">def compute_inverse_mat(dynamical_matrix_nn, dyn_matrix_dim, m_cols_rows):</div><div class=""> # Set up the matrix B </div><div class=""> B = PETSc.Mat().create(comm=comm)</div><div class=""> B.setSizes((dyn_matrix_dim,m_cols_rows))</div><div class=""> B.setType("dense")</div><div class=""> B.setFromOptions()</div><div class=""> B.setUp()</div><div class=""> Rstart, Rend = B.getOwnershipRange()</div><div class=""> print("rank, size, start_frame, end_frame \n", rank, " / ", size, " / ", Rstart, " / ", Rend)</div><div class=""> # Fill the matrix B </div><div class=""> if (Rstart < m_cols_rows) and (Rend <= m_cols_rows):</div><div class=""> <span class="Apple-tab-span" style="white-space:pre"> </span>for i in range(Rstart, Rend):</div><div class=""> <span class="Apple-tab-span" style="white-space:pre"> </span>B[i,i] = 1</div><div class=""> if (Rstart < m_cols_rows) and (Rend >= m_cols_rows):</div><div class=""> <span class="Apple-tab-span" style="white-space:pre"> </span>for i in range(Rstart, m_cols_rows):</div><div class=""> <span class="Apple-tab-span" style="white-space:pre"> </span>B[i,i] = 1</div><div class=""> # Assemble B </div><div class=""> B.assemble()</div><div class=""><br class=""></div><div class=""> # Set up matrix x</div><div class=""> x = PETSc.Mat().create(comm=comm)</div><div class=""> x.setSizes((dyn_matrix_dim,m_cols_rows))</div><div class=""> x.setType("dense")</div><div class=""> x.setFromOptions()</div><div class=""> x.setUp()</div><div class=""> x.assemble()</div><div class=""><br class=""></div><div class=""> # Create the linear solver and set various options</div><div class=""> ksp = PETSc.KSP().create(comm=comm)</div><div class=""> # This implements a stub method that applies ONLY the preconditioner.</div><div class=""> ksp.setType("preonly")</div><div class=""> # Set the matrix associated with the linear system </div><div class=""> ksp.setOperators(dynamical_matrix_nn, dynamical_matrix_nn)</div><div class=""> # Define the preconditioner object</div><div class=""> pc = ksp.getPC()</div><div class=""> # Set the preconditioner to LU-factorization</div><div class=""> pc.setType("lu")</div><div class=""> # Use MUMPS </div><div class=""> pc.setFactorSolverType("mumps")</div><div class=""> # Prepares for the use of the preconditioner </div><div class=""> pc.setFromOptions()</div><div class=""> pc.setUp()</div><div class=""> # Sets up the inernal data structures for the later use of an iterative solver </div><div class=""> ksp.setFromOptions()</div><div class=""> ksp.setUp()</div><div class=""><br class=""></div><div class=""> # Get the factorized dynamical matrix </div><div class=""> fac_dyn_matrix = pc.getFactorMatrix()</div><div class=""><br class=""></div><div class=""> # Compute part of the inverse by solving</div><div class=""> # A*X=B </div><div class=""> greens_function = fac_dyn_matrix.matSolve(B,x)</div><div><blockquote type="cite" class=""><div class="">Am 01.10.2019 um 10:09 schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>>:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" style="caret-color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none;" class=""><div dir="ltr" class="">On Tue, Oct 1, 2019 at 4:07 AM Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" class="">griesser.jan@googlemail.com</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr" class=""><div class="">Hey Matt,</div>Can you elaborate a little bit on your idea for calculating the inverse matrix ?<span class="Apple-converted-space"> </span><br class=""></div></blockquote><div class=""><br class=""></div><div class="">It is exactly what you were doing before, except you use KSP with</div><div class=""><br class=""></div><div class=""> -ksp_type preonly -pc_type lu -pc_mat_solver_package mumps</div><div class=""><br class=""></div><div class="">and then MatMatSolve on the identity matrix.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr" class="">Greetings Jan </div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Mo., 30. Sept. 2019 um 17:50 Uhr schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank" class="">knepley@gmail.com</a>>:<br class=""></div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr" class="">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.<div class=""><br class=""></div><div class=""> Matt</div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Sep 30, 2019 at 11:47 AM Smith, Barry F. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>> wrote:<br class=""></div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><br class=""> 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.<br class=""><br class=""><br class=""><br class="">> On Sep 30, 2019, at 10:13 AM, Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank" class="">griesser.jan@googlemail.com</a>> wrote:<br class="">><span class="Apple-converted-space"> </span><br class="">> 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.<span class="Apple-converted-space"> </span><br class="">><span class="Apple-converted-space"> </span><br class="">> Am Mo., 30. Sept. 2019 um 16:57 Uhr schrieb Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank" class="">bsmith@mcs.anl.gov</a>>:<br class="">><span class="Apple-converted-space"> </span><br class="">> 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<span class="Apple-converted-space"> </span><br class="">><span class="Apple-converted-space"> </span><br class="">> Barry<br class="">><span class="Apple-converted-space"> </span><br class="">><span class="Apple-converted-space"> </span><br class="">> > On Sep 30, 2019, at 9:44 AM, Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank" class="">griesser.jan@googlemail.com</a>> wrote:<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > Is the MatMumpsGetInverse also wrapped to the python version in PETSc4py ? If yes is there any example for using it ?<span class="Apple-converted-space"> </span><br class="">> > My other question is related to the LU factoriation (<a href="https://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix" rel="noreferrer" target="_blank" class="">https://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix</a>).<span class="Apple-converted-space"> </span><br class="">> > Is the LU factorization only possible for sequential Aij matrices ? I read in the docs that this is the case for ordering.<span class="Apple-converted-space"> </span><br class="">> > After setting up my matrix A, B and x i tried:<br class="">> > r, c = dynamical_matrix_nn.getOrdering("nd")<br class="">> > fac_dyn_matrix = dynamical_matrix_nn.factorLU(r,c)<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > resulting in an error:<br class="">> > [0] No support for this operation for this object type<br class="">> > [0] Mat type mpiaij<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > Am Fr., 27. Sept. 2019 um 16:26 Uhr schrieb Zhang, Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank" class="">hzhang@mcs.anl.gov</a>>:<br class="">> > See ~petsc/src/mat/examples/tests/ex214.c on how to compute selected entries of inv(A) using mumps.<br class="">> > Hong<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > On Fri, Sep 27, 2019 at 8:04 AM Smith, Barry F. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>> wrote:<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > 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.<br class="">> ><span class="Apple-converted-space"> </span><br class="">> > Barry<br class="">> ><span class="Apple-converted-space"> </span><br class="">> ><span class="Apple-converted-space"> </span><br class="">> ><span class="Apple-converted-space"> </span><br class="">> > > On Sep 27, 2019, at 3:34 AM, Jan Grießer via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>> wrote:<br class="">> > ><span class="Apple-converted-space"> </span><br class="">> > > Hi all,<br class="">> > > 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!<br class="">> > > Greetings Jan<span class="Apple-converted-space"> </span><br class="">> ><span class="Apple-converted-space"> </span><br class="">><span class="Apple-converted-space"> </span><br class=""><br class=""></blockquote></div><br clear="all" class=""><div class=""><br class=""></div>--<span class="Apple-converted-space"> </span><br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></blockquote></div></blockquote></div><br clear="all" class=""><div class=""><br class=""></div>--<span class="Apple-converted-space"> </span><br class=""><div dir="ltr" class="gmail_signature"><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a></div></div></div></div></div></div></div></div></div></blockquote></div><br class=""></div></body></html>