How efficient is your Heston (or similar) PDE solver?

Last year I wrote an overly long post comparing the performance of an old Heston PDE solver of mine with those available in QuantLib. I was surprised to see that despite not using ADI techniques, it still seemed to produce very accurate solutions in a few seconds, just like QuantLib's ADI solvers. Published papers from the last 10 years or so also indicate that the calculation time in order to get high-accuracy prices from this 2D PDE is a few seconds give or take (on 5-10 yrs old CPUs). I've even seen a paper from last year presenting a "very efficient", state-of-the-art finite element method that gets Heston vanilla prices in 30 (!) seconds.
This year though I've worked a lot more on such solvers and I now realise that those run times are clearly not what one should be expecting, nor aim for with (fine-tuned) implementations. So just how fast should your numerical solver of the Heston (or similar) PDE be? What kind of run time should be expected for high-accuracy solutions? The answer is: a few milliseconds.  Back in that Aug 2016 post I used some test cases from various papers and I'll use those again here, all puts:
​​

	κ	η	σ	ρ	rd	rf	T	K	νo
Case 0	5	0.16	0.9	0.1	0.1	0	0.25	10	0.0625
Case 1	1.5	0.04	0.3	-0.9	0.025	0	1	100	0.0625
Case 2	3	0.12	0.04	0.6	0.01	0.04	1	100	0.09
Case 3	0.6067	0.0707	0.2928	-0.7571	0.03	0	3	100	0.0625
Case 4	2.5	0.06	0.5	-0.1	0.0507	0.0469	0.25	100	0.0625
Case 5	3	0.04	0.01	-0.7	0.05	0	0.25	100	0.09

My current "state-of-the-art" solver uses second-order spatial discretization and the Hundsdorfer-Verwer ADI scheme. For the results below I used the solver "as is", so no fine-tuning for each case, everything (grid construction) decided automatically by the solver. In order to give an idea of the accuracy achieved overall for each case I plot the solution error in the asset direction (so across the moneyness spectrum). The plots are cut-off where the option value becomes too small, the actual grids used by the solver extend further to the right. I am showing results for two different grid resolutions  (NS x NV x NT), grid A (60 x 40 x 30) and grid B (100 x 60 x 50). The timings were taken on an i7-920 PC from 2009 in single-threaded mode. Obviously one can expect at least double single-thread speed from a modern high-spec machine. (ADI can be parallelized as well with almost no extra effort, with a parallel efficiency of about 80% achieved using basic OpenMP directives). The errors in each case are calculated by comparing with the exact (semi-analytic) values obtained using QuantLib at the highest precision. Note that I am using the relative errors which are harder to bring down as the option value tends to zero. But when one wants to use the PDE solver as the pricing engine for a calibration, then it is important to price far out-of-the-money options with low relative errors in order to get accurate implied volatilities and properly capture any smile behavior. The present solver can indeed be used to calibrate the Heston model (or any other model in this family) accurately in less than a minute and in many cases in just a few seconds.

So, if you have a Heston PDE solver then this is an efficiency reference you can use to benchmark against. Price the 6 options above and compare error plots and timings with those below. Let me know what you get. If you're getting larger errors than I do with the above resolutions, I'll tell you why that is!

There is of course another important quality for a PDE solver, robustness. The scheme I used here (H-V) is fairly robust, but can still produce some spurious oscillations when one uses too low NT/NS. I may expand this post testing other schemes in the near future.