ex3.py: runs in serial but not in parallel

[Francis Poulin]

I am very new to petsc/slepc, python and petsc4py/slepc4py and am trying to learn about each of them but admit I still have a ways to go. I have a couple of questions I'd like to ask.

1) First is about ex3.py that comes as a demo with slepc4py. I run it in serial and it works fine. But when I try and run it in parallel I get the following error. I have not changed it and it is the only example that doesn't seem to run. Is there a problem with how my system is set up? If you have any suggestions as to what I should try and I am happy to give it a go.

All of the other examples seem to run in serial and in parallel but when I try fpoulin@vortex:~/software/slepc4py-3.3.1/demo$ mpiexec -np 2 python ex3.py Symmetric Eigenproblem (matrix-free), N=1024 (32x32 grid) Traceback (most recent call last): File "ex3.py", line 101, in <module> Traceback (most recent call last): File "ex3.py", line 101, in <module> main() File "ex3.py", line 98, in main main() File "ex3.py", line 98, in main solve_eigensystem(A) File "ex3.py", line 64, in solve_eigensystem solve_eigensystem(A) File "ex3.py", line 64, in solve_eigensystem E.solve() E.solve() File "EPS.pyx", line 1039, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:16594) File "EPS.pyx", line 1039, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:16594) File "libpetsc4py.pyx", line 834, in libpetsc4py.MatMult_Python (src/libpetsc4py/libpetsc4py.c:8664) File "libpetsc4py.pyx", line 834, in libpetsc4py.MatMult_Python (src/libpetsc4py/libpetsc4py.c:8664) File "ex3.py", line 37, in mult File "ex3.py", line 37, in mult xx = x[...].reshape(m,n) xx = x[...].reshape(m,n) ValueError: ValueError: total size of new array must be unchanged total size of new array must be unchanged

2) This is more of a slepc question and maybe should make this separate. If so please let me know and I will send this elsewhere. I want to understand ex3.py because I want to compute the eigenvalues of a 2D system. The example is for a regular eigenvalue problem but one problem I want to solve is for a generalized eigenvalue problem: A x = lambda B x.

From what I've read it seems that matrix free methods are what are recommended, which is why I want to figure out this example. I understand how to set up the A part of it and suspect that I need to solve a system to essentially isolate the lambda*x. If there are any examples of how to do this in slepc4py or slepc I would be happy to look at those in details.

Thank in advance and sorry for the bother.

Best regards, Francis

[Jose E. Roman]

1) This example is not intended to run in parallel, it should fail gracefully. It is the equivalent of SLEPc's ex3.c, which checks that the number of processes is 1. If you want to run an example that uses a shell matrix in parallel then the shell matrix must support parallelism, that is, it must implement all communication required to do the matrix-vector multiplication in parallel.

2) Shell matrices are good for problems that can be solved only using matrix-vector products. That is the case for standard eigenvalue problems A x = lambda x, but for generalized problems A x = lambda B x the situation is different because linear systems must be solved within the eigensolver iteration, and this complicates things a lot.

If possible, it is better to just represent the matrix explicitly, computing its coefficients as in ex2.py

[Francis Poulin]

Hello Jose,

Thank you very much for your responses.

1) I had actually looked over the code to try and understand it and I should have appreciated it was only serial. My apologies for my carelessness. I will try and be more careful in the future before I start asking questions.

2) Actually, my first thought was to use sparse matrices but I read some news groups and thought that shell matrices might be the better way to go. If you think that using sparse matrices is easier than I will start taking that route. But out of curiosity, are there any examples floating around for how to solve generalized eigenvalue problems using shell matrices?

3) Thanks to 2) I now have a small additional question. Because I am wanting to discretize a Partial Differential Equation in two-dimensions I need to use a kronecker product. I could write my own if need be but in hopes of saving some time is there anything in slepc4py? In python there is numpy.kron that should do exactly what I need. If not I am prepared to write my own but thought I would check since creating this in serial should be easy but in parallel sounds a little tricky.

Best regards, Francis

[Jose E. Roman]

There is no example of generalized eigenproblem using shell matrices. In slepc4py there is nothing for Kronecker products but you can use numpy.kron if it helps in computing the coefficients to be stored locally by each process.

[Francis Poulin]

Thank you.

I just found a link that talked about how to convert arrays from numpy to petsc. So I think I will build the sparse matrix in numpy and then convert it to AIJ type.

http://lists.mcs.anl.gov/pipermail/petsc-users/2012-May/013379.html

Thanks again for the help.

[Francis Poulin]

We are trying to translate ex3.py into parallel and we are having difficulties. It seems that we should define a vector, share it between the codes and define ghost cells. We are having some problems in figuring out how to do this. Could someone perhaps help us to get started?

Thanks in advance.

[Jose E. Roman]

See poisson2d.py in petsc4py

[Francis Poulin]

Thank you.

We are looking at the code for poisson2d.py to figure out how to parallelize the example. However, before we implement that, we are playing with the matrix-free version of our code and having some problems.

For the matrix version, we can use shift and invert to find the eigenvalue in a particular range (usually with a non-zero imaginary part) but when we try applying the same procedure to the matrix-free method we get an error saying that the matrix cannot be duplicated.

File "EPS.pyx", line 1039, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:16594) petsc4py.PETSc.Error: error code 83 [0] EPSSolve() line 90 in src/eps/interface/solve.c [0] EPSSetUp() line 215 in src/eps/interface/setup.c [0] STSetUp() line 294 in src/st/interface/stsolve.c [0] STSetUp_Sinvert() line 138 in src/st/impls/sinvert/sinvert.c [0] STMatGAXPY_Private() line 363 in src/st/interface/stsolve.c [0] MatDuplicate() line 4063 in src/mat/interface/matrix.c [0] MatDuplicate_Python() line 2510 in petsc4py-3.4/src/libpetsc4py/libpetsc4py.c [0] [0] method duplicate()

I presume that means we cannot use shift-invert in the matrix-free approach? If yes, is there another technique that we should consider to find particular eigenvalues (possibly complex) with a matrix-free approach?

[Jose E. Roman]

Read section 3.4.2 of SLEPc's users guide.