slepc4py with the fft

[Francis Poulin]

I am attaching a simple python code that uses a matrix free method to compute the eigenvalues of the simple differential operator:

-d^2 u/dx^2 = lambda u

First, a LinearOperator is defined that evaluates the left hand side using the FFT in scipy.

Second, the eigenvalues of this operator are computed using eigs in scipy.

I very much would like to learn how to solve this problem, and a bunch of other problems similar to it in spirit, using slepc4py.

Question: Is it possible to adapt this code to solve for the eigenvalues of the LinearOperator using slepc4py? Or how do I have to modify it in order to be able to use slepc4py?

[Jose E. Roman]

Please have a look at example ex3.py under directory demo. You can use it right away, just adapting the class Laplacian2D to use your function fd2spec instead of laplace2d. Please note that the example assumes a symmetric operator; for a general non-symmetric operator one should replace ProblemType.HEP by ProblemType.NHEP. Let us know if you have any difficulties.

[Francis Poulin]

Thank you for all the help. That does sound like it will be a small modification.

Before I try the fft, which is only 1D at this time, I thought i would simplify ex1.py to a 1d Laplacian, to make it easier to compare and build on.

I have attached the code above, ex3_1d.py, which does not work but gives me the following error. It would appear I am giving mult 4 instead of 3 arguments.

Do you have some advice as to how I might fix this?

Traceback (most recent call last): File "ex3_1d.py", line 96, in <module> main() File "ex3_1d.py", line 93, in main solve_eigensystem(A) File "ex3_1d.py", line 60, in solve_eigensystem E.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) File "libpetsc4py/libpetsc4py.pyx", line 858, in libpetsc4py.MatMult_Python (src/libpetsc4py/libpetsc4py.c:9984) TypeError: mult() takes exactly 3 arguments (4 given)

[Jose E. Roman]

You have to replace the last line of mult but do not change anything else. In particular, you have modified the number of arguments of mult, and U must still be a 2-dimensional array.

[Francis Poulin]

Ah, I think I understand. We don't want a vector but we want an (m,1) 2-dimensional array. That makes sense.

I am happy to say that I have the 1d version of ex3.py working. I will now start testing the fft and let you know when I have something working.

[Francis Poulin]

I modified my function a little bit and copied it into my modified example.

The good news is that it does run.

The problem is that now I get a k value of 0.

Do you think there is something wrong with my function?

If not, then I have to learn how to modify the slepc parameters in order to converge to the eigenvalues with the smallest magnitude, which is what is converged to in ex3.py

Thanks again for your help.

def D2fft(U, x, f):         # - d^2/dx^2

    # Grid spacing
    m, n = x.shape

    #  Define wavenumber (frequency)
    L = 1.0

    K2 = np.zeros(([m,1]))
    K2[:,0] = (2*np.pi/L*np.hstack([range(0,m/2+1), range(-m/2+1,0)]))**2 

    f = ifft(K2*fft(x)).real

[Jose E. Roman]

The original ex3.py computes largest eigenvalues.

I don't know if your function is wrong or not. I am not familiar with FFT.

[Francis Poulin]

I wanted to let you know that I did get the fft working. Thanks a lot for your help. I am uploading it here in case others might benefit from this sample.

I should point out that I made a few minor changes to the function laplacian1d from ex3.py

1. I divided by the grid spacing since I want to find the eigenvalues of the 1D Laplacian operator
2. In this Finite Difference formulation, I think that, because there are boundary values that are set to zero, the grid spacing should be hx = 1./(m+1)
3. I think that plotting the eigenvalues with the smallest magnitude is perhaps better since these are the ones that are best resolved.
4. The eigenvalues in the FFT version are largest by a factor of 4 because there we don't impose any boundary conditions but just look for periodic solution. That is why 0 is a valid eigenvalue.

I hope this helps someone.

python ex3_fft.py -eps_nev 5 -eps_smallest_magnitude

[Jose E. Roman]

Great