GPU support for `steadystate_fourier`

[albertomercurio]

The steadystate_fourier function doesn't support GPU arrays.

The support can be implemented by following the same method as steadystate.

One can take advantage of the following functions

• _dense_similar: Creates the equivalent dense version of the input array. Uses multiple dispatch to generate a GPU array if necessary.
• allowed_setindex!: Allows to set the index of an array. Use this only if the operation is done only once.

[aarontrowbridge]

if it makes sense, which it looks like it does, to use CUDSS.jl in LinearSolve.jl to do this I'm happy to take a wack at it this weekend

i didn't look at the work @Fe-r-oz did already to see if this approach makes sense or not

[albertomercurio]

Hi @aarontrowbridge,

As far as I know, CUDSS.jl allows to solve direct methods, right?

Here I think that we should use iterative methods, as the matrix is usually large.

If we make it work with iterative methods (which should be easier that using CUDSS.jl I think), then we can think to implement the direct solver, but maybe in another PR.

You can take the steadystate function as an example, and also the last part of this section.

[aarontrowbridge]

hey @albertomercurio sorry i didn't get back right away, i was traveling by train across Europe end of last week, then home by plane, and then to new york by train again and didn't have remote access to my desktop which has an nvidia GPU, while my mac laptop does not. i just realized today is the last day, i'm gonna see what i can do with it today but i'm going to be unavailable almost the entire day, attending an event in new york -- so it might be best to go with @Fe-r-oz's solution if it's in a good place

i do think that CUDSS.jl is a good approach to this problem, it is a direct method, but it's made for solving large, sparse systems. i also have an ulterior motive to get familiar with CUDSS.jl for my own purposes. it should also be very easy to implement, as it just requires that the passed matrix is a CuSparseMatrix format

[albertomercurio]

Hi @aarontrowbridge, no worries!

I think we can ask for an extension of the deadline up to Monday 16th. I'm ok if we implement the direct method. Just be sure to use LinearSolve.jl, so that one can uses it for iterative solvers in the future. See then steadystate case indeed.

[aarontrowbridge]

@albertomercurio that sounds great! and yes it uses the same solver interface. it might be the case that getting it working for one will be that it automatically works for the other too.