I am porting a 2d CFD code written in Octave/Matlab to fortran. The domain is periodic, so the scheme is based on FT. The following matrix, 'laplacian':
0 -1 -4 -9 -4 -1
-1 -2 -5 -10 -5 -2
-4 -5 -8 -13 -8 -5
-9 -10 -13 -18 -13 -10
-4 -5 -8 -13 -8 -5
-1 -2 -5 -10 -5 -2
represents the Laplacian operator for a FT on a 6 by 6 grid. I want the inverse matrix, even though Laplacian is singular. In matlab/octave, 'inv(laplacian)' returns all 'Inf', however '1./laplacian' returns the correct answer (albeit the (1,1) element, returned as Inf, has to be set to zero).
The question is how to translate the second form using LAPACK. My usual matrix inversion sequence 'DGETRF/DGETRI' fails with info=4, no surprise. There are some twenty other DxxTRF. Does anyone know what might stand a chance of doing what Octave does?
If you are doing what I think you are then you want to multiply different wavenumbers of the Fourier transform by different coefficients which are derived from the eigenvalues of the Laplacian operator.
Something like
lambda(kx, ky, kz) = (kx**2 + ky**2 + kz**2)
Notice the 1, 4, and 9 in your "matrix". They are these squares kx**2.
That is NOT a matrix inversion, that is really just dividing 1.0 by numbers which are written in a form of a table. The table looks like a matrix, because your code is 2D, so you have just lambda(kx, ky).
The whole actual matrix for the Laplacian operator would be very big (N times N, where N=nx*ny in 2D and N=nx*ny*nz in 3D) and would have the lambdas on the diagonal and zeros everywhere else. The inverse matrix would have 1./labdas on the diagonal. So your operation is matrix inversion in a sense, but in a different sense than you thought.
What you than do is
FT(kx,ky,kz) = FT(kx,ky,kz) / (kx**2 + ky**2 + kz**2)
which can be also written as
FT = FT * inv_laplacian
where
inv_laplacian = 1. / laplacian
where laplacian is the coefficient (eigenvalue) array.
That is NOT a matrix inversion, it is just dividing 1.0 by numbers.
Now, what the do with the 0 eigenvalue? If you do not enable floating point exception trapping, and I advice you against enabling them, then you just do:
FT = FT * inv_laplacian
Because inv_laplacian(0,0,0) is 0, FT(0,0,0) was divided by 0 and is undefined (NaN or similar). You can just set it to whatever you want. The meaning of FT(0,0,0) is the mean value of your field, which is arbitrary. So do just
FT(0,0,0) = 0 !or any number you want
and that's it.
BTW I have seen even extreme practices in a real world scientific code such as:
for i =0, nkx-1
for j =0, nkx-1
for k =0, nkx-1
FT(i, j, k) = FT(i, j, k) / (cos(i*ax) + cos(j*ay) + cos(k*az)
end
end
end
it took AGES to compute. It is very similar to your case, your eigenvalues are just not cosines but squares.
The point is that the coefficients are CONSTANT in time and SEPARABLE.
One can just compute once:
lambdax(i) = i**2 !or cos(ax*i)
lambday(i) = j**2 !or cos(ay*j)
lambdaz(i) = k**2 !or cos(az*k)
and then do
FT(kx,ky,kz) = FT(kx,ky,kz) / (lambdax(kx) + lambday(ky) + lambdaz(kz))
You can have a look at the source code of my fast Poisson solver PoisFFT to see an example https://github.com/LadaF/PoisFFT/blob/master/src/poisfft-solvers-inc.f90 and find the appropriate boundary conditions. Your case is likely PoisFFT_Solver2D_FullPeriodic on line 152.