I have to solve a non-linear ODE of the shape $$\partial_zA=f(A)$$ with $f$ a non-linear function and $A$ a matrix/vector with >1e6 variables (i.e. $A$ is a matrix with >1000x1000 entries). For each calculation of the right hand side of the system (i.e. $f(A)$) I have to calculate roughly five matrix-matrix multiplications $B\cdot A$ with $B$ a constant matrix and roughly five FFTs over each row of $A$. $A$ is using complex numbers.
To solve this system I intended to use a solver such as odeint from boost. I also intended to run the calculation on a GPU, to use the increased speed of matrix-matrix-multiplications and batch FFTs. I had to find out, though, that odeint does not support thrust::complex as value type. Therefore, my current approach is to copy $A$ to the GPU in each step, and return $f(A)$ afterwards to the solver. This is rather costly, and therefore I would like to avoid it.
Thus, are there other strategies (apart from writing my own solver) which could help solving this PDE in reasonable time? Should I still try to go for the GPU, or rather try things like MPI (or a hybrid approach)? Moreover, I noticed that for several sets of parameters the equation behaves rather stiff. Are there options for solving the system in that case?
